text:
using BilevelJuMP
using Test, MathOptInterface, JuMP, Dualization
# using MathOptFormat

const MOI  = MathOptInterface
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges
const MOIT = MathOptInterface.Test

# TODO
# add JUMPExtension test

struct Config
    atol::Float64
    rtol::Float64
    bound_hint::Bool
    start_value::Bool
    function Config(;
        atol = 1e-6,
        rtol = 1e-6,
        bound_hint = false,
        start_value = false,
        )
        new(atol, rtol, bound_hint, start_value)
    end
end

config = Config()
CONFIG_1 = Config(atol = 1e-1, rtol = 1e-2)
CONFIG_2 = Config(atol = 1e-2, rtol = 1e-2)
CONFIG_3 = Config(atol = 1e-3, rtol = 1e-3)
CONFIG_3_start = Config(atol = 1e-3, rtol = 1e-3, start_value = true)
CONFIG_3_hint = Config(atol = 1e-3, rtol = 1e-3, bound_hint = true)
CONFIG_4 = Config(atol = 1e-4, rtol = 1e-4)
CONFIG_5 = Config(atol = 1e-5, rtol = 1e-5)

OptModeType = NamedTuple{(:opt, :mode),Tuple{Any,Any}}

solvers = OptModeType[]
solvers_cached = OptModeType[]
solvers_sos = OptModeType[]
solvers_unit = OptModeType[]
solvers_indicator = OptModeType[]
solvers_quad = OptModeType[]
solvers_bin_exp = OptModeType[]
solvers_sos_quad = OptModeType[]
solvers_nlp = OptModeType[]
solvers_nlp_sd = OptModeType[]
solvers_nlp_sd_e = OptModeType[]
solvers_nlp_sd_i = OptModeType[]
solvers_nlp_lowtol = OptModeType[]
solvers_sos_quad_bin = OptModeType[]
solvers_fa_quad_bin = OptModeType[]
solvers_fa_quad_bin_mixed = OptModeType[]
solvers_fa = OptModeType[]
solvers_fa2 = OptModeType[]
solvers_complements = OptModeType[]

include("solvers/ipopt.jl")
include("solvers/scip.jl")
# DONE
# include("solvers/cbc.jl")
# include("solvers/gurobi.jl")
# include("solvers/xpress.jl")
# include("solvers/knitro.jl")
# include("solvers/gams.jl")
# include("solvers/couenne.jl")
# include("solvers/bonmin.jl")
# include("solvers/baron.jl")
# TODO
# include("solvers/alpine.jl") # require NLP from JuMP
# include("solvers/pavito.jl") # require NLP from JuMP
# include("solvers/juniper.jl") # require NLP from JuMP
# include("solvers/path.jl") # require LCP method

include("moi.jl")
include("jump.jl")
include("jump_unit.jl")
include("mibs.jl")
include("jump_nlp.jl")

@testset "BilevelJuMP tests" begin

@testset "nlp" begin
    jump_nlp_01(IPO_OPT, mode = BilevelJuMP.ProductMode(1e-8))
    jump_nlp_02(IPO_OPT, mode = BilevelJuMP.ProductMode(1e-8))
    jump_nlp_03(IPO_OPT, mode = BilevelJuMP.ProductMode(1e-8))
    jump_nlp_04(IPO_OPT, mode = BilevelJuMP.ProductMode(1e-8))
end

@testset "Unit" begin
    jump_display()
    jump_objective()
    jump_bounds()
    jump_attributes()
    mixed_mode_unit()
    jump_constraints()
    jump_variables()
    variables_unit()
    jump_no_cb()
    constraint_unit()
    for solver in solvers_unit
        invalid_lower_objective(solver.opt, solver.mode)
        jump_display_solver(solver.opt, solver.mode)
        invalid_optimizer(solver.opt, solver.mode)
        jump_objective_solver(solver.opt, solver.mode)
        jump_attributes_solver(solver.opt, solver.mode)
    end
end

@testset "Simple LP" begin
    for solver in solvers_cached
        moi_01(solver.opt)
        moi_02(solver.opt, solver.mode)
        moi_03(solver.opt, solver.mode)
    end
end

@testset "Simple BLP JuMP" begin
    for solver in solvers_nlp
        jump_01(solver.opt, solver.mode, CONFIG_3)
        jump_01vec(solver.opt, solver.mode, CONFIG_3)
        jump_02(solver.opt, solver.mode) # numerical isntability in ipopt
        jump_03(solver.opt, solver.mode, CONFIG_3_start)
        jump_03_vec(solver.opt, solver.mode, CONFIG_3_start)
        jump_04(solver.opt, solver.mode, CONFIG_3_start)
        jump_05(solver.opt, solver.mode)
        jump_3SAT(solver.opt, solver.mode, CONFIG_3)
        jump_06(solver.opt, solver.mode)
        jump_06_sv(solver.opt, solver.mode)
        jump_07(solver.opt, solver.mode, CONFIG_2)
        jump_08(solver.opt, solver.mode, CONFIG_3_start)
        jump_09a(solver.opt, solver.mode)
        jump_09b(solver.opt, solver.mode)
        jump_11a(solver.opt, solver.mode)
        jump_11b(solver.opt, solver.mode)
        jump_12(solver.opt, solver.mode)
        jump_14(solver.opt, solver.mode)
    end
    for solver in solvers_sos
        jump_01_mixed(solver.opt)
    end
    for solver in solvers_sos
        jump_01(solver.opt, solver.mode)#
        jump_02(solver.opt, solver.mode)#
        jump_03(solver.opt, solver.mode)
        jump_04(solver.opt, solver.mode)
        jump_05(solver.opt, solver.mode)#
        jump_3SAT(solver.opt, solver.mode)
        jump_06(solver.opt, solver.mode)#
        jump_06_sv(solver.opt, solver.mode)
        jump_07(solver.opt, solver.mode)#
        jump_08(solver.opt, solver.mode)#
        jump_09a(solver.opt, solver.mode) # fail on cbc positive SOS
        jump_09b(solver.opt, solver.mode)
        jump_11a(solver.opt, solver.mode)#
        jump_11b(solver.opt, solver.mode)
        jump_12(solver.opt, solver.mode)#
        jump_14(solver.opt, solver.mode)
        #
        jump_16(solver.opt, solver.mode)
    end
    for solver in solvers_bin_exp
        jump_01(solver.opt, solver.mode, CONFIG_3_hint)
        jump_01vec(solver.opt, solver.mode, CONFIG_3_hint)
        jump_02(solver.opt, solver.mode, CONFIG_3_hint)
        jump_03(solver.opt, solver.mode, CONFIG_3_hint)
    end
    for solver in solvers_fa
        jump_01(solver.opt, solver.mode, CONFIG_3_hint)
        jump_01vec(solver.opt, solver.mode, CONFIG_3_hint)
        jump_02(solver.opt, solver.mode, CONFIG_3_hint)
        jump_03(solver.opt, solver.mode, CONFIG_3_hint)
    end
    for solver in solvers_fa2
        jump_01(solver.opt, solver.mode, CONFIG_3)
        jump_01vec(solver.opt, solver.mode, CONFIG_3)
        jump_02(solver.opt, solver.mode, CONFIG_3)
        jump_03(solver.opt, solver.mode, CONFIG_3)
    end
    for solver in solvers_indicator
        jump_01(solver.opt, solver.mode)
        jump_01vec(solver.opt, solver.mode)
        jump_02(solver.opt, solver.mode) # fail cbc - pass xpress 8.9
        jump_03(solver.opt, solver.mode)
    end
    for solver in solvers_complements # requires KNITRO
        # jump_01(solver.opt, solver.mode, CONFIG_3)
        # jump_01vec(solver.opt, solver.mode, CONFIG_3)
        # jump_02(solver.opt, solver.mode) # numerical isntability in ipopt
        jump_03(solver.opt, solver.mode, CONFIG_3_start)
        jump_03_vec(solver.opt, solver.mode, CONFIG_3_start)
        jump_04(solver.opt, solver.mode, CONFIG_3_start)
        # jump_05(solver.opt, solver.mode)
        jump_3SAT(solver.opt, solver.mode)
        jump_06(solver.opt, solver.mode)
        jump_06_sv(solver.opt, solver.mode)
        # jump_07(solver.opt, solver.mode, CONFIG_2)
        jump_08(solver.opt, solver.mode, CONFIG_3_start)
        jump_09a(solver.opt, solver.mode)
        # jump_09b(solver.opt, solver.mode)
        jump_11a(solver.opt, solver.mode)
        jump_11b(solver.opt, solver.mode)
        # jump_12(solver.opt, solver.mode)
        # jump_14(solver.opt, solver.mode)
    end
end

@testset "JuMP quad" begin
    for solver in solvers_quad
        jump_quad_01_a(solver.opt, solver.mode)
        jump_quad_01_b(solver.opt, solver.mode)
        jump_quad_01_c(solver.opt, solver.mode)
        jump_13_quad(solver.opt, solver.mode)
    end
end

@testset "Princeton Handbook Linear" begin
    for solver in vcat(solvers_nlp, solvers_nlp_sd)
        jump_HTP_lin01(solver.opt, solver.mode, CONFIG_3_start)
        # jump_HTP_lin02(solver.opt, solver.mode, CONFIG_4)
        jump_HTP_lin03(solver.opt, solver.mode)
        jump_HTP_lin03_vec(solver.opt, solver.mode)
        jump_HTP_lin04(solver.opt, solver.mode)
        jump_HTP_lin05(solver.opt, solver.mode) # broken on cbc linux on julia 1.0 and 1.2 but not 1.1 see: https://travis-ci.org/joaquimg/BilevelJuMP.jl/builds/619335351
        jump_HTP_lin06(solver.opt, solver.mode)
        jump_HTP_lin07(solver.opt, solver.mode, CONFIG_2)
        jump_HTP_lin09(solver.opt, solver.mode)
        # jump_HTP_lin10(solver.opt, solver.mode)
    end
    for solver in solvers_nlp
        jump_HTP_lin02(solver.opt, solver.mode)
        jump_HTP_lin10(solver.opt, solver.mode)
    end
    for solver in solvers_nlp_sd
        jump_HTP_lin02(solver.opt, solver.mode, CONFIG_2)
    end
    for solver in solvers_nlp_sd_i
        # jump_HTP_lin08(solver.opt, solver.mode, CONFIG_1)
        jump_HTP_lin10(solver.opt, solver.mode, CONFIG_4)
    end
    for solver in solvers_nlp_sd_e
        # TODO add dual start
        # jump_HTP_lin08(solver.opt, solver.mode, CONFIG_4)
        # jump_HTP_lin10(solver.opt, solver.mode)
    end
    for solver in solvers_sos
        jump_HTP_lin01(solver.opt, solver.mode)
        jump_HTP_lin02(solver.opt, solver.mode)
        jump_HTP_lin03(solver.opt, solver.mode) # failing cbc
        jump_HTP_lin04(solver.opt, solver.mode)
        jump_HTP_lin05(solver.opt, solver.mode) # broken on cbc linux on julia 1.0 and 1.2 but not 1.1 see: https://travis-ci.org/joaquimg/BilevelJuMP.jl/builds/619335351
        jump_HTP_lin06(solver.opt, solver.mode)
        jump_HTP_lin07(solver.opt, solver.mode)
        jump_HTP_lin08(solver.opt, solver.mode)
        jump_HTP_lin09(solver.opt, solver.mode)
        jump_HTP_lin10(solver.opt, solver.mode)
    end
    for solver in solvers_complements # requires KNITRO
        jump_HTP_lin01(solver.opt, solver.mode, CONFIG_3_start)
        jump_HTP_lin02(solver.opt, solver.mode, CONFIG_4)
        # jump_HTP_lin03(solver.opt, solver.mode)
        # jump_HTP_lin03_vec(solver.opt, solver.mode)
        # jump_HTP_lin04(solver.opt, solver.mode)
        jump_HTP_lin05(solver.opt, solver.mode) # broken on cbc linux on julia 1.0 and 1.2 but not 1.1 see: https://travis-ci.org/joaquimg/BilevelJuMP.jl/builds/619335351
        jump_HTP_lin06(solver.opt, solver.mode)
        # jump_HTP_lin07(solver.opt, solver.mode, CONFIG_2)
        # jump_HTP_lin09(solver.opt, solver.mode)
        jump_HTP_lin10(solver.opt, solver.mode)
    end
end

@testset "Princeton Handbook Quadratic" begin
    for solver in solvers_nlp
        jump_HTP_quad01(solver.opt, solver.mode)
        jump_HTP_quad02(solver.opt, solver.mode)
        jump_HTP_quad04(solver.opt, solver.mode, CONFIG_3)
        jump_HTP_quad05(solver.opt, solver.mode)
        jump_HTP_quad06(solver.opt, solver.mode, CONFIG_3)
        # jump_HTP_quad06b(solver.opt, solver.mode) # TODO weird results
        jump_HTP_quad07(solver.opt, solver.mode)
        jump_HTP_quad08(solver.opt, solver.mode) # not PSD
    end
    for solver in solvers_nlp
        jump_HTP_quad03(solver.opt, solver.mode)
        jump_HTP_quad09(solver.opt, solver.mode)
    end
    for solver in solvers_sos_quad
        jump_HTP_quad01(solver.opt, solver.mode, CONFIG_5)
        jump_HTP_quad02(solver.opt, solver.mode)
        jump_HTP_quad04(solver.opt, solver.mode)
        jump_HTP_quad05(solver.opt, solver.mode)
        jump_HTP_quad06(solver.opt, solver.mode)
        jump_HTP_quad06b(solver.opt, solver.mode)
        jump_HTP_quad07(solver.opt, solver.mode)
        jump_HTP_quad08(solver.opt, solver.mode) # not PSD
    end
    for solver in solvers_sos
        jump_HTP_quad03(solver.opt, solver.mode)
        jump_HTP_quad09(solver.opt, solver.mode)
    end
end

@testset "Ferris educational" begin
    for solver in solvers
        jump_jointc1(solver.opt, solver.mode)
        jump_jointc2(solver.opt, solver.mode)
    end
end

@testset "Bard Efficient Point Example" begin
    for solver in solvers_sos
        jump_EffPointAlgo(solver.opt, solver.mode)
    end
end

@testset "SemiRandom" begin
    for solver in solvers_sos
        jump_SemiRand(solver.opt, solver.mode)
    end
end

@testset "DNE" begin
    for solver in solvers
        jump_DTMP_01(solver.opt, solver.mode)
    end
end

@testset "DNE - modifications" begin
    for solver in solvers
        jump_DTMP_01_mod1(solver.opt, solver.mode)
        jump_DTMP_01_mod2_error(solver.opt, solver.mode)
        jump_DTMP_01_mod3vec(solver.opt, solver.mode)
        jump_DTMP_01_mod3(solver.opt, solver.mode)
    end
end

@testset "equilibrium" begin
    for solver in solvers_nlp
        jump_conejo2016(solver.opt, solver.mode)
        jump_fanzeres2017(solver.opt, solver.mode)
        # jump_eq_price(solver.opt, solver.mode)
    end
    for solver in solvers_nlp_sd
        jump_conejo2016(solver.opt, solver.mode)
        # jump_fanzeres2017(solver.opt, solver.mode)
        # jump_eq_price(solver.opt, solver.mode)
    end
    for solver in solvers_sos_quad_bin
        jump_conejo2016(solver.opt, solver.mode, config, bounds = true) # fail travis on cbc
        # jump_fanzeres2017(solver.opt, solver.mode)
        jump_eq_price(solver.opt, solver.mode) # fail travis on cbc
    end
    for solver in solvers_fa_quad_bin
        jump_conejo2016(solver.opt, solver.mode, config, bounds = true)
        # jump_fanzeres2017(solver.opt, solver.mode)
        jump_eq_price(solver.opt, solver.mode)
    end
end

@testset "Fruits" begin
    for solver in solvers_fa2
        jump_fruits(solver.opt, solver.mode, CONFIG_4, 0.05)
        jump_fruits(solver.opt, solver.mode, CONFIG_4, 0.09)
    end
end

@testset "Bilevel Conic JuMP NLP" begin
    for solver in solvers_nlp_lowtol
        jump_conic01(solver.opt, solver.mode)
        jump_conic02(solver.opt, solver.mode)
        jump_conic03(solver.opt, solver.mode)
        jump_conic04(solver.opt, solver.mode)
    end
end

@testset "Bilevel Conic JuMP MIP" begin
    for solver in solvers_fa_quad_bin_mixed
        # @time jump_conic01(solver.opt, solver.mode, config, bounds = true)
        @time jump_conic02(solver.opt, solver.mode, config, bounds = true)
        @time jump_conic03(solver.opt, solver.mode, config, bounds = true)
        @time jump_conic04(solver.opt, solver.mode, config, bounds = true)
    end
end

end


text:
using  Random, Test, WaveFD, LinearAlgebra

nthreads = 12
nz,ny,nx = 101, 200, 301
dz,dy,dx = 25.0f0,25.0f0, 25.0f0
nsponge = 40
freqQ = 5.0f0
qMin = 0.1f0
qInterior = 100.0f0
dtmod = 0.002f0

@testset "2D, serial, freesurface=$(freesurface)" for freesurface in (false,true)
    w1 = zeros(Float32,nz,nx)
    w2 = zeros(Float32,nz,nx)
    WaveFD.setup_q_profile_2D_serial!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    w2 .= WaveFD.DtOmegaInvQ(Prop2DAcoIsoDenQ_DEO2_FDTD(;
        freesurface=freesurface, nz=nz, nx=nx, dz=dz, dx=dx, nsponge=nsponge, nthreads=nthreads,
        dt=dtmod, freqQ=freqQ, qMin=qMin, qInterior=qInterior))
    @test w1 ≈ w2
end

@testset "2D, threaded, freesurface=$(freesurface)" for freesurface in (false,true)
    w1 = zeros(Float32,nz,nx)
    w2 = zeros(Float32,nz,nx)
    WaveFD.setup_q_profile_2D_threaded!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    w2 .= WaveFD.DtOmegaInvQ(Prop2DAcoIsoDenQ_DEO2_FDTD(;
        freesurface=freesurface, nz=nz, nx=nx, dz=dz, dx=dx, nsponge=nsponge, nthreads=nthreads,
        dt=dtmod, freqQ=freqQ, qMin=qMin, qInterior=qInterior))
    @test w1 ≈ w2
end

@testset "2D, serial == threaded, freesurface=$(freesurface)" for freesurface in (false,true)
    w1 = zeros(Float32,nz,nx)
    w2 = zeros(Float32,nz,nx)
    WaveFD.setup_q_profile_2D_serial!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    WaveFD.setup_q_profile_2D_threaded!(w2, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    @test w1 ≈ w2
end

@testset "3D, serial, freesurface=$(freesurface)" for freesurface in (false,true)
    w1 = zeros(Float32,nz,ny,nx)
    w2 = zeros(Float32,nz,ny,nx)
    WaveFD.setup_q_profile_3D_serial!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    w2 .= WaveFD.DtOmegaInvQ(Prop3DAcoIsoDenQ_DEO2_FDTD(;
        freesurface=freesurface, nz=nz, ny=ny, nx=nx, dz=dz, dy=dy, dx=dx, nsponge=nsponge, nthreads=nthreads,
        dt=dtmod, freqQ=freqQ, qMin=qMin, qInterior=qInterior))
    @test w1 ≈ w2
end

@testset "3D, threaded, freesurface=$(freesurface)" for freesurface in (false,true)
    w1 = zeros(Float32,nz,ny,nx)
    w2 = zeros(Float32,nz,ny,nx)
    WaveFD.setup_q_profile_3D_threaded!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    w2 .= WaveFD.DtOmegaInvQ(Prop3DAcoIsoDenQ_DEO2_FDTD(;
        freesurface=freesurface, nz=nz, ny=ny, nx=nx, dz=dz, dy=dy, dx=dx, nsponge=nsponge, nthreads=nthreads,
        dt=dtmod, freqQ=freqQ, qMin=qMin, qInterior=qInterior))
    @test w1 ≈ w2
end

@testset "3D, serial == threaded, freesurface=$(freesurface)" for freesurface in (false,true)
    w1 = zeros(Float32,nz,ny,nx)
    w2 = zeros(Float32,nz,ny,nx)
    WaveFD.setup_q_profile_3D_serial!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    WaveFD.setup_q_profile_3D_threaded!(w2, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    @test w1 ≈ w2
end

@testset "2D, unphysical q values result in array set to zero, freesurface=$(freesurface)" for freesurface in (false,true)
    qMin = 0f0
    qInterior = 0f0
    w1 = ones(Float32,nz,nx)
    w2 = ones(Float32,nz,nx)
    w3 = ones(Float32,nz,nx)
    WaveFD.setup_q_profile_2D_serial!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    WaveFD.setup_q_profile_2D_threaded!(w2, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    w3 .= WaveFD.DtOmegaInvQ(Prop2DAcoIsoDenQ_DEO2_FDTD(;
        freesurface=freesurface, nz=nz, nx=nx, dz=dz, dx=dx, nsponge=nsponge, nthreads=nthreads,
        dt=dtmod, freqQ=freqQ, qMin=qMin, qInterior=qInterior))
    for w in (w1, w2, w3)
        @test w ≈ zeros(Float32,nz,nx)
    end
end

@testset "3D, unphysical q values result in array set to zero, freesurface=$(freesurface)" for freesurface in (false,true)
    qMin = 0f0
    qInterior = 0f0
    w1 = ones(Float32,nz,ny,nx)
    w2 = ones(Float32,nz,ny,nx)
    w3 = ones(Float32,nz,ny,nx)
    WaveFD.setup_q_profile_3D_serial!(w1, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    WaveFD.setup_q_profile_3D_threaded!(w2, freesurface, nsponge, dtmod, freqQ, qMin, qInterior)
    w3 .= WaveFD.DtOmegaInvQ(Prop3DAcoIsoDenQ_DEO2_FDTD(;
            freesurface=freesurface, nz=nz, ny=ny, nx=nx, dz=dz, dy=dy, dx=dx, nsponge=nsponge, nthreads=nthreads,
            dt=dtmod, freqQ=freqQ, qMin=qMin, qInterior=qInterior))
    for w in (w1, w2, w3)
        @test w ≈ zeros(Float32,nz,ny,nx)
    end
end

nothing

ext: jl

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
function complex_test()
    c1 = 2 + 3im
    c2 = 4 + 5im
    c3 = c1 + c2
    @assert(c3 == (6 + 8im))
    c4 = c1 + 3
    @assert(c4 == (5 + 3im))
    c5 = c1 + 4im
    @assert(c5 == (2 + 7im))
    c6 = c3 - 2.3
    @assert(c6 == (3.7 + 8im))
end

if abspath(PROGRAM_FILE) == @__FILE__
    complex_test()
    println("OK")
end


text:
include("required.jl")

## Conditional arguments
# outcome = "raw"
# outcome = "bart"
# save_figures = true

# Make sure "outcome" is defined
if !(@isdefined outcome) || !(outcome in ["raw", "bart"])
    @error "'outcome' must be either 'raw' or 'bart'"
else
    @info "'outcome' currently set to '$(outcome)'"
end

## Load distribution data for all species

# Load species names
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame)
spenames = filter(:type => ==("species"), glossary).full_name

# Load distributions
distributions = [
    geotiff(
        SimpleSDMPredictor, joinpath("data", "raster", "distributions_$(outcome).tif"), i
    ) for i in eachindex(spenames)
]

## Extract subareas

# Northeast subarea
coords_NE = (left=-80.0, right=-60.0, bottom=40.0, top=50.0)
distributions_NE = [d[coords_NE] for d in distributions]

# Southwest subarea
coords_SW = (left=-120.0, right=-100.0, bottom=30.0, top=40.0)
distributions_SW = [d[coords_SW] for d in distributions]

# Get Ymatrices
Y_NE = Ymatrix(distributions_NE)
Y_SW = Ymatrix(distributions_SW)

# Get Y matrices for observed sites only
Yobs_NE = Ymatrix(distributions_NE; observed=true)
Yobs_SW = Ymatrix(distributions_SW; observed=true)

# Richness
richness_NE = richness(Y_NE, distributions_NE[1])
richness_SW = richness(Y_SW, distributions_SW[1])

# Absolute LCBD values
lcbd_NE = lcbd(Y_NE, distributions_NE[1]; relative=false)
lcbd_SW = lcbd(Y_SW, distributions_SW[1]; relative=false)

# Total beta
beta_NE = beta_total(Y_NE)
beta_SW = beta_total(Y_SW)

## Subarea figures

# Choose plotting function according to outcome
if outcome == "raw"
    plotfct = :plot_layer # slower, but background is nice for sparse data
else
    plotfct = :plot # faster, sufficient for continuous data
end

# Functions to add total beta diversity in rectangle box
rectangle(w, h, x, y) = Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])
function rectangle!(w, h, x, y, textstring, textsize)
    plot!(
        rectangle(w, h, x, y),
        color=:white, legend=:none,
        annotations=(x + (w / 2), y + (h / 2), text(textstring, textsize, :center)),
    )
end

# Function to produce combined plots
function plot_lcbd_relationship(
    richness_layer, lcbd_layer, beta_total; scale=true, scaling_value=1, kw...
)
    # Scale values for a nicer result
    if scale
        # Set scaling value automatically (optional, set scaling_value=... otherwise)
        if scaling_value == 1
            scaling_factor = lcbd_layer |> maximum |> log10 |> abs |> ceil |> Int
            scaling_value = 10^scaling_factor
        end
        lcbd_layer = rescale(lcbd_layer, extrema(lcbd_layer) .* scaling_value)
    end

    # LCBD map subpanel
    p1 = eval(plotfct)(
        lcbd_layer;
        c=:viridis,
        colorbar_title="LCBD value (x $(format(scaling_value, commas=true)))",
        clim=extrema(lcbd_layer)
    )

    # Relationship subpanel
    p2 = histogram2d(
        richness_layer,
        lcbd_layer;
        c=:cividis,
        bins=40,
        colorbar_title="Number of sites",
        xlabel="Richness",
        ylabel="LCBD value (x $(format(scaling_value, commas=true)))",
        xlim=(1, 50),
        ylim=(-Inf, Inf),
    )
    vline!([median(richness_layer)], label=:none, linestyle=:dash, c=:grey)
    hline!([median(lcbd_layer)], label=:none, linestyle=:dash, c=:grey)

    # Add total beta value in rectangle box (fitted to dimensions)
    lmin, lmax = extrema(lcbd_layer)
    lrange = lmax - lmin
    rectangle!(
        16.0,
        0.15 * lrange,
        33.0,
        lmax - 0.2 * lrange,
        "BDtot = $(round(beta_total; digits=3))",
        7,
    )

    # Combine in subpanels in single plot
    l = @layout [a{0.5w} b{0.5w,0.9h}]
    p = plot(
        p1,
        p2;
        layout=l,
        size=(900, 300),
        # bottommargin = [0.0mm 5.0mm],
        rightmargin=[5.0mm 0mm],
        leftmargin=[5.0mm 5.0mm],
        kw...
    )
    return p
end

# Create subarea figures
resNEtr = plot_lcbd_relationship(richness_NE, lcbd_NE, beta_NE)
resSWtr = plot_lcbd_relationship(richness_SW, lcbd_SW, beta_SW; scaling_value=1000)

# Fix axis dimensions
yticks!(resNEtr[1], 40:2:50)
yticks!(resSWtr[1], 30:2:40)
xticks!(resNEtr[1], -80:4:-60)
xticks!(resSWtr[1], -120:4:-100)
if outcome == "bart"
    local cmax = 450
    resNEtr[2][:clims] = (-Inf, cmax)
    resSWtr[2][:clims] = (-Inf, cmax)
elseif outcome == "raw"
    local cmax = 285
    resNEtr[2][:clims] = (-Inf, cmax)
    resSWtr[2][:clims] = (-Inf, cmax)
end

# Combine subarea figures in 4 panel figure
combined_plot = plot(
    resNEtr,
    resSWtr;
    layout=grid(2, 1),
    size=(900, 600),
    bottommargin=0.0mm,
    title=["a) Northeast subregion" "" "b) Southwest subregion" ""],
    dpi=600,
)

# Export figures
# save_figures = true
if (@isdefined save_figures) && save_figures == true
    savefig(combined_plot, joinpath("fig", outcome, "05_$(outcome)_subareas.png"))
end

## Repeat for different subareas

function plot_subareas(
    coords,
    initial_distributions;
    display_coords=coords,
    transform=true,
    kw...
)
    # Get analysis values
    ds = [d[coords] for d in initial_distributions]
    y = Ymatrix(ds)
    r = richness(y, ds[1])
    l = lcbd(y, ds[1]; transform=transform)
    bt = beta_total(y)

    # Plot subareas
    p = plot_lcbd_relationship(r, l, bt; kw...)

    return p
end

# Set coordinates for initial subarea
left = -71.0; right = -64.0; bottom = 46.0; top = 50.0;
coords_subarea = (left=left, right=right, bottom=bottom, top=top)

# Test plot initial subarea
p = plot_subareas(
    coords_subarea,
    distributions;
    formatter=[f -> "$(Int(f))" :plain],
    title=["LCBD" "Relationship"]
)

## Create expanding GIF

# Set initial coordinates
left = -71.0; right = -64.0; bottom = 46.0; top = 50.0;
coords_subarea = (left=left, right=right, bottom=bottom, top=top)

# Set values to control spatial expansion
asp_ratio = 92.60 / 60.75
dim_ratio = (top - bottom)/(right - left)

# Get increasing subareas coordinates
subarea_coords = []
@time while left > -145.0 + asp_ratio && bottom > 20.0 + asp_ratio * dim_ratio
    global left -= asp_ratio
    global bottom -= asp_ratio * dim_ratio
    bbox = (left=left, right=right, bottom=bottom, top=top)
    push!(subarea_coords, bbox)
end

# Plot all expanding subareas
subarea_plots = []
for sc in subarea_coords
    local p = plot_subareas(
        sc,
        distributions;
        formatter=[f -> "$(Int(round(f, digits=0)))" :plain],
        title=["LCBD" "Relationship"],
        dpi=200
    )
    push!(subarea_plots, p)
end

# Create GIF
anim = @animate for p in subarea_plots[Not(1)]
    plot(p)
end
gif(anim, fps=3)

# Export GIF
# save_figures = true
if (@isdefined save_figures) && save_figures == true
    gif(anim, joinpath("fig", outcome, "05_$(outcome)_scaling.gif"); fps=3)
end

## 3 extents comparison

# Extract LCBD & relationship subplots for first, middle & last GIF plots
mid_ind = median(1:length(subarea_plots)) |> round |> Int64
ps = subarea_plots[[1, mid_ind, end]]

# Combine 3 scales in single figure
p = plot(
    deepcopy(ps)...;
    dpi=600,
    layout=(3, 1),
    size=(900, 960),
    title=["a) Regional extent" "" "b) Intermediate extent" "" "c) Continental extent" ""],
    titleloc=:left,
    bottommargin=-2.0mm,
)

# Tweak axes
yticks!(p[3], 34:3:50)
if outcome == "bart"
    yticks!(p[6], 1.5:0.5:5.0)
elseif outcome == "raw"
    yticks!(p[4], 0.5:0.5:2.5)
    yticks!(p[6], 2:1:9)
end
p

# Export figure
# save_figures = true
if (@isdefined save_figures) && save_figures == true
    savefig(p, joinpath("fig", outcome, "05_$(outcome)_extents.png"))
end

## Scaling medians figure

# Create empty elements
richness_medians = []
lcbd_medians = []
lcbd_mins = []
lcbd_maxs = []
lcbd_abs_medians = []
beta_values = []
gamma_values = []

# Get analysis values for all subareas
for sc in subarea_coords
    local ds = [d[sc] for d in distributions]
    local y = Ymatrix(ds)
    local r = richness(y, ds[1])
    local l = lcbd(y, ds[1]; relative=false)
    local l_abs = lcbd(y, ds[1]; relative=false)
    local bt = beta_total(y)
    local g = gamma(y)

    push!(richness_medians, median(r))
    push!(lcbd_medians, median(l))
    push!(lcbd_mins, minimum(l))
    push!(lcbd_maxs, maximum(l))
    push!(lcbd_abs_medians, l_abs)
    push!(beta_values, bt)
    push!(gamma_values, g)
end

# Check values
richness_medians
lcbd_medians
beta_values
gamma_values

# Transform to relative values (to plot all variables on same axis)
medians_df = DataFrame(
    idx=eachindex(richness_medians),
    richness=richness_medians ./ maximum(richness_medians),
    lcbd=lcbd_medians ./ maximum(lcbd_medians),
    beta=beta_values ./ maximum(beta_values),
    gamma=gamma_values ./ maximum(gamma_values)
)

# Plot values across scales (step-by-step plots)
medians_p1 = plot(
    medians_df.idx,
    medians_df.richness;
    label="Median richness",
    lw=2,
    xlabel="Subarea extent",
    ylabel="Subarea value (relative to maximum)",
    legend=:bottomright,
    xticks=:none,
    ylim=(0, 1),
    top_margin=mm
)
medians_p2 = plot!(deepcopy(medians_p1), medians_df.lcbd, label="Median LCBD", lw=2)
medians_p3 = plot!(deepcopy(medians_p2), medians_df.beta, label="Total beta diversity", lw=2)
medians_p4 = plot!(deepcopy(medians_p3), medians_df.gamma, label="Gamma diversity", lw=2)

# Export figure
# save_figures = true
if (@isdefined save_figures) && save_figures == true
    # Final plot
    savefig(medians_p4, joinpath("fig", outcome, "05_$(outcome)_medians.png"))
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
#! /usr/bin/julia

# Rosetta Code, Arithmetic-geometric mean

function agm{T<:FloatingPoint,U<:Integer}(x::T, y::T, e::U=5)
    0 < y && 0 < y && 0 < e || throw(DomainError())
    err = e*eps(x)
    (g, a) = extrema([x, y])
    while err < (a - g)
        ap = a
        a = 0.5*(a + g)
        g = sqrt(ap*g)
    end
    return a
end

x = 1.0
y = 1.0/sqrt(2.0)

println("Using literal-precision float numbers:")
println("    agm(", x, ",", y, ") = ", agm(x, y))

println()
println("Using half-precision float numbers:")
x = float16(x)
y = float16(y)
println("    agm(", x, ",", y, ") = ", agm(x, y))

println()
println("Using ", get_bigfloat_precision(), "-bit float numbers:")
x = BigFloat(1.0)
y = x/sqrt(BigFloat(2.0))
println("    agm(", x, ",", y, ") = \n    ", agm(x, y))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using CompScienceMeshes

"""
    Creates a predicate that checks that 2 edges belonging to `small` and `big`
    lie on the same plane and the normals of its surrounding cells points in the
    same direction
"""
function normal_cpredicate(small, big)
    @assert dimension(small) == 2
    @assert dimension(big) == 2

    ss = skeleton(small,1)
    sb = skeleton(big,1)

    D1 = connectivity(ss,small)
    D2 = connectivity(sb ,big)

    function pred(e1,e2)
        #TODO:Modify to take simplex rather than edge indices
        edgefaces1 = D1[:,e1]
        edgefaces2 = D2[:,e2]
        faceid1 = edgefaces1.nzind
        faceid2 = edgefaces2.nzind

        N1 = length(faceid1)
        N2 = length(faceid2)
        narray1 = Vector{CompScienceMeshes.StaticArrays.SArray{
                Tuple{3},Float64,1,3}}()
        narray2 = copy(narray1)

        for i = 1:N1
            push!(narray1,chart(small,small.faces[faceid1[i]]).normals[1])
        end
        @assert all(y->y==narray1[1],narray1) "egde on small mesh found to lie on junction"

        for i = 1:N2
            push!(narray2,chart(big,big.faces[faceid2[i]]).normals[1])
        end

        @assert all(y->y==narray2[1],narray2) "egde on big mesh found to lie on junction"

        return narray1[1] == narray2[1]
    end
end

"""
    Enforces extra constraint on edges in `small` and `big` that overlap with
    edges in `pInt` during `embedding` to ensure that both edges in `small` and
    `big` lie on planes with normals pointing in the same direction.
"""
function reembedding(small, big, pInt)
    @assert dimension(small) == 2
    @assert dimension(big) == 2
    @assert dimension(pInt) == 2

    ss = skeleton(small,1)
    sb = skeleton(big,1)
    sp = CompScienceMeshes.interior(pInt)

    epred = normal_cpredicate(small,big)
    overlaps = overlap_gpredicate(sp)
    mpred(e1,e2) = overlaps(chart(sb,cells(sb)[e2])) ? epred(e1,e2) : true
    #applies constraint during embedding only on edges that overlap with pInt
    CompScienceMeshes.embedding(ss,sb,predicate=mpred)
end

isclosed(m) = (numcells(boundary(m)) == 0)

plate = meshrectangle(1.0,1.0,1.0)
G01 = CompScienceMeshes.SComplex2D(plate)
PT = weld(G01, -G01, seam=boundary(G01));

A1 = reembedding( G01, PT, G01)
A2 = reembedding(-G01, PT, G01)

@show A1
@show A2


text:
let
  x = @match :(2+3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
  end
  @test x == (2, 3)
end

let
  x = @match :(2-3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
  end
  @test x == (3, 2)
end

let
  x = @match :(2/3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
  end
  @test x == nothing
end

let
  x = @match :(2/3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
    _ => :default
  end
  @test x == :default
end

let
  ex = :(mutable struct Foo
           x::Int
           y
         end)
  @capture(ex, mutable struct T_ fields__ end)
  @test T == :Foo
  @test fields == [:(x::Int), :y]

  @capture(ex, mutable struct T_ fields__ end)
  @test T == :Foo
  @test fields == [:(x::Int), :y]
end

let
  ex = :(f(x))
  @capture(ex, f_(xs__))
  @test f == :f
  @test xs == [:x]
end

let
  ex = :(f(x, y, z))
  @capture(ex, f_(x_, xs__))
  @test f == :f
  @test x == :x
  @test xs == [:y, :z]
end

let
  ex = quote
    function foo(a, b)
      return a+b
    end
  end
  @assert @capture(shortdef(ex), f_(args__) = body_)
end

let
  ex = :(a = b)
  @capture(ex, a_ = b_)
  @test (a, b) == (:a, :b)
end

let
  ex = :(f(a = b))
  @capture(ex, f(a_ = b_))
  @test (a, b) == (:a, :b)
  @capture(ex, f(x_))
  @test isexpr(x, :kw)
end

let
  ex = :(@foo(a,b))
  @capture(ex, @foo(a_,b_))
  @test (a, b) == (:a, :b)
end


text:
using Documenter
using SolverLogging

makedocs(
    sitename = "SolverLogging.jl",
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        ansicolor=true
    ),
    pages = [
        "Introduction" => "index.md",
        "API" => "api.md",
        "Examples" => "examples.md"
    ]
)

deploydocs(
    repo = "github.com/bjack205/SolverLogging.jl.git",
    devbranch = "master"
)


text:
using Logging
using Plots
# using XfromProjections.curve_utils
# using XfromProjections.snake_forward
# using XfromProjections
using Printf
using LinearAlgebra
include("./utils.jl")
include("../../phantoms/sperm_phantom.jl")
include("../snake_forward.jl")
include("../snake.jl")
include("../curve_utils.jl")

#Radius function
r(s) = 1.0
bins = collect(-38.0:0.125:38.0)
cwd = @__DIR__
frame_nr = 1
nangles = 1
angles = [π/2]
ang = angles[1]

images, tracks = get_sperm_phantom(11,grid_size=0.1)
tail_length = curve_lengths(tracks[end])[end]

#determine outline from skeleton
outline, normals = get_outline(tracks[frame_nr], r)
sinogram = parallel_forward(outline, angles, bins)

num_points = 30

template = get_straight_template(sinogram[:,1], r, [0.0 0.0], angles[1], num_points,bins)

plt = plot(template[:,1], template[:,2], aspect_ratio=:equal, label="template", framestyle=:none, color=:black, legend=false)
plot!(outline[:,1], outline[:,2], label = "target", color=:green)

w_u = ones(num_points+1)
w_u[1] = 0.0
w_u[2] = 0.0
w_l = ones(num_points+1)
w_l[1] = 0.0

stepsize = 0.1
best_residual = Inf
best_recon = deepcopy(template)

recon1 = recon2d_tail(deepcopy(best_recon),r,angles,bins,sinogram,5000, 0.0, stepsize, 1, w_u, zeros(num_points+1))
recon2 = recon2d_tail(deepcopy(best_recon),r,angles,bins,sinogram,5000, 0.0, stepsize, 1, zeros(num_points+1), w_l)

initial1 = deepcopy(recon1)
initial2 = deepcopy(recon2)
for flip_pt=1:num_points
    recon1_flipped = flip(initial1,flip_pt,angles[1])
    recon2_flipped = flip(initial2,flip_pt,angles[1])
    #mirror and reconstruct with weights on both sides
    recon1 = recon2d_tail(deepcopy(recon1_flipped),r,angles,bins,sinogram,100, 0.0, stepsize, 1, w_u, zeros(num_points+1))
    recon2 = recon2d_tail(deepcopy(recon2_flipped),r,angles,bins,sinogram,100, 0.0, stepsize, 1, zeros(num_points+1), w_l)

    global best_residual, best_recon[:,:], residual1, residual2 = try_improvement(best_residual, recon1, recon2, angles[1], bins, sinogram, best_recon, tail_length, r)
end

@info "keeping the best parts and restarting reconstruction"
recon_best = keep_best_parts(residual1, deepcopy(best_recon), ang, bins, 3, num_points, tail_length, sinogram[:,1], r)
recon1 = recon2d_tail(deepcopy(recon_best),r,[ang],bins,sinogram,5000, 0.0, stepsize, 1, w_u, zeros(num_points+1))
recon2 = recon2d_tail(deepcopy(recon_best),r,[ang],bins,sinogram,5000, 0.0, stepsize, 1, zeros(num_points+1), w_l)
best_residual,best_recon[:,:], residual1, residual2 = try_improvement(best_residual, recon1, recon2, ang, bins, sinogram, best_recon, tail_length, r)

initial1 = deepcopy(recon1)
initial2 = deepcopy(recon2)
for flip_pt=1:num_points
    recon1_flipped = flip(initial1,flip_pt,angles[1])
    recon2_flipped = flip(initial2,flip_pt,angles[1])
    #mirror and reconstruct with weights on both sides
    recon1 = recon2d_tail(deepcopy(recon1_flipped),r,angles,bins,sinogram,100, 0.0, stepsize, 1, w_u, zeros(num_points+1))
    recon2 = recon2d_tail(deepcopy(recon2_flipped),r,angles,bins,sinogram,100, 0.0, stepsize, 1, zeros(num_points+1), w_l)

    global best_residual, best_recon[:,:], residual1, residual2 = try_improvement(best_residual, recon1, recon2, angles[1], bins, sinogram, best_recon, tail_length, r)
end

recon1 = flip(best_recon,1,ang)
recon2 = deepcopy(best_recon)
recon2d_tail(deepcopy(recon1),r,angles,bins,sinogram,1, 0.0, stepsize, 1, w_u, zeros(num_points+1), doplot=true)
recon2d_tail(deepcopy(recon2),r,angles,bins,sinogram,1, 0.0, stepsize, 1, zeros(num_points+1), w_l, doplot=true)
plot!(recon1[:,1],recon1[:,2], label="upper", color=:red, linewidth=2)
plot!(recon2[:,1],recon2[:,2], label="lower", color=:blue, linewidth=2)

outline1, normals = get_outline(recon1, r)
s1 = parallel_forward(outline1, angles, bins)

outline2, normals = get_outline(recon2, r)
s2 = parallel_forward(outline2, angles, bins)

y_vals = cat(bins,s1,dims=2)
y_vals = rotate_points(translate_points(y_vals, [0.0 20.0]), angles[1])
plot_vals = y_vals[10:end-250,:]
plot!(plot_vals[:,1], plot_vals[:,2],label="FP upper", fill = (0, 0.2, :red), color=:white)

y_vals = cat(bins,s2,dims=2)
y_vals = rotate_points(translate_points(y_vals, [0.0 20.0]), angles[1])
plot_vals = y_vals[10:end-250,:]
plot!(plot_vals[:,1], plot_vals[:,2],label="FP lower", fill = (0, 0.2, :blue), color=:white)

y_vals = cat(bins,sinogram,dims=2)
y_vals = rotate_points(translate_points(y_vals, [0.0 20.0]), angles[1])
plot_vals = y_vals[10:end-250,:]
plot!(plot_vals[:,1], plot_vals[:,2],label="sinogram")

cwd = @__DIR__
savepath = normpath(joinpath(@__DIR__, "result"))
!isdir(savepath) && mkdir(savepath)
cd(savepath)
savefig("mirror_fixed")

cd(cwd)

plot!()

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
# #+HTML:

# #+HTML:

# #+HTML: MPS Contraction
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Matrix%20Product%20States][Matrix Product States:3]]
"""
    Base.:(*)(ψ′::Adjoint{T, MPS{L, T}}, ϕ::MPS{L, T}) where {L, T}
representing
    •--(b¹ b¹)--•--(b² b²)--•--(b³ b³)--•
    |           |           |           |
    σ′¹         σ′²         σ′³         σ′⁴
    σ′¹         σ′²         σ′³         σ′⁴
    |           |           |           |
    •--(a¹ a¹)--•--(a² a²)--•--(a³ a³)--•
"""
function Base.:(*)(ψ′::Adjoint{T, MPS{L, T}}, ϕ::MPS{L, T}) where {L, T}
    ψ = ψ′.parent

    M   = ϕ.tensors[1]
    M̃dg = dg(ψ.tensors[1])

    @tensor cont[b₁, a₁] := M̃dg[b₁, 1, σ₁] * M[1, a₁, σ₁]

    for i in 2:L-1
        M   = ϕ.tensors[i]
        M̃dg = dg(ψ.tensors[i])

        @tensor cont[bᵢ, aᵢ] := M̃dg[bᵢ, bᵢ₋₁, σᵢ] * cont[bᵢ₋₁, aᵢ₋₁] * M[aᵢ₋₁, aᵢ, σᵢ]
    end
    M   = ϕ.tensors[L]
    M̃dg = dg(ψ.tensors[L])

    @tensor M̃dg[1, bᴸ⁻¹, σᴸ] * cont[bᴸ⁻¹, aᴸ⁻¹] * M[aᴸ⁻¹, 1, σᴸ]
end
# Matrix Product States:3 ends here

# #+HTML: MPO Contraction
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Matrix%20Product%20Operators][Matrix Product Operators:2]]
"""
    Base.:(*)(O::MPO, ψ::MPS)
representing

    σ¹          σ²          σ³          σ⁴
    |           |           |           |
    •--(b¹ b¹)--•--(b² b²)--•--(b³ b³)--•
    |           |           |           |
    σ′¹         σ′²         σ′³         σ′⁴
    σ′¹         σ′²         σ′³         σ′⁴
    |           |           |           |
    •--(a¹ a¹)--•--(a² a²)--•--(a³ a³)--•
"""
function Base.:(*)(O::MPO{L, T}, ψ::MPS{L, T}) where {L, T}
    tensors = Array{T,3}[]
    for i in 1:L
        W = O.tensors[i]
        M = ψ.tensors[i]

        @reduce N[(bᵢ₋₁, aᵢ₋₁), (bᵢ, aᵢ), σᵢ] :=  sum(σ′ᵢ) W[bᵢ₋₁, bᵢ, σᵢ, σ′ᵢ] * M[aᵢ₋₁, aᵢ, σ′ᵢ]

        push!(tensors, N)
    end
    MPS{L, T}(tensors)
end

"""
    Base.:(*)(O1::MPO, O2::MPO)
representing

    σ¹          σ²          σ³          σ⁴
    |           |           |           |
    •--(b¹ b¹)--•--(b² b²)--•--(b³ b³)--•
    |           |           |           |
    σ′′¹        σ′′²        σ′′³        σ′′⁴
    σ′′¹        σ′′²        σ′′³        σ′′⁴
    |           |           |           |
    •--(a¹ a¹)--•--(a² a²)--•--(a³ a³)--•
    |           |           |           |
    σ′¹         σ′²         σ′³         σ′⁴
"""
function Base.:(*)(O1::MPO{L, T}, O2::MPO{L, T}) where {L, T}
    tensors = Array{T,4}[]
    for i in 1:L
        W1 = O1.tensors[i]
        W2 = O2.tensors[i]

        @reduce V[(bᵢ₋₁, aᵢ₋₁), (bᵢ, aᵢ), σᵢ, σ′ᵢ] :=  sum(σ′′ᵢ) W1[bᵢ₋₁, bᵢ, σᵢ, σ′′ᵢ] * W2[aᵢ₋₁, aᵢ, σ′′ᵢ, σ′ᵢ]

        push!(tensors, V)
    end
    MPO{L, T}(tensors)
end

"""
    Base.:(*)(ψ::Adjoint{T,MPS{L,T}}, O::MPO) where {L,T}
representing

    •--(a¹ a¹)--•--(a² a²)--•--(a³ a³)--•
    |           |           |           |
    σ′¹         σ′²         σ′³         σ′⁴
    σ′¹         σ′²         σ′³         σ′⁴
    |           |           |           |
    •--(b¹ b¹)--•--(b² b²)--•--(b³ b³)--•
    |           |           |           |
    σ¹          σ²          σ³          σ⁴
"""
function Base.:(*)(ψ′::Adjoint{T,MPS{L,T}}, O::MPO{L, T}) where {L,T}
    ψ = ψ′.parent
    tensors = Array{T,3}[]
    Ws = dg.(reverse(O.tensors))
    for i in 1:L
        W = Ws[i]
        M = ψ.tensors[i]

        @reduce N[(bᵢ₋₁, aᵢ₋₁), (bᵢ, aᵢ), σᵢ] :=  sum(σ′ᵢ) W[bᵢ₋₁, bᵢ, σᵢ, σ′ᵢ] * M[aᵢ₋₁, aᵢ, σ′ᵢ]
        push!(tensors, N)
    end
    adjoint(MPS{L, T}(tensors))
end
# Matrix Product Operators:2 ends here


text:
using ISA #References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

fs = 10000
t = collect(0.0:1/fs:0.1)

a₀(t) = exp(-(t)^2)
ω₀(t) = 1000
φ₀ = 0
𝐶₀ = Tuple([a₀,ω₀,φ₀])

a₁(t) = exp(-1000.0(t-0.05)^2)
ω₁(t) = 2*π*1000 + 1500*sin(50t)
φ₁ = π
𝐶₁ = Tuple([a₁,ω₁,φ₁])

a₂(t) = 0.8*cos(2t)
ω₂(t) = 2*π*(-500 + 5000t)
φ₂ = π
𝐶₂ = Tuple([a₂,ω₂,φ₂])

𝑆 = [𝐶₀,𝐶₁,𝐶₂]
z = AMFMmodel(𝑆)

p1 = isaPlot3d(z, t, backend="Plotly")

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
module InsuranceContracts
import BitemporalPostgres
import SearchLight: DbId
import Base: @kwdef
export Contract, ContractRevision, ContractPartnerRef, ContractPartnerRefRevision, ProductItem, ProductItemRevision, ProductItemTariffRef, ProductItemTariffRefRevision, ProductItemPartnerRef, ProductItemPartnerRefRevision
using BitemporalPostgres

"""
Contract

  a contract component of a bitemporal entity

"""
@kwdef mutable struct Contract <: BitemporalPostgres.Component
  id::DbId = DbId()
  ref_history :: DbId = InfinityKey
  ref_version :: DbId = InfinityKey
end

"""
ContractRevision

  a revision of a contract component

"""
@kwdef mutable struct ContractRevision <: BitemporalPostgres.ComponentRevision
  id::DbId = DbId()
  ref_component :: DbId = InfinityKey
  ref_validfrom::DbId = InfinityKey
  ref_invalidfrom::DbId = InfinityKey
  description::String = ""
end

"""
ProductItem

  a productitem component of a contract component

"""
@kwdef mutable struct ProductItem <: BitemporalPostgres.SubComponent
  id::DbId = DbId()
  ref_history :: DbId = InfinityKey
  ref_version :: DbId = InfinityKey
  ref_super :: DbId = InfinityKey
end

"""
ProductItemRevision

  a revision of a productitem component

"""
@kwdef mutable struct ProductItemRevision <: BitemporalPostgres.ComponentRevision
  id::DbId = DbId()
  ref_component :: DbId = InfinityKey
  position::Integer = 0 :: Int64
  ref_validfrom::DbId = InfinityKey
  ref_invalidfrom::DbId = InfinityKey
  description::String = ""
end

"""
ContractPartnerRef

  a partner reference of a contract component, i.e. policy holder, premium payer

"""
@kwdef mutable struct ContractPartnerRef <: BitemporalPostgres.SubComponent
  id::DbId = DbId()
  ref_history :: DbId = InfinityKey
  ref_version :: DbId = InfinityKey
  ref_super :: DbId = InfinityKey
end

"""
ContractPartnerRefRevision

  a revision of a contract's partner reference

"""
@kwdef mutable struct ContractPartnerRefRevision <: BitemporalPostgres.ComponentRevision
  id::DbId = DbId()
  ref_component :: DbId = InfinityKey
  ref_validfrom::DbId = InfinityKey
  ref_invalidfrom::DbId = InfinityKey
  description::String = ""
  ref_partner::DbId = InfinityKey
end

"""
ProductItemTariffRef

  a reference to a tariff of a productitem

"""
@kwdef mutable struct ProductItemTariffRef <: BitemporalPostgres.SubComponent
  id::DbId = DbId()
  ref_history :: DbId = InfinityKey
  ref_version :: DbId = InfinityKey
  ref_super :: DbId = InfinityKey
end

"""
ProductItemTariffRefRevision

  a revision of a productitem's reference to a tariff

"""
@kwdef mutable struct ProductItemTariffRefRevision  <: BitemporalPostgres.ComponentRevision
  id::DbId = DbId()
  ref_component :: DbId = InfinityKey
  ref_validfrom::DbId = InfinityKey
  ref_invalidfrom::DbId = InfinityKey
  description::String = ""
  ref_tariff::DbId = InfinityKey
end

"""
ProductItemPartnerRef

  a reference to a partner of a productitem, i.e. insured person

"""
@kwdef mutable struct ProductItemPartnerRef <: BitemporalPostgres.SubComponent
  id::DbId = DbId()
  ref_history :: DbId = InfinityKey
  ref_version :: DbId = InfinityKey
  ref_super :: DbId = InfinityKey
end

"""
ProductItemPartnerRefRevision

  a revision of a productItem's partner reference

"""
@kwdef mutable struct ProductItemPartnerRefRevision  <: BitemporalPostgres.ComponentRevision
  id::DbId = DbId()
  ref_component :: DbId = InfinityKey
  ref_validfrom::DbId = InfinityKey
  ref_invalidfrom::DbId = InfinityKey
  description::String = ""
  ref_partner::DbId = InfinityKey
end

end

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
# Compared to cudnnAddTensor!(copy(a),b), cudnnOpTensor is ~50% faster on
# (14,14,256,32)+(1,1,256,1), ~50% slower on (1,1,100,100)+(1,1,100,1) Unlike cudnnAddTensor
# it supports all broadcasting shapes up to ndims=5 as described in the documentation.

"""
    cudnnOpTensor(x1, x2; op, compType, nanOpt, alpha1, alpha2)
    cudnnOpTensor(x1, x2, d::cudnnOpTensorDescriptor; alpha1, alpha2)
    cudnnOpTensor!(y, x1, x2; op, compType, nanOpt, alpha1, alpha2, beta)
    cudnnOpTensor!(y, x1, x2, d::cudnnOpTensorDescriptor; alpha1, alpha2, beta)

Return the result of the specified broadcasting operation applied to `x1` and `x2`.
Optionally `y` holds the result and `d` specifies the operation. Each dimension of the input
tensor `x1` must match the corresponding dimension of the destination tensor `y`, and each
dimension of the input tensor `x2` must match the corresponding dimension of the destination
tensor `y` or must be equal to 1. Keyword arguments:

* `alpha1=1, alpha2=1, beta=0` are used for scaling, i.e. `y .= beta*y .+ op.(alpha1*x1, alpha2*x2)`

Keyword arguments used when `cudnnOpTensorDescriptor` is not specified:

* `op = CUDNN_OP_TENSOR_ADD`, ADD can be replaced with MUL, MIN, MAX, SQRT, NOT; SQRT and NOT performed only on x1; NOT computes 1-x1
* `compType = (eltype(x1) <: Float64 ? Float64 : Float32)`: Computation datatype (see cudnn docs for available options)
* `nanOpt = CUDNN_NOT_PROPAGATE_NAN`: NAN propagation policy. The other option is `CUDNN_PROPAGATE_NAN`.
"""
cudnnOpTensor, cudnnOpTensor!

# Public methods:
cudnnOpTensor(x1,x2; o...)    = cudnnOpTensorWithDefaults(x1,x2; o...)
cudnnOpTensor!(y,x1,x2; o...) = cudnnOpTensorWithDefaults(x1,x2; y, o...)
cudnnOpTensor(x1,x2,d::cudnnOpTensorDescriptor; o...)    = cudnnOpTensorWithDefaults(x1,x2; opTensorDesc=d, o...)
cudnnOpTensor!(y,x1,x2,d::cudnnOpTensorDescriptor; o...) = cudnnOpTensorWithDefaults(x1,x2; y, opTensorDesc=d, o...)

# Private method:
function cudnnOpTensorWithDefaults(
    x1, x2;
    y = similar(x1),
    op::cudnnOpTensorOp_t = CUDNN_OP_TENSOR_ADD,
    compType::DataType = (eltype(x1) <: Float64 ? Float64 : Float32),
    nanOpt::cudnnNanPropagation_t = CUDNN_NOT_PROPAGATE_NAN,
    opTensorDesc::cudnnOpTensorDescriptor = cudnnOpTensorDescriptor(op, cudnnDataType(compType), nanOpt),
    alpha1::Real = 1,
    alpha2::Real = 1,
    beta::Real = 0,
    x1Desc::cudnnTensorDescriptor = cudnnTensorDescriptor(x1),
    x2Desc::cudnnTensorDescriptor = cudnnTensorDescriptor(x2),
    yDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(y)
)
    @assert ndims(x1) <= 5
    @assert size(y) == size(x1)
    @assert all(size(x2,i) == size(x1,i) || size(x2,i) == 1 for i in 1:ndims(x2))
    T = eltype(x1)
    alpha1, alpha2, beta = scalingParameter(T,alpha1), scalingParameter(T,alpha2), scalingParameter(T,beta)
    cudnnOpTensorAD(x1, x2; opTensorDesc, alpha1, x1Desc, alpha2, x2Desc, beta, yDesc, y)
end

# AD method: This method aids gradient definition, please do not remove!
function cudnnOpTensorAD(x1, x2; opTensorDesc, alpha1, x1Desc, alpha2, x2Desc, beta, yDesc, y)
    cudnnOpTensor(handle(), opTensorDesc, alpha1, x1Desc, x1, alpha2, x2Desc, x2, beta, yDesc, y)
    return y
end

# Deprecated:
function cudnnOpTensor(op::cudnnOpTensorOp_t,
                       A::DenseCuArray{T,N}, B::DenseCuArray{T,N}, C::DenseCuArray{T,N};
                       alpha1=true, alpha2=true, beta=false) where {T,N}
    @warn "cudnnOpTensor(op,A,B,C) is deprecated, please use one of the methods in `@doc cudnnOpTensor`." maxlog=1
    cudnnOpTensorWithDefaults(A, B; y=C, op, alpha1, alpha2, beta)
end


text:
#Get p-values for each number of isoforms for gene. Will this work for global?
function FindpValues(InputDict, mean_val, threshold)

    OutputDict = Dict()

    for (key, value) in InputDict

        #Remove p-values for insufficiently continuous distributions
        if length(unique(value)) < threshold
            OutputDict[key] = missing
        else
            pValue = sum(x -> x > mean_val, value) / length(value)
            OutputDict[key] = pValue
        end
    end

    return OutputDict
end

#Get p-value given an array and mean value
function GlobalpValues(meanArr, mean_val, threshold)

    if length(unique(meanArr)) < threshold
        pValue = missing
    else
        pValue = sum(x -> x > mean_val, meanArr) / length(meanArr)
    end

    return pValue
end

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
struct ArgLooper  <: AbstractArgIter
    dict_list::Vector{Vector{String}}
    runs_iter::AbstractArray
    stable_arg::Vector{String}
    done::Bool
    run_name::String
    ArgLooper(list, stable_arg, runs_iter, run_name="--run") = new(list, runs_iter, stable_arg, false, run_name)
end

function make_arguments(iter::ArgLooper, state)
    arg_list = Vector{String}()
    arg_list = [iter.dict_list[state[2][1]]; [iter.run_name, string(iter.runs_iter[state[2][2]])]; iter.stable_arg]
    return arg_list
end

function Base.iterate(iter::ArgLooper)
    state = (1, [1, 1])
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function next_state(iter::ArgLooper, _state)
    state = _state[2]
    n_state = _state[1]

    state[2] += 1
    if state[2] > length(iter.runs_iter)
        state[2] = 1
        state[1] += 1
    end
    return (n_state+1, state)

end

function Base.iterate(iter::ArgLooper, state)
    if state[2][1] > length(iter.dict_list)
        return nothing
    end
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function Base.length(iter::ArgLooper)
    return length(iter.dict_list)*length(iter.runs_iter)
end


text:
# This file was generated, do not modify it. # hide
PCA = @load PCA
pca_pipe = Standardizer() |> PCA(maxoutdim=2)
pca = machine(pca_pipe, Xtrain)
fit!(pca)
W = transform(pca, Xtrain)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
struct TensorIntegrationGrid{T1,T2}
    quads::T1
    quad_indices::T2

    function TensorIntegrationGrid(dim_, number_polys::Tuple, vars; wType = Float64)
        @assert length(number_polys) == length(vars)

        quads_temp = []

        for i in Base.OneTo(length(vars))
            @assert number_polys[i] > 0
            quad_list = []
            var = vars[i]
            if var isa Uniform
                @assert var.a ≈ -1.0 && var.b ≈ 1.0
                push!(quads_temp, GaussLegendreQuad(number_polys[i]; wType))
            elseif var isa Normal
                # Using physicist's Hermite polynomial
                @assert var.μ ≈ 0.0 && var.σ ≈ sqrt(2.0) / 2
                push!(quads_temp, GaussHermiteQuad(number_polys[i]; wType))
            else
                throw("Not implemented")
            end
        end

        quads = Tuple(quads_temp)
        quad_indices = CartesianIndices(Tuple(number_polys))

        return new{typeof(quads),typeof(quad_indices)}(quads, quad_indices)
    end

end

max_stars_repo_name: d-netto/PolyChaosODE

max_issues_repo_name: d-netto/PolyChaosODE

max_forks_repo_name: d-netto/PolyChaosODE


text:
## Simple PH model with Weibull baseline hazard

struct WeibullPHData <: AbstractData
    id::Vector{Int}
    entry::Vector{Float64}
    time::Vector{Float64}
    status::Vector{Bool}
    z::Matrix{Float64}
end

WeibullPHData() = WeibullPHData(Vector{Int}(undef, 0), Vector{Float64}(undef, 0), Vector{Float64}(undef, 0),
              Vector{Bool}(undef, 0), Matrix{Float64}(undef, 0, 0))

mutable struct WeibullPHParam{T <: Real} <: AbstractParam
    shape::T
    scale::T
    beta::Vector{T}
end

struct WeibullPHFixedParam <: AbstractFixedParam
    n::Int
    tmax::Float64
    max_entry::Float64
end

struct WeibullPHPrior <: AbstractPrior
    shape::ContinuousUnivariateDistribution
    scale::ContinuousUnivariateDistribution
    beta::MultivariateDistribution
end

WeibullPHPrior(k=1) = WeibullPHPrior(Gamma(1, 1), Gamma(1, 1), MvNormal(diagm(0 => 10*ones(k))))

## exponential model
WeibullPHModel() = Model(WeibullPHParam(1.0, 1.0, [0.0]), WeibullPHFixedParam(100, 55.0, 48.0), ph_exp_1)

WeibullPHModel(param, fparam, gen_data=ph_exp_1) = Model(param, fparam, gen_data)

loglik(data::WeibullPHData, param::WeibullPHParam) = log_likelihood(param.shape, param.scale, param.beta, data)

struct WeibullPHFit <: AbstractFit
    shape::Vector{Float64} ## Weibull shape
    scale::Vector{Float64} ## Weibull scale
    beta::Matrix{Float64} ## hazard ratio
    accept::Float64 ## acceptance rates
    prior::Prior
end

WeibullPHFit() = WeibullPHFit(Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Float64}(undef,0), Vector{Float64}(undef,0), Prior())

function slice(i::Int, x::WeibullPHFit)
    if length(x.shape) == 0
        sample_prior(1, x.prior)[1]
    else
        WeibullPHParam(x.shape[i], x.scale[i], x.beta[i,:])
    end
end

function sample_prior(B::Int, prior::WeibullPHPrior)
    shape = rand(prior.shape, B)
    scale = rand(prior.scale, B)
    beta = rand(prior.beta, B)
    [WeibullPHParam(shape[b], scale[b], beta[:,b]) for b in 1:B]
end

## return prior distribution (before sampling) or posterior sample (after sampling)
## of treatment effect parameter
function get_outcome(x::WeibullPHFit)
    if size(x.beta, 1) == 0
        x.prior.beta
    else
        x.beta
    end
end

function predict_data(par::WeibullPHParam, fpar::WeibullPHFixedParam, data0::WeibullPHData, caltime::Float64, nevents=NaN)
    n = fpar.n
    tmax = fpar.tmax
    max_entry = n / 15 ##fpar.max_entry

    beta = par.beta
    shape = par.shape
    scale = par.scale
    n0 = length(data0.time)

    J = length(par.beta)
    group = Matrix{Float64}(undef,n,J)
    z = sample(0:J, n)

    for j in 1:J
        group[:,j] = (z .== j) .- 0.5
    end

    if n < n0
        error("requested sample size must be at least as large as sample size of old data")
    end

    entry = Vector{Float64}(undef, n)
    chaz = zeros(Float64, n)

    if caltime >= max_entry
        if n > n0
            error("max_entry has passed, cannot recruit more patients")
        end
    else
        copyto!(entry, rand(Uniform(caltime, max_entry), n))
    end

    left = caltime .- entry

    ## keep old data
    if n0 > 0
        for j in 1:J
            group[1:n0,j] = data0.z[:,j]
        end
        entry[1:n0] = data0.entry
        chaz[1:n0] = scale*exp.(data0.z*beta).*data0.time.^shape
        left[1:n0] = caltime .- data0.entry
        keep = data0.time .< left[1:n0]
        chaz[1:n0][keep] .= 0.0
    end

    U = .-log.(1 .-rand(Uniform(0, 1), n)) .+ chaz
    T = (U ./ (exp.(group*beta).*scale)).^(1/shape)

    time = min.(T, tmax)
    status = T .< tmax

    max_id = 0
    id = zeros(Int, n)

    ## keep old data
    if n0 > 0
        keep = data0.time .< left[1:n0]

        if sum(keep) > 0
            keep2 = Vector{Bool}(undef, n)
            fill!(keep2, false)
            keep2[1:n0] = keep

            time[keep2] = data0.time[keep]
            status[keep2] = data0.status[keep]
        end

        id[1:n0] = data0.id
        max_id = maximum(data0.id)
    end

    if n > n0
        id[(n0+1):end] = (max_id+1):(max_id+n-n0)
    end

    WeibullPHData(id, entry, time, status, group)
end

function ph_exp_1(par::WeibullPHParam, fpar::WeibullPHFixedParam)
    n = fpar.n

    z = zeros(Float64, n, 1)
    z[:,1] = rand(Binomial(1, 0.5), n) .- 0.5

    entry = rand(Uniform(0, fpar.max_entry), n)
    U = .-log.(1 .-rand(Uniform(0, 1), n))
    T = (U./(par.scale .* exp.(z*par.beta))).^(1/par.shape)
    ##T = U.^(1/par.shape).*exp.(-par.beta*z).*par.scale
    ##C = min.(rand(Exponential(fpar.tmax), n), fpar.tmax)
    C = min.(rand(Exponential(fpar.tmax), n), fpar.tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:fpar.n

    return Data(id, entry, time, status, z)
end

function ph_weibull_1(par, fpar)
    n = fpar.n

    z = rand(Binomial(1, 0.5), n) .- 0.5

    ## 25 patients per month
    ##max_entry = n / 25
    entry = rand(Uniform(0, fpar.max_entry), n)
    ##U = -log.(1 .-rand(Uniform(0, 1), n))
    U = 1 .- rand(Uniform(0, 1), n)

    shape = 1.25
    scale = log(2)/22.5

    T = [quantile(Gamma(6, 2), 1 .- exp.(log(U[i]).*exp.(.-par.beta .* z[i]))) for i in 1:n]

    ##C = min.(rand(Exponential(fpar.tmax), n), fpar.tmax)
    C = min.(rand(Uniform(0, 3*36), n), fpar.tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:fpar.n

    Data(id, entry, time, status, z)
end

function select_data(t::Float64, data::WeibullPHData; entry_max=t)
    sel = data.entry .<= entry_max
    entry2 = data.entry[sel]
    time2 = data.time[sel]
    status2 = data.status[sel]

    cutoff = max.(t .- entry2, 0)

    status2 = status2 .& (time2 .<= cutoff)
    time2 = min.(time2, cutoff)

    Data(data.id[sel], entry2, time2, status2, data.z[sel,:])
end

S(time, z, lambda, beta) = exp.(-(time ./ (scale .* exp.(-z*beta)))^shape)

function log_likelihood(shape, scale, beta, data)
    time = data.time
    status = data.status
    psi = exp.(data.z*beta)

    if (shape .<= 0.0) | (scale .<= 0.0)
        -Inf
    else
        sum(status.*(log(shape) + log(scale) .+ log.(psi) .+ (shape-1).*log.(time)) .- (scale .* psi .* time.^shape))
    end
end

function laplace_approx(data::WeibullPHData, param::WeibullPHParam, prior::WeibullPHPrior)
    f(x) = -log_likelihood(exp(x[1]), exp(x[2]), x[3:end], data) - logpdf(prior.shape, exp(x[1])) - x[1] - logpdf(prior.scale, exp(x[2])) - x[2] - logpdf(prior.beta, x[3:end])

    res = optimize(f, [log(param.shape), log(param.scale), param.beta...])
    x = Optim.minimizer(res)
    est_log_shape = x[1]
    est_log_scale = x[2]
    est_beta = x[3:end]

    ## hessian of -loglik
    H = ForwardDiff.hessian(f, x)

    H_chol = cholesky!(H; check=false)
    if !issuccess(H_chol)
        [DiracPM(est_beta[1])]
    else
        L = inv(H_chol.U')
        ##MvNormal(est_beta, (L'*L)[3:end,3:end])
        [Normal(est_beta[1], sqrt((L'*L)[end,end]))]
    end
end

function sample_param(shape, scale, beta, data, V, t, prior)
    target(x) = log_likelihood(exp(x[1]), exp(x[2]), x[3:end], data) + logpdf(prior.shape, exp(x[1])) + x[1] + logpdf(prior.scale, exp(x[2])) + x[2] + logpdf(prior.beta, x[3:end])
    AM_RR_single(target, [shape, scale, beta...], V, t)
end

function sample_param(shape, scale, data, V, t, prior)
    target(x) = log_likelihood(exp(x[1]), exp(x[2]), prior.beta.x, data) + logpdf(prior.shape, exp(x[1])) + x[1] + logpdf(prior.scale, exp(x[2])) + x[2]
    AM_RR_single(target, [shape, scale], V, t)
end

function get_approx(param::WeibullPHParam, new_data::WeibullPHData, prior::Prior)
    laplace_approx(new_data, param, prior)
    ##x = fit(10000, new_data; prior=prior)
    ##x.beta
    ##Normal(mean(x.beta), sqrt(var(x.beta)))
end

function fit(M::Int, data::WeibullPHData; prior::WeibullPHPrior=WeibullPHPrior(), warmup::Int=Int(floor(M/2)), use_NUTS=false)
    if M == 0 | length(data.time) == 0
        return WeibullPHFit(Vector{Float64}(undef, 0), Vector{Float64}(undef, 0), Matrix{Float64}(undef, 0, 0), 0.0, prior)
    end

    J = length(prior.beta)
    beta = Matrix{Float64}(undef,M+1,J)
    beta[1,:] = zeros(Float64, J)
    shape = ones(Float64, M+1)
    scale = ones(Float64, M+1)

    if prior.beta isa MultivariateDiracPM
        V = diagm(0 => ones(Float64, 2))
        mb = [log(shape[1]), log(scale[1])]
    else
        V = diagm(0 => ones(Float64, 2+J))
        mb = [log(shape[1]), log(scale[1]), beta[1,:]...]
    end

    eps = 0.001*V
    attempt = 0
    accept = 0
    index = (warmup+1):M

    for t in 1:M
        attempt += 1
        if prior.beta isa MultivariateDiracPM
            x, ac = sample_param(log(shape[t]), log(scale[t]), data, V, t, prior)
            beta[t+1,:] = prior.beta.x
        else
            x, ac = sample_param(log(shape[t]), log(scale[t]), beta[t,:], data, V, t, prior)
            beta[t+1,:] = x[3:end]
        end

        shape[t+1] = exp(x[1])
        scale[t+1] = exp(x[2])

        new_mb = (x + mb*t) / (t+1)
        V = updateVar(t+1, x, new_mb, mb, V, eps)
        mb = new_mb
        accept += ac
    end

    WeibullPHFit(shape[index], scale[index], beta[index,:], accept/attempt, prior)
end

function fit(chains::Int, M::Int, data::WeibullPHData; prior::WeibullPHPrior=WeibullPHPrior(), warmup::Int=Int(floor(M/2)), use_NUTS=false)
    @distributed vcat for k in 1:chains
        fit(M, data; warmup=warmup, prior=prior, use_NUTS=use_NUTS)
    end
end

function merge_chains(a::WeibullPHFit, b::WeibullPHFit)
    WeibullPHFit(vcat(a.shape, b.shape), vcat(a.scale, b.scale), vcat(a.beta, b.beta), 0.5*(a.accept .+ b.accept), a.prior)
end

merge_chains(a) = a
merge_chains(a...) = foldr(merge_chains, a)

function summary(x::Vector{WeibullPHFit})
    chains = length(x)

    z = merge_chains(x...)
    df = summary(z)

    if chains > 1
        len = Int(floor(length(x[1].shape)/2))
        df.Rhat = [potential_scale_reduction([z.beta[(1+(k-1)*len):(k*len), j] for k in 1:(2*chains) for j in 1:size(z.beta,2)]...),
                   potential_scale_reduction([z.shape[(1+(k-1)*len):(k*len)] for k in 1:(2*chains)]...),
                   potential_scale_reduction([z.scale[(1+(k-1)*len):(k*len)] for k in 1:(2*chains)]...)]
    end

    df
end

function summary(x::WeibullPHFit)
    J = size(x.beta,2)

    DataFrame(variable=[repeat(["beta"], inner=J)..., "shape", "scale"],
              mean=[mean(x.beta,dims=1)[1,:]..., mean(x.shape), mean(x.scale)],
              variance=[var(x.beta,dims=1)[1,:]..., var(x.shape), var(x.scale)],
              neff=[[effective_sample_size(x.beta[:,j]) for j in 1:J]..., effective_sample_size(x.shape), effective_sample_size(x.scale)])
end

function sim(B::Int, M::Int; warmup=Int(floor(M/2)), use_NUTS=false,
             par::WeibullPHParam=WeibullPHParam(1.2, 1.5, [log(0.77)]),
             fpar::WeibullPHFixedParam=WeibullPHFixedParam(100, 55, 48))
    res = Array{Float64, 2}(undef, B, 4)

    J = length(par.beta)
    beta = Matrix{Float64}(undef, B,J)
    var_beta = Matrix{Float64}(undef,B,J)
    cpb = Matrix{Bool}(undef, B,J)

    shape = Vector{Float64}(undef, B)
    scale = Vector{Float64}(undef, B)
    var_shape = Vector{Float64}(undef, B)
    var_scale = Vector{Float64}(undef, B)
    cpsh = Vector{Bool}(undef, B)
    cpsc = Vector{Bool}(undef, B)

    covered(x, v, b) = ((x - 1.96*sqrt(v)) < b) & (b < (x + 1.96*sqrt(v)))

    prior = Prior(Gamma(1, 1), Gamma(1, 1), MultivariateDiracPM([0.0]))

    for b in 1:B
        data = ph_exp_1(par, fpar)
        data = select_data(fpar.tmax, data)
        x = fit(M, data, prior=prior, warmup=warmup, use_NUTS=use_NUTS)

        ##@show x.accept
        ##@info effective_sample_size(x.beta), effective_sample_size(x.shape), effective_sample_size(x.scale)

        beta[b,:] = mean(x.beta, dims=1)
        var_beta[b,:] = var(x.beta, dims=1)
        cpb[b,:] = [(quantile(x.beta[:,j], 0.025) < par.beta[j]) .& (quantile(x.beta[:,j], 0.975) > par.beta[j]) for j in 1:length(par.beta)]

        shape[b] = mean(x.shape)
        var_shape[b] = var(x.shape)
        cpsh[b] = covered(shape[b], var_shape[b], par.shape)

        scale[b] = mean(x.scale)
        var_scale[b] = var(x.scale)
        cpsc[b] = covered(scale[b], var_scale[b], par.scale)
    end

    DataFrame(variable=["shape", "scale", repeat(["beta"], inner=J)...],
              mean=[mean(shape), mean(scale), mean(beta,dims=1)[1,:]...],
              sse=sqrt.([var(shape), var(scale), var(beta,dims=1)[1,:]...]),
              ese=sqrt.([mean(var_shape), mean(var_scale), mean(var_beta,dims=1)[1,:]...]),
              cp=[mean(cpsh), mean(cpsc), mean(cpb,dims=1)[1,:]...])
end

function weibull_test_laplace_approx()
    par = WeibullPH.Param(1.5, 1.0, [log(0.78)])
    fpar = WeibullPH.FixedParam(400, 55.0, 48.0)
    data = WeibullPH.ph_exp_1(par, fpar)

    x = WeibullPH.fit(10000, data)
    la = WeibullPH.laplace_approx(data, par, WeibullPH.Prior())

    @show WeibullPH.summary(x)
    @show (mean(la), var(la))
end

function weibull_test_chains(M1, M2, chains)
    model = WeibullPH.WeibullPHModel()

    data = model.generate_data(model.param, model.fparam)
    x = WeibullPH.fit(chains, M1, data, use_NUTS=false)
    summary(x)
end


text:
"""
        matrix2graph(sparse_matrix, [partition_by_rows::Bool=true])

A utility function to generate a graph from input
sparse matrix, columns are represented with vertices
and 2 vertices are connected with an edge only if
the two columns are mutually orthogonal.
"""
function matrix2graph(sparse_matrix::SparseMatrixCSC{<:Number, Int}, partition_by_rows::Bool=true)

    dropzeros!(sparse_matrix)
    (rows_index, cols_index, _) = findnz(sparse_matrix)

    ncols = size(sparse_matrix, 2)
    nrows = size(sparse_matrix, 1)

    num_vtx = partition_by_rows ? nrows : ncols

    inner = SimpleGraph(num_vtx)

    if partition_by_rows
        @inbounds for i in eachindex(rows_index)
            cur_row = rows_index[i]
            for j in 1:(i-1)
                next_row = rows_index[j]
                if cols_index[i] == cols_index[j] && cur_row != next_row
                    add_edge!(inner, cur_row, next_row)
                end
            end
        end
    else
        @inbounds for i in eachindex(cols_index)
            cur_col = cols_index[i]
            for j in 1:(i-1)
                next_col = cols_index[j]
                if rows_index[i] == rows_index[j] && cur_col != next_col
                    add_edge!(inner, cur_col, next_col)
                end
            end
        end
    end
    return VSafeGraph(inner)
end


text:
@testset "Greedy Coloring" begin

    g3 = StarGraph(10)

    for g in testgraphs(g3)
        for op_sort in (true, false)
            C = @inferred(greedy_color(g, reps=5, sort_degree=op_sort))
            @test C.num_colors == 2
        end
    end

    g4 = PathGraph(20)
    g5 = CompleteGraph(20)

    for graph in [g4, g5]
        for g in testgraphs(graph)
            for op_sort in (true, false)
                C = @inferred(greedy_color(g, reps=5, sort_degree=op_sort))

                @test C.num_colors <= maximum(degree(g))+1
                correct = true
                for e in edges(g)
                    C.colors[src(e)] == C.colors[dst(e)] && (correct = false)
                end
                @test correct
            end
        end
    end
end


text:
module includes

import Random
import REPL
import TOML
using Dates

"""
note: works on mac & windows
safe_path credits to: the Julia Programming Team (from [Pkg.jl/src/manifest.jl](https://github.com/JuliaLang/Pkg.jl/blob/master/src/manifest.jl)):
(ordered as-is by github.com) [10 members]:

Kristoffer Carlsson     @KristofferC
Ian Butterworth         @IanButterworth #PRO
David Varela            @00vareladavid
Tim Holy                @timholy #PRO [neuroscientist, professor, julia tutor]
Gunnar Farnebäck        @GunnarFarneback
Fredrik Ekre            @fredrikekre #PRO [postdoc]
Dilum Aluthge           @DilumAluthge
Tamas K. Papp           @tpapp
David Widmann           @devmotion #PRO
Shuhei Kadowaki         @aviatesk

"""
# turn /-paths into \-paths on Windows
"""safe_path
source:
https://github.com/JuliaLang/Pkg.jl/blob/b4da4946735fe4c7c6bb3e2bc16af95b4e76e487/src/manifest.jl#L57

helps in resolving uris' slashes for windows users

"""
function safe_path(path::String)
    if Sys.iswindows() && !isabspath(path)
        path = joinpath(split(path, "/")...) # what if: if we can do all urls as "\": thus we'll only Modify for windows users while other OSs can communicate too worry-free  -newly standardized
    end
    return path
end
paths = ["functions", "functions/Generators", "Distributions", "NeuralNet/Networks", "src/DataStructures/Arrays/Array.jl", "tests"] #Idea: all string paths in here # used in a for loop - in  a custom iterable function #parallizable
#TODO: Iterate all locations in paths
#include(safe_path(".")) #why refeing to thyself?
#ERROR:PermissionDenied #TODO: find a solution on windows pc to solve previlages in julia
include(safe_path("functions"))
include(safe_path("functions/Generators"))
include(safe_path("Distributions"))
include(safe_path("NeuralNet/Networks"))
include(safe_path("tests"))
include(safe_path("src/constants.jl"))
include(safe_path("src/DataStructures/Arrays/Array.jl"))

"""
to be called upon path error
"""
#= #systemError: Permission Denied (on Windows)#TODO: check common # solution?  #nice try !, needs further checking to see why it defaults to last catch

function pathRevise(path::String) # TRY TO FIX (possible) url typos
    strings = split(path, "/")
    try
        strings[1] = Upper(strings[1])#Gatcha#1: try Uppercasing the first letter  #(input is a miniscule, while actual directory is Majiscule)
        path = joinpath(strings...)
    catch
    end
    #=for i in enumerate length(strings)
        #TODO: do minor changes to path - i.e. uppercase other words #TODO: how to spot others , first? [Hint: names are CamelToed!] #possible-fix #available

    end=#
end
paths = ["functions", "functions/Generators", "Distributions", "NeuralNet/Networks", "src/DataStructures/Arrays/Array.jl", "tests"]
#TODO: Iterate all locations in paths
#
try

    #include(safe_path(".")) #why refeing to thyself?
    include(safe_path("functions"))
    include(safe_path("functions/Generators"))
    include(safe_path("Distributions"))
    include(safe_path("NeuralNet/Networks"))
    include(safe_path("tests"))
    include(safe_path("src/constants.jl"))
    include(safe_path("src/DataStructures/Arrays/Array.jl"))

catch
    try
        strings = split(path, "/")
        n = length(strings)
        for i in enumerate(n)

            strings[1] = Upper(strings[1])#try Uppercasing the first letter
            path = joinpath(strings...)
        end
    catch

        println("path does not exist")
    end
end

end
=#
"""

DEPOT_PATH

  A stack of "depot" locations where the package manager, as well as Julia's code loading mechanisms, look for package registries, installed packages, named
  environments, repo clones, cached compiled package images, and configuration files. By default it includes:

    1. ~/.julia where ~ is the user home as appropriate on the system;

    2. an architecture-specific shared system directory, e.g. /usr/local/share/julia;

    3. an architecture-independent shared system directory, e.g. /usr/share/julia.

  So DEPOT_PATH might be:

  [joinpath(homedir(), ".julia"), "/usr/local/share/julia", "/usr/share/julia"]

  The first entry is the "user depot" and should be writable by and owned by the current user. The user depot is where: registries are cloned, new package
  versions are installed, named environments are created and updated, package repos are cloned, newly compiled package image files are saved, log files are
  written, development packages are checked out by default, and global configuration data is saved. Later entries in the depot path are treated as read-only
  and are appropriate for registries, packages, etc. installed and managed by system administrators.

  DEPOT_PATH is populated based on the JULIA_DEPOT_PATH environment variable if set.

  DEPOT_PATH contents
  =====================

  Each entry in DEPOT_PATH is a path to a directory which contains subdirectories used by Julia for various purposes. Here is an overview of some of the
  subdirectories that may exist in a depot:

    •  clones: Contains full clones of package repos. Maintained by Pkg.jl and used as a cache.

    •  compiled: Contains precompiled *.ji files for packages. Maintained by Julia.

    •  dev: Default directory for Pkg.develop. Maintained by Pkg.jl and the user.  <--- ok

    •  environments: Default package environments. For instance the global environment for a specific julia version. Maintained by Pkg.jl.

    •  logs: Contains logs of Pkg and REPL operations. Maintained by Pkg.jl and Julia.

    •  packages: Contains packages, some of which were explicitly installed and some which are implicit dependencies. Maintained by Pkg.jl.

    •  registries: Contains package registries. By default only General. Maintained by Pkg.jl.

  See also: JULIA_DEPOT_PATH, and Code Loading.
"""

depots() = Base.DEPOT_PATH

logdir(depot = depots1()) = joinpath(depot, "logs")
devdir(depot = depots1()) = get(ENV, "JULIA_PKG_DEVDIR", joinpath(depot, "dev"))
envdir(depot = depots1()) = joinpath(depot, "environments")
const UPDATED_REGISTRY_THIS_SESSION = Ref(false)
const OFFLINE_MODE = Ref(false)

function depots1()
    d = depots()
    isempty(d) && Pkg.Types.pkgerror("no depots found in DEPOT_PATH")
    return d[1]
end

function pkg_server()
    server = get(ENV, "JULIA_PKG_SERVER", "https://pkg.julialang.org")
    isempty(server) && return nothing
    startswith(server, r"\w+://") || (server = "https://$server")
    return rstrip(server, '/')
end

#--- For globally overriding in e.g. tests
const DEFAULT_IO = Ref{Union{IO,Nothing}}(nothing)
stderr_f() = something(DEFAULT_IO[], stderr)
stdout_f() = something(DEFAULT_IO[], stdout)
const PREV_ENV_PATH = Ref{String}("") #a goo rule of thumb

can_fancyprint(io::IO) = (io isa Base.TTY) && (get(ENV, "CI", nothing) != "true") # fancyprint: <interesting>

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
function update_y!(
  y,
  nbset,  # tuple with the following elements.
    # nodes: array of Strings containing independent node names.
    # num_nodes: number of nodes.
    # branches: array of struct Branch.
    # num_branches: sum of branch phases.
    # num_dyns: number of dynamic branches.
  h,  # time-step size.
  x1  # value of the unknown vector x at the previous time step.
  )  #=

updates the vector "y", which is given as the first argument, for a transient simulation. This update is done only once for each time step before entering a Newton-Raphson iteration. In this function, only source and dynamic equations are updated. Nonlinear equations are updated in the function update_Fy!().

Author: Taku Noda
Started on Aug. 29, 2020
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

( nodes, num_nodes, branches, num_branches, num_dyns ) = nbset

noff_v = num_nodes + num_branches
noff_d = num_nodes + 2*num_branches
h2 = 0.5*h

for branch in branches
  if branch.type == "E"
    vec = eval( branch.val )
    y[ noff_v + branch.p_ns ] = vec[ 1 ]
  elseif branch.type == "J"
    vec = eval( branch.val )
    y[ num_nodes + branch.p_ns ] = vec[ 1 ]
  elseif branch.type == "L"
    nph = branch.nph
    range = 0:( nph - 1 )
    v_range = ( noff_v + branch.p_vars ) .+ range
    d_range = ( noff_d + branch.p_de ) .+ range
    y[ d_range ] = h2*x1[ v_range ] + x1[ d_range ]
  elseif branch.type == "C"
    nph = branch.nph
    range = 0:( nph - 1 )
    i_range = ( num_nodes + branch.p_vars ) .+ range
    d_range = ( noff_d + branch.p_de ) .+ range
    y[ d_range ] = h2*x1[ i_range ] + x1[ d_range ]
  elseif branch.type == "NL"
    dpos = noff_d + branch.p_de
    y[ dpos ] = h2*x1[ noff_v + branch.p_vars ] + x1[ dpos ]
  elseif branch.type == "NC"
    dpos = noff_d + branch.p_de
    y[ dpos ] = h2*x1[ num_nodes + branch.p_vars ] + x1[ dpos ]
  end
end

return nothing

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
# # Example 3 : Tiled Domain
#
# This example first reads in Global Ocean variables, generated by the [MITgcm](https://mitgcm.readthedocs.io/en/latest/index.html), which are split into a collection of subdomain arrays ( _tiles_ ) using `MeshArrays.jl`. It then writes them to a collection of `NetCDF` files ( _nctiles_ ) using `NCTiles.jl`.

# ### Packages & Helper Functions
#

using NCTiles;

p=dirname(pathof(NCTiles))
fil = joinpath(p, "../examples/helper_functions.jl")
include(fil);

# ### File Paths & I/O Back-End
#

inputs=NCTiles.NCTILES_TESTCASES
NCTiles.ensure_testcases_installed()

outputs = joinpath(tempdir(),"NCTILES_TESTCASES_OUTPUT/ex3/")
if ~ispath(outputs); mkpath(outputs); end

nc=NCTiles.NCDatasets # I/O Back-End

# ## Process Global Ocean Variables
#
# Here we process a two-dimensional field (`ETAN`), two three-dimensional tracer fields (`THETA`, `SALT`), and the three components of a vector field (`UVELMASS`, `VVELMASS` and `WVELMASS`). In each case, `flds` denotes the lazy representation of the processing chain (incl. all needed metadata) which the `write` function instantiates (i.e. outputs to files).
#
# _Note: on a `C-grid` these components are staggered in space._

writedir=outputs
readme = readlines(inputs)

# ### 2D example
(flds,savename,readme)=prep_nctiles_native(inputs,"state_2d_set1","ETAN",Float32)
write(flds,savename,README=readme);

# ### 3D examples
(flds,savename,readme)=prep_nctiles_native(inputs,"state_3d_set1","THETA",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"state_3d_set1","SALT",Float32);
write(flds,savename,README=readme);

# ### 3D vector field examples
(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","UVELMASS",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","VVELMASS",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","WVELMASS",Float32);
write(flds,savename,README=readme);


text:
using Printf

function dsdemo3()
    println("Modulator realization and scaling")
    println()

    fig = figure(1)
    fig.clf()
    fig.canvas.set_window_title("NTF")

    order = 5
    OSR = 42
    opt = 1
    H = synthesizeNTF(order, OSR, opt)
    a, g, b, c = realizeNTF(H)
    b = [b[1]; zeros(length(b)-1)]  # Use a single feed-in for the input

    plotPZ(H)

    println("Unscaled modulator")
    print("   DAC feedback coefficients = ")
    for i = 1:order
        @printf(" %.6f", a[i])
    end
    print("\n   resonator feedback coefficients = ")
    for i = 1:order÷2
        @printf(" %.6f", g[i])
    end
    println()
    println()

    println("Calculate the state maxima.")
    ABCD = stuffABCD(a, g, b, c)
    u = range(0, 0.6, length=30)
    N = 10000
    T = ones(N)
    maxima = zeros(order, length(u))
    for i = 1:length(u)
        ui = u[i]
        v, xn, xmax, = simulateDSM(ui * T, ABCD)
        maxima[:,i] = xmax[:]
        if any(xmax .> 1e2)
            umax = ui
            u = u[1:i]
            maxima = maxima[:,1:i]
            break
        end
    end

    fig = figure(2)
    fig.clf()
    fig.canvas.set_window_title("Simulated State Maxima")

    for i = 1:order
        semilogy(u, maxima[i,:], linestyle="--", linewidth=1, marker="o", markerfacecolor="none", markeredgewidth=0.5)
    end
    grid()
    grid(which="minor", linestyle="dotted")
    axis([0, 0.6, 1e-4, 10])
    xlabel("DC input")

    println("paused")
    readline()

    println("Calculate the scaled coefficients.")

    scaleABCD(ABCD, 2, 0, 1, NaN, NaN, 10000)
    ABCDs, umax, = scaleABCD(ABCD, 2, 0, 1, NaN, NaN, 10000)
    as, gs, bs, cs = mapABCD(ABCDs)

    @printf("Scaled modulator, umax=%.2f\n", umax)
    print("   DAC feedback coefficients = ")
    for i = 1:order
        @printf(" %.6f", as[i])
    end
    print("\n   resonator feedback coefficients = ")
    for i = 1:order÷2
        @printf(" %.6f", gs[i])
    end
    print("\n   interstage coefficients = ")
    for i = 1:order
        @printf(" %.6f", cs[i])
    end
    print("\n   feed-in coefficients = ")
    for i = 1:order
        @printf(" %.6f", bs[i])
    end
    println()
    println()

    println("Calculate the state maxima.")
    u = range(0, umax, length=30)
    N = 10000
    T = ones(N)
    maxima = zeros(order, length(u))
    for i = 1:length(u)
        ui = u[i]
        v, xn, xmax, = simulateDSM(ui * T, ABCDs)
        maxima[:,i] = xmax[:]
        if any(xmax .> 1e2)
            umax = ui
            u = u[1:i]
            maxima = maxima[:,1:i]
            break
        end
    end

    fig = figure(2)
    fig.clf()

    for i = 1:order
        semilogy(u, maxima[i,:], linestyle="--", linewidth=1, marker="o", markerfacecolor="none", markeredgewidth=0.5)
    end
    grid()
    xlabel("DC input")
    grid(which="minor", linestyle="dotted")
    axis([0, 0.6, 4e-2, 4])

    fig.canvas.manager.window.raise_()

    return nothing
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
"""
    CharXRay

Represents a specific known characteristic X-ray in a specific element.
"""
struct CharXRay
    z::Int
    transition::Transition
    function CharXRay(z::Int, transition::Transition)
        @assert charactericXRayAvailable(
            z,
            transition.innershell.index,
            transition.outershell.index
        ) "$z $transition unknown!"
        "$(symbol(element(z))) $(transition) does not occur."
        return new(z, transition)
    end
end

Base.hash(cxr::CharXRay, h::UInt)::UInt = hash(cxr.z, hash(cxr.transition, h))

Base.isequal(cxr1::CharXRay, cxr2::CharXRay) =
    isequal(cxr1.z, cxr2.z) && isequal(cxr1.transition, cxr2.transition)

Base.isless(cxr1::CharXRay, cxr2::CharXRay) =
    isequal(cxr1.z, cxr2.z) ? isless(cxr1.transition, cxr2.transition) : cxr1.z < cxr2.z

characteristic(elm::Element, tr::Transition) = CharXRay(z(elm), tr)

Base.show(io::IO, cxr::CharXRay) = print(io, element(cxr.z).symbol, " ", cxr.transition)

"""
    inner(cxr::CharXRay)

Return the inner AtomicSubShell associated with the specified CharXRay.
"""
inner(cxr::CharXRay) = AtomicSubShell(cxr.z, cxr.transition.innershell)

"""
    outer(cxr::CharXRay)

Return the outer AtomicSubShell associated with the specified CharXRay.
"""
outer(cxr::CharXRay) = AtomicSubShell(cxr.z, cxr.transition.outershell)

"""
    transition(cxr::CharXRay)

Return the Transition structure associated with the specified CharXRay.
"""
transition(cxr::CharXRay) = cxr.transition

"""
    shell(cxr::CharXRay)

Return the shell, Shell(1), Shell(2), ... associated with the inner shell.
"""
shell(cxr::CharXRay) = shell(cxr.transition)

"""
    element(cxr::CharXRay)

Return the element for this CharXRay.
"""
element(cxr::CharXRay) = PeriodicTable.elements[cxr.z]

"""
    normweight(elm::Element, tr::Transition, overvoltage = 4.0)::Float64

Return the nominal line strength for the specified transition in the specified element.
The strength differs from the weight in that the weight is normalized to the most intense line in the shell.
"""
normweight(elm::Element, tr::Transition, overvoltage = 4.0) =
    has(elm, tr) ? normweight(characteristic(elm, tr), overvoltage) : 0.0

ν(cxr::CharXRay) = energy(cxr) / plancksConstant
ω(cxr::CharXRay) = 2π * ν(cxr)

"""
    λ(cxr::CharXRay)

X-ray wavelength in cm.
"""
λ(cxr::CharXRay) = hc / energy(cxr)
λ(energy::Real) = hc / energy
"""
    wavenumber(cxr::CharXRay)

X-ray wavenumber in cm¯¹.
"""
wavenumber(cxr::CharXRay) = 1.0 / λ(cxr)
wavenumber(energy::Real) = 1.0 / λ(energy)

"""
    weight(cxr::CharXRay)

The line weight of the specified characteristic X-ray relative to the other lines from the same element in the
same shell.  The most intense line is normalized to unity.
"""
function weight(cxr::CharXRay, overvoltage = 4.0)::Float64
    return strength(cxr) / maxweight(inner(cxr))
end

const __normWeights = Dict{AtomicSubShell, Float64}()

"""
    normweight(cxr::CharXRay)

The line weight of the specified characteristic X-ray with the sum of the
weights in a subshell equals to unity.
"""
function normweight(cxr::CharXRay, overvoltage = 4.0)::Float64
    ass = inner(cxr)
    if !haskey(__normWeights, ass)
        safeSS(elm, tr) = has(elm, tr) ? strength(elm, tr) : 0.0
        __normWeights[ass] = sum(safeSS(element(ass), tr2) for tr2 in transitionsbysubshell[subshell(ass)])
    end
    return strength(cxr) / __normWeights[ass]
end

"""
    brightest(elm::Element, shell)

Return the brightest transition among the shell of transitions for the
specified element.  (group="K"|"Ka"|"Kb"|"L" etc)
"""
brightest(elm::Element, transitions) = brightest(characteristic(elm, transitions))

"""
    characteristic(elm::Element, iter::Tuple{Vararg{Transition}}, minweight=0.0, maxE=1.0e6)

"""
characteristic(
    elm::Element,
    iter::Tuple{Vararg{Transition}},
    minweight = 0.0,
    maxE = 1.0e6,
) = characteristic(
    elm,
    iter,
    cxr -> (weight(cxr) > minweight) && (energy(inner(cxr)) <= maxE),
)
characteristic(
    elm::Element,
    iter::AbstractVector{Transition},
    minweight = 0.0,
    maxE = 1.0e6,
) = characteristic(
    elm,
    iter,
    cxr -> (weight(cxr) > minweight) && (energy(inner(cxr)) <= maxE),
)

"""
    characteristic(elm::Element, iter::Tuple{Vararg{Transition}}, filterfunc::Function)
    characteristic(elm::Element, iter::AbstractVector{Transition}, filterfunc::Function)
    characteristic(elm::Element, iter::AbstractVector{Transition}, minweight = 0.0, maxE = 1.0e6)
    characteristic(elm::Element, iter::Tuple{Vararg{Transition}}, minweight = 0.0, maxE = 1.0e6)
    characteristic(ass::AtomicSubShell)

The collection of available CharXRay for the specified `Element` or `AtomicSubShell`.  `filterfunc(...)`
    * ` maxE` is compared to the edge energy.
    * `minWeight` is compared to the weight

Example:

    characteristic(n"Fe",ltransitions,0.01)
    characteristic(n"Fe",ltransitions,cxr->energy(cxr)>700.0)
    characteristic(n"Fe L3")
"""
function characteristic(
    elm::Element,
    iter::Tuple{Vararg{Transition}},
    filterfunc::Function,
)::Vector{CharXRay}
    res = CharXRay[]
    for tr in iter
        if has(elm, tr)
            cxr=characteristic(elm, tr)
            filterfunc(cxr) && push!(res, cxr)
        end
    end
    res
end

characteristic(
    elm::Element,
    iter::AbstractVector{Transition},
    filterfunc::Function,
)::Vector{CharXRay} = map(
    tr -> characteristic(elm, tr),
    filter(tr -> has(elm, tr) && filterfunc(characteristic(elm, tr)), collect(iter)),
)

characteristic(ass::AtomicSubShell) =
    characteristic(element(ass), filter(tr -> inner(tr) == ass.subshell, alltransitions))

"""
    splitbyshell(cxrs)

Splits a collection of CharXRay into a dictionary where the key is the inner
AtomicSubShell and the values are a vector of CharXRay for that inner shell.
"""
function splitbyshell(cxrs)
    res = Dict{AtomicSubShell,Vector{CharXRay}}()
    for cxr in cxrs
        shell = inner(cxr)
        if haskey(res, shell)
            push!(res[shell], cxr)
        else
            res[shell] = [cxr]
        end
    end
    return res
end

brightest(cxrs::Vector{CharXRay})::CharXRay = cxrs[findmax(weight.(cxrs))[2]]

"""
    name(cxrs::AbstractVector{CharXRay}, byfamily=false)

An abbeviated name for a collection of CharXRay.
"""
function name(cxrs::AbstractVector{CharXRay}, byfamily::Bool=false)::String
    if byfamily
        if all(transition(cxr) in kalpha for cxr in cxrs)
            "Kα"
        elseif all(transition(cxr) in kbeta for cxr in cxrs)
            "Kβ"
        elseif all(transition(cxr) in ktransitions for cxr in cxrs)
            "K"
        elseif all(transition(cxr) in lalpha for cxr in cxrs)
            "Lα"
        elseif all(transition(cxr) in lbeta for cxr in cxrs)
            "Lβ"
        elseif all(transition(cxr) in ltransitions for cxr in cxrs)
            "L"
        elseif all(transition(cxr) in malpha for cxr in cxrs)
            "Mα"
        elseif all(transition(cxr) in mbeta for cxr in cxrs)
            "Mβ"
        elseif all(transition(cxr) in mtransitions for cxr in cxrs)
            "M"
        elseif all(transition(cxr) in ntransitions for cxr in cxrs)
            "N"
        else
            "Unknown"
        end
    else
        res = []
        elms = Set{Element}(element.(cxrs))
        for elm in elms
            for sh in Set(shell.(filter(cxr -> element(cxr) == elm, cxrs)))
                fc = filter(cxr -> (element(cxr) == elm) && (shell(cxr) == sh), cxrs)
                br, cx = brightest(fc), length(fc)
                other = cx > 2 ? "others" : "other"
                push!(res, cx > 1 ? "$(br) + $(cx-1) $(other)" : "$(br)")
            end
        end
        return join(res, ", ")
    end
end

function Base.show(io::IO, cxrs::AbstractVector{CharXRay})
    print(io, name(cxrs))
end

function NeXLUncertainties.asa(::Type{DataFrame}, cxrs::AbstractVector{CharXRay})
    cc = sort(cxrs)
    return DataFrame(
        XRay = cc,
        Inner = inner.(cc),
        Outer = outer.(cc),
        Energy = energy.(cc),
        Strength = strength.(cc),
        Weight = weight.(cc),
    )
end


text:
# SDE.jl

"""
    abstract type AbstactSDE <: Any

subtypes: SDE and SDEWithConstantVariance.

Both subtypes implement a stochastic differential equation
dX_t=θ(X_t)dt + σ(X_t)dW_t.
"""
abstract type AbstractSDE end

"""
    struct SDE{S<:Function, T<:Function} <: AbstractSDE
      b::S
      σ::T
      beginvalue::Float64
      endtime::Float64
      dt::Float64
    end

Implements a stochastic differential equation dX_t=b(X_t)dt+σ(X_t)dW_t on time
interval [0,endtime], with W_t a Brownian motion. Both b and σ are functions
R→R. Beginvalue X_0=beginvalue and dt>0 is the precision with which we
discretise.

See also: SDEWithConstantVariance.

# Example
```julia
using PyPlot
clf()
sde = SDE(sin, x -> 2 + sin(x), 0, 10, 0.001)
X = rand(sde)
#
plot(X)
```
"""
struct SDE{S<:Function, T<:Function} <: AbstractSDE
  b::S
  σ::T
  beginvalue::Float64
  endtime::Float64
  dt::Float64

  function SDE(b::S, σ::T, beginvalue, endtime, dt) where {S<:Function, T<:Function}
    dt > 0 || error("The time discretisation (dt) should be positive")
    endtime > 0 || error("The endtime should be positive")
    new{S, T}(b, σ, beginvalue, endtime, dt)
  end
end

"""
    struct SDEWithConstantVariance{S<:Function} <: AbstractSDE
        b::S
        σ::Float64
        beginvalue::Float64
        endtime::Float64
        dt::Float64
    end

Implements the stochastic differential equation dX_t=b(X_t)dt+ σdW_t, with
b:R→R a function and σ a real number. Is a special case of SDE type with field σ
set equal to the lambda expression x->σ, but is faster.

# Example
```julia
using PyPlot
clf()
#
sde = SDEWithConstantVariance(cos, 1, 0, 10, 0.001)
X = rand(sde)
plot(X)
```
"""
struct SDEWithConstantVariance{S<:Function} <: AbstractSDE
    b::S
    σ::Float64
    beginvalue::Float64
    endtime::Float64
    dt::Float64

    function SDEWithConstantVariance(b::T, σ, beginvalue, endtime, dt) where
            T<:Function
        dt > 0 || error("The time discretisation (dt) should be positive")
        endtime > 0 || error("The endtime should be positive")
        new{T}(b, σ, beginvalue, endtime, dt)
    end
end

# Extends Distributions.rand for SDEWithConstantVariance types.
"""
    rand(sde::SDEWithConstantVariance)
    rand(sde::SDE)

Returns a SamplePath object which represents a sample path from an
AbstractSDE subtype sde. From time 0.0 to time sde.endtime, discretised with
precision sde.dt.

# Examples

```julia
using PyPlot
clf()
sde = SDEWithConstantVariance(sin, 1.0, 0.0, 10.0, 0.01)
X = rand(sde)
plot(X)
```

```julia
using PyPlot
clf()
sde = SDE(x->0.5*x, identity, 1, 1, 0.001)
X = rand(sde)
plot(X)
```
"""
function rand(sde::SDEWithConstantVariance)::SamplePath
  timeinterval = 0.0:sde.dt:sde.endtime
  lengthoftimeinterval = length(timeinterval)
  # Brownian motion increments multiplied with variance.
  increments = sde.σ * rand(Normal(0, sqrt(sde.dt)), lengthoftimeinterval)
  samplevalues = Array{Float64}(lengthoftimeinterval)
  samplevalues[1] = sde.beginvalue
  prevXval = sde.beginvalue
  for k in 2:lengthoftimeinterval
    samplevalues[k] = prevXval + sde.b(prevXval)*sde.dt + increments[k]
    prevXval = samplevalues[k]
  end
  return SamplePath(timeinterval, samplevalues)
end

# Extends Distributions.rand for SDE types.
function rand(sde::SDE)::SamplePath
  timeinterval = 0.0:sde.dt:sde.endtime
  lengthoftimeinterval = length(timeinterval)
  BMincrements = rand(Normal(0, sqrt(sde.dt)), lengthoftimeinterval)
  samplevalues = Array{Float64}(lengthoftimeinterval)
  samplevalues[1] = sde.beginvalue
  prevXval = sde.beginvalue
  for k in 2:lengthoftimeinterval
    samplevalues[k] = prevXval + sde.b(prevXval)*sde.dt +
      sde.σ(prevXval)*BMincrements[k]
    prevXval = samplevalues[k]
  end
  return SamplePath(timeinterval, samplevalues)
end