text:
include("deflump.jl")

filename = "methanecomplete.kpp"

f = open(filename)
lines = [replace(ln,r"[;\h]+","") for ln in split(replace(readstring(f),r":\h*\n",":"),"\n")]
eqns = [split(eq,r"[\}=:]") for eq in filter(x->ismatch(r"\{[\d\.]+\}", x), lines)]

specs = Set(vcat(pmap(getspecs,eqns)...))
nspecs = length(specs)

#rate number dictionary
rated = Set([r[4] for r in eqns])
r_rated = Dict(zip(UInt32(1):UInt32(length(rated)),rated))
rated = Dict(zip(rated,UInt32(1):UInt32(length(rated))))

#spec dicionary
specd = Dict(zip(specs,UInt32(1):UInt32(length(specs))))
rspec = Dict(zip(UInt32(1):UInt32(length(specs)),specs))

@eval @everywhere rated=$rated
@eval @everywhere specd=$specd

#each eq numbered
eqnum = pmap(numerate,eqns)
@eval @everywhere eqnum=$eqnum

###############################

lumplist = [["CH3O2","CH3O2NO2"],["HCHO","CH3O"]]
#lumplist = [["CH3O2","CH3O2NO2"]]
#lumplist = [["CL"]]#,["HCHO","CH3O"]]

pspecs = [specd[x] for x in collect(Base.flatten(lumplist))]
parr = pmap(production_array,pspecs)
parr = Dict(zip(pspecs,parr))

lumpdict = []
lumpspec = []
  ppair = []

counter = 1
for pairing in enumerate(lumplist)
  for j in pairing[2]
    push!(lumpdict, [j,pairing[1]] )
    push!(lumpspec, [j,counter] )
    counter +=1
  end
end

  @eval @everywhere lumpspec =  Dict(lumpspec)
lumpdict = Dict(lumpdict)

for element in enumerate(lumplist)

  @eval @everywhere pnum = $element[1]
  pairing = element[2]
  p(pairing)

  pair = [specd[i] for i in pairing]

  #a-> b :r  if l->l ignore . Therefore ignore and cancel in normal eqnum
  p(pair)

  push!(ppair, "\n! LMP$(pnum): "*join(["$(lumpspec[i]) -> $i" for i in pairing], ", ")*"\n")

  for n in 1:length(pair)
    i = pair[n]
    temp = "0"
    for j in parr[i]
      s = ""
      self = false
      for k in j[1][1]
        if k != "1"
          s = s * k *"*"
        end
      end
      for k in j[2][1]
        #if k in pair
          #self = true
          #break
        #end
        try
          s=s * "C(ind_LMP$(lumpdict[rspec[k]]))*LC_$(lumpspec[rspec[k]])_*"
        catch
          s=s * "C(ind_"* rspec[k] * ")*"
        end
      end
      s=s * r_rated[j[3]]
      #if self == false
        temp = temp*"+"*s*""
      #end
    end

    thisspec = rspec[i]
    thisno = lumpspec[thisspec]

    #temp = replace(tostr(sympy.N(sympy.expand(replace(temp,r"(\d)D([\+\-\d])",s"\1E\2")),3)),"+","&\n      +")
    temp = replace(temp,"+","&\n      +")

    temp = replace(temp,r"EXP\(([-\d\.]+)/TEMP\)",s"ET(\1)")
    temp = replace(temp,r"LC_(\d)_",s"LC(\1)")
    push!(ppair, "PA($(thisno))="*temp)

    push!(ppair, "LC($(thisno)) =  PA($(thisno))/("*join(["PA($(lumpspec[rspec[e]]))" for e in pair],"+")*")")

  end

  @eval @everywhere pgroup=$pairing

  eqns = pmap(replace_match,eqns)

end

include("writelump.jl")

x= 1

max_stars_repo_name: wolfiex/MechLump

max_issues_repo_name: wolfiex/MechLump

max_forks_repo_name: wolfiex/MechLump


text:
module MahjongTiles

# define the suits and all the tiles
include("tiles.jl")

# from the tiles, create decks of tiles for games
include("decks.jl")

# check "hands"
include("hands.jl")

end

ext: jl


text:
x = 3

y = 2x

typeof(y)

f(x) = 2 + x

f

f(10)

function g(x, y)
    z = x + y
    return z^2
end

g(1, 2)

let s = 0
    for i in 1:10
        s += i    # Equivalent to s = s + i
    end
    s
end

typeof(1:10)

function mysum(n)
    s = 0
    for i in 1:n
        s += i
    end
    return s
end

mysum(100)

a = 3

a < 5

if a < 5
    "small"
else
    "big"
end

v = [1, 2, 3]

typeof(v)

v[2]

v[2] = 10

v2 = [i^2 for i in 1:10]

M = [1 2
     3 4]

typeof(M)

zeros(5, 5)

zeros(Int, 4, 5)

[i + j for i in 1:5, j in 1:6]


text:
"""
A linear Gaussian CPD, always returns a Normal

	Assumes that target and all parents can be converted to Float64 (ie, are numeric)

	P(x|parents(x)) = Normal(μ=a×parents(x) + b, σ)
"""
mutable struct LinearGaussianCPD <: CPD{Normal}
    target::NodeName
    parents::NodeNames
    a::Vector{Float64}
    b::Float64
    σ::Float64
end
LinearGaussianCPD(target::NodeName, μ::Float64, σ::Float64) = LinearGaussianCPD(target, NodeName[], Float64[], μ, σ)

name(cpd::LinearGaussianCPD) = cpd.target
parents(cpd::LinearGaussianCPD) = cpd.parents
nparams(cpd::LinearGaussianCPD) = length(cpd.a) + 2
function (cpd::LinearGaussianCPD)(a::Assignment)

    # compute A⋅v + b
    μ = cpd.b
    for (i, p) in enumerate(cpd.parents)
        μ += a[p]*cpd.a[i]
    end

    Normal(μ, cpd.σ)
end
(cpd::LinearGaussianCPD)() = (cpd)(Assignment()) # cpd()
(cpd::LinearGaussianCPD)(pair::Pair{NodeName}...) = (cpd)(Assignment(pair)) # cpd(:A=>1)

function Distributions.fit(::Type{LinearGaussianCPD},
    data::DataFrame,
    target::NodeName;
    min_stdev::Float64=0.0, # an optional minimum on the standard deviation
    )

    # no parents

    arr = data[!,target]
    eltype(arr) <: Real || error("fit LinearGaussianCPD requrires target to be numeric")

    μ = convert(Float64, mean(arr))
    σ = convert(Float64, stdm(arr, μ))
    σ = max(σ, min_stdev)

    LinearGaussianCPD(target, NodeName[], Float64[], μ, σ)
end
function Distributions.fit(::Type{LinearGaussianCPD},
    data::DataFrame,
    target::NodeName,
    parents::NodeNames;
    min_stdev::Float64=0.0, # an optional minimum on the standard deviation
    )

    if isempty(parents)
        return fit(LinearGaussianCPD, data, target, min_stdev=min_stdev)
    end

    # ---------------------
    # pull parental dataset
    # 1st row is all of the data for the 1st parent
    # 2nd row is all of the data for the 2nd parent, etc.

    nparents = length(parents)
    X = Array{Float64}(undef, nrow(data), nparents+1)
    for (i,p) in enumerate(parents)
        arr = data[!,p]
    	for j in 1 : nrow(data)
            X[j,i] = convert(Float64, arr[j])
    	end
    end
    X[:,end] .= 1.0

    y = convert(Vector{Float64}, data[!,target])

    # --------------------
    # solve the regression problem
    #   β = (XᵀX)⁻¹Xᵀy
    #
    #     X is the [nsamples × nparents+1] data matrix
    #     where the last column is 1.0
    #
    #     y is the [nsamples] vector of target values
    #
    # NOTE: this will fail if X is not full rank

    β = (X'*X)\(X'*y)

    a = β[1:nparents]
    b = β[end]
    σ = max(std(y), min_stdev)

    LinearGaussianCPD(target, parents, a, b, σ)
end


text:
@defcomp welfare begin
    regions = Index()

    PERIODU = Variable(index=[time, regions]) # One period utility function
    CEMUTOTPER = Variable(index=[time, regions]) # Period utility
    REGCUMCEMUTOTPER = Variable(index=[time, regions]) # Cumulative period utility
    REGUTILITY = Variable(index=[regions]) # Welfare Function
    UTILITY = Variable()

    CPC = Parameter(index=[time, regions])
    l = Parameter(index=[time, regions]) # Level of population and labor
    elasmu = Parameter(index=[regions]) # Elasticity of marginal utility of consumption
    rr = Parameter(index=[time, regions]) # Average utility social discount rate
    scale1 = Parameter(index=[regions]) # Multiplicative scaling coefficient
    scale2 = Parameter(index=[regions]) # Additive scaling coefficient
    alpha = Parameter(index=[time, regions])

    function run_timestep(p, v, d, t)

        #Define function for PERIODU #NEED TO ADD IF STATEMENT LIKE IN JUMP MODEL OR IS THAT ONLY ISSUES WHEN ELASMU = 1.0?
        for r in d.regions
            if p.elasmu[r]==1.
                v.PERIODU[t,r] = log(p.CPC[t,r]) * p.alpha[t,r]
            else
                v.PERIODU[t,r] = ((1. / (1. - p.elasmu[r])) * (p.CPC[t,r])^(1. - p.elasmu[r]) + 1.) * p.alpha[t,r]
            end
        end

        #Define function for CEMUTOTPER
        for r in d.regions
            if t.t != 60
                v.CEMUTOTPER[t,r] = v.PERIODU[t,r] * p.l[t,r] * p.rr[t,r]
            else
                v.CEMUTOTPER[t,r] = v.PERIODU[t,r] * p.l[t,r] * p.rr[t,r] / (1. - ((p.rr[t-1,r] / (1. + 0.015)^10) / p.rr[t-1,r]))
            end
        end

        #Define function for REGCUMCEMUTOTPER
        for r in d.regions
            if is_first(t)
                v.REGCUMCEMUTOTPER[t,r] = v.CEMUTOTPER[t,r]
            else
                v.REGCUMCEMUTOTPER[t,r] = v.REGCUMCEMUTOTPER[t-1, r] + v.CEMUTOTPER[t,r]
            end
        end

        if t.t == 60
            #Define function for REGUTILITY
            for r in d.regions
                v.REGUTILITY[r] = 10 * p.scale1[r] * v.REGCUMCEMUTOTPER[t,r] + p.scale2[r]
            end
            #Define function for UTILITY
            v.UTILITY = sum(v.REGUTILITY[:])
        end
    end
end


text: this is just a test .jl file


text:
using Assignment07
using Test

@testset "Assignment07" begin

@testset "Using Strings" begin

    @testset "normalizeDNA" begin
        @test normalizeDNA("aatgn") == "AATGN"
        @test_throws Exception normalizeDNA("ZCA")
        @test_throws Exception normalizeDNA(42)
        c = normalizeDNA('C')
        @test c == "C"
        @test typeof(c) == String
    end # normalizeDNA

    @testset "composition" begin
        seq = rand(['A','T','G','C','N'], 20) |> join
        bc = composition(seq)
        @test bc isa Dict

        @test bc['A'] == count(x-> x == 'A', seq)
        @test bc['C'] == count(x-> x == 'C', seq)
        @test bc['G'] == count(x-> x == 'G', seq)
        @test bc['T'] == count(x-> x == 'T', seq)
        @test bc['N'] == count(x-> x == 'N', seq)

        bc = composition(lowercase(seq))

        @test bc['A'] == count(x-> x == 'A', seq)
        @test bc['C'] == count(x-> x == 'C', seq)
        @test bc['G'] == count(x-> x == 'G', seq)
        @test bc['T'] == count(x-> x == 'T', seq)
        @test bc['N'] == count(x-> x == 'N', seq)
    end # composition

    @testset "gc_content" begin
        @test gc_content("ANTG") == 0.25
        @test gc_content("cccggg") * 100 == 100.0
        @test gc_content("ATta") == 0.0
        @test_throws Exception gc_content("ATty")
    end # gc_content

    @testset "complement" begin
        @test complement("ATTAN") == "TAATN"
        @test complement("gcta") == "CGAT"
        @test complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception complement("ABC")
    end # complement

    @testset "reverse_complement" begin
        @test reverse_complement("ATTAN") == "NTAAT"
        @test reverse_complement("gcta") == "TAGC"
        @test reverse_complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception reverse_complement("ABC")
    end # reverse_complement

    @testset "parse_fasta" begin
        testpath = normpath(joinpath(@__DIR__, "..", "data"))
        genomes = joinpath(testpath, "cov2_genomes.fasta")
        ex1_path = joinpath(testpath, "ex1.fasta")
        ex2_path = joinpath(testpath, "ex2.fasta")

        ex1 = parse_fasta(ex1_path)
        @test ex1 isa Tuple
        @test all(x-> x isa String, ex1[1])
        @test all(x-> x isa String, ex1[2])

        @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
        @test ex1[2] == ["AATTATAGC", "CGCCCCCCAGTCGGATT"]

        @test_throws Exception parse_fasta(ex2_path)

        cov2 = parse_fasta(genomes)
        @test length(cov2[1]) == 8
        @test length(cov2[2]) == 8
    end #parse_fasta

end # strings

# @testset "Using BioSequences" begin

#     @testset "normalizeDNA" begin
#         @test normalizeDNA("aatgn") == dna"AATGN"
#         @test_throws Exception normalizeDNA("ZCA")
#         @test_throws Exception normalizeDNA(42)
#         c = normalizeDNA('C')
#         @test c == dna"c"
#         @test c isa LongSequence
#     end #  normalizeDNA

#     @testset "gc_content" begin
#         @test gc_content(dna"ANTG") == 0.25
#         @test gc_content(dna"cccggg") * 100 == 100.0
#         @test gc_content(dna"ATta") == 0.0
#     end #  composition

#     @testset "composition" begin
#         seq = rand(['A','T','G','C','N'], 20) |> join |> LongDNASeq
#         bc = composition(seq)

#         @test bc[DNA_A] == count(==(DNA_A), collect(seq))
#         @test bc[DNA_C] == count(==(DNA_C), collect(seq))
#         @test bc[DNA_G] == count(==(DNA_G), collect(seq))
#         @test bc[DNA_T] == count(==(DNA_T), collect(seq))
#         @test bc[DNA_N] == count(==(DNA_N), collect(seq))
#     end #  gc_content

#     @testset "complement" begin
#         @test complement(dna"ATTAN") == dna"TAATN"
#         @test complement(dna"gcta") == dna"CGAT"
#         @test complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  complement

#     @testset "reverse_complement" begin
#         @test reverse_complement(dna"ATTAN") == dna"NTAAT"
#         @test reverse_complement(dna"gcta") == dna"TAGC"
#         @test reverse_complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  reverse_complement

#     @testset "parse_fasta" begin
#         testpath = normpath(joinpath(@__DIR__, "..", "data"))
#         genomes = joinpath(testpath, "cov2_genomes.fasta")
#         ex1_path = joinpath(testpath, "ex1.fasta")
#         ex2_path = joinpath(testpath, "ex2.fasta")

#         ex1 = parse_fasta(ex1_path)
#         @test ex1 isa Tuple
#         @test all(x-> x isa String, ex1[1])
#         @test all(x-> x isa LongSequence, ex1[2])

#         @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
#         @test ex1[2] == [dna"AATTATAGC", dna"CGCCCCCCAGTCGGATT"]

#         @test_throws Exception parse_fasta(ex2_path)

#         cov2 = parse_fasta(genomes)
#         @test length(cov2[1]) == 8
#         @test length(cov2[2]) == 8
#     end # parse_fasta

# end # BioSequences

end # Assignment07


text:
include(joinpath(@__DIR__, "types.jl"));

struct NonInteractingEnergy <: Energy end

(::Energy)(::EAPChain) = error("Not yet implemented.");

@inline function (::NonInteractingEnergy)(chain::EAPChain)
  return (sum(chain.us) - dot(end_to_end(chain), [chain.Fx; 0.0; chain.Fz]));
end

struct InteractingEnergy <: Energy end

@inline function (::InteractingEnergy)(chain::EAPChain)
  return (sum(chain.us) + U_interaction(chain)
          - dot(end_to_end(chain), [chain.Fx; 0.0; chain.Fz]));
end

struct IsingEnergy <: Energy end

@inline function (::IsingEnergy)(chain::EAPChain)
  return (sum(chain.us) + U_Ising(chain)
          - dot(end_to_end(chain), [chain.Fx; 0.0; chain.Fz]));
end

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "mdbtools"
version = v"0.8.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/cyberemissary/mdbtools.git", "b753ff36a0f1d88ae8a300ed6712f4aa2ddb7d08")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd mdbtools/
autoreconf -if
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-man
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(exclude = [Windows(:i686), Windows(:x86_64)])

# The products that we will ensure are always built
products = [
    ExecutableProduct("mdb-count", :mdb_count),
    ExecutableProduct("mdb-ver", :mbd_ver),
    ExecutableProduct("mdb-prop", :mdb_prop),
    ExecutableProduct("mdb-sql", :mdb_sql),
    ExecutableProduct("mdb-tables", :mdb_tables),
    ExecutableProduct("mdb-schema", :mdb_schema),
    ExecutableProduct("mdb-parsecsv", :mdb_parsecsv),
    ExecutableProduct("mdb-header", :mdb_header),
    ExecutableProduct("mdb-export", :mdb_export),
    ExecutableProduct("mdb-hexdump", :mdb_hexdump),
    LibraryProduct("libmdbsql", :libmdbsql),
    LibraryProduct("libmdb", :libmdb),
    ExecutableProduct("mdb-array", :mdb_array)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Glib_jll", uuid="7746bdde-850d-59dc-9ae8-88ece973131d"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: kalmarek/Yggdrasil

max_issues_repo_name: kalmarek/Yggdrasil

max_forks_repo_name: kalmarek/Yggdrasil


text:
using POMDPs, Crux, Flux, POMDPGym, Random, Distributions, POMDPPolicies
Crux.set_function("isfailure", POMDPGym.isfailure)

struct MvNormal2
    μ
    logσ
    μ_gpu
    logσ_gpu
    d
    MvNormal2(μ, σ) = new(μ, Base.log.(σ), μ|>gpu, Base.log.(σ)|>gpu, MvNormal(cpu(μ), cpu(σ)))
end

function Base.rand(d::MvNormal2, sz::Int...)
    r = rand(d.d, sz...)
    if length(r) == 1
        return r[1]
    end
    r
end

function Distributions.logpdf(d::MvNormal2, x)
    if device(x)==gpu
        Crux.gaussian_logpdf(d.μ_gpu, d.logσ_gpu, x)
    else
        Crux.gaussian_logpdf(d.μ, d.logσ, x)
    end
end

# Define the disturbance distribution based on a normal distribution
px = MvNormal2([0f0, 0f0], [0.2f0, 0.2f0])

# Construct the MDP
function build_mdp()
    Random.seed!(1)
    rp() = rand(Uniform(-1, 1))
    rad() = abs(rand(Uniform(0.05, 0.15)))

    failures = Dict(Circle(Vec2(rp(),rp()), rad())=>-0.01f0 for i=1:1)
    rewards = Dict(failures..., Circle(Vec2(-0.5f0,0.5f0), 0.2f0)=>0.9f0)
    mdp = AdditiveAdversarialMDP(ContinuumWorldMDP(rewards=rewards, vel_thresh=1f0, discount=1.0), px)
end
mdp = build_mdp()
render(mdp)
S = state_space(mdp)

# construct the models for the protagonist
QSA() = ContinuousNetwork(Chain(Dense(S.dims[1] + 2, 128, relu), Dense(128, 128, relu), Dense(128, 1)))
A() = ContinuousNetwork(Chain(Dense(S.dims[1], 128, relu), Dense(128, 128, relu), Dense(128, 2, tanh)))
Protag() = ActorCritic(A(), DoubleNetwork(QSA(), QSA()))

# And the models for the antagonist
Pf() = ContinuousNetwork(Chain(Dense(S.dims[1] + 2, 128, relu), Dense(128, 128, relu), Dense(128, 1, sigmoid)))
function G()
    base = Chain(Dense(S.dims[1], 128, relu), Dense(128, 128, relu))
    mu = ContinuousNetwork(Chain(base..., Dense(128, 2)))
    logΣ = ContinuousNetwork(Chain(base..., Dense(128, 2)))
    GaussianPolicy2(mu, logΣ, true)
end
Antag() = ActorCritic(G(), Pf())

# Solve with a traditional RL approach
TD3_params() = (N=100000, S=S, buffer_size=Int(1e5), buffer_init=1000)

𝒮_td3 = TD3(;π=ActorCritic(A(), DoubleNetwork(QSA(), QSA())), TD3_params()...)
π_td3 = solve(𝒮_td3, mdp)

# Solver with adversarial approach
ARL_params() = (TD3_params()..., desired_AP_ratio=1.0)

# solve with IS
𝒮_isarl = ISARL_Continuous(;π=AdversarialPolicy(π_td3, Antag()) |> gpu, ARL_params()..., px)
π_isarl = solve(𝒮_isarl, mdp)

# solve with RARL
𝒮_rarl = Crux.RARL_Continuous(;π=AdversarialPolicy(Protag(), Protag(), GaussianNoiseExplorationPolicy(0.2f0)) |>gpu, ARL_params()...)
π_rarl = solve(𝒮_rarl, mdp)

Dvanilla = steps!(Sampler(mdp, π_td3, S=S, required_columns=[:fail]), Nsteps=10000)
sum(Dvanilla[:fail])

Disarl = steps!(Sampler(mdp, protagonist(π_isarl), S=S, required_columns=[:fail]), Nsteps=10000)
sum(Disarl[:fail])

Drarl = steps!(Sampler(mdp, protagonist(π_rarl), S=S, required_columns=[:fail]), Nsteps=10000)
sum(Drarl[:fail])

gif(mdp, protagonist(π_isarl), "isarl.gif", Neps=10)
gif(mdp, π_td3, "td3.gif", Neps=100)


text:
using .JSort
import YAML

function yamlvalue(yaml::Dict{String, Any}, key::String, default)
    if key ∉ keys(yaml)
        default
    else
        param = yaml[key]
    end
end

struct Parameters
    root::String
    readpath::String
    savepath::String
    #datafiles::Vector{String}
    isotope::String
    gates::Dict{Symbol, Gate}
    calibrator::Dict{Symbol, AbstractCalibrator}
    percent::Float64
    badchannels::Set{Int}
    qkinz::Vector{Poly{Float32}}
end

function Parameters(path::AbstractString)
    parameters::Dict{String,Any} = YAML.load_file(path)
    readkeys = Set{String}()
    function get(key::String, default)
        push!(readkeys, key)
        yamlvalue(parameters, key, default)
    end
    root     = splitdir(path)[1]
    savepath = get("save path", root)
    if startswith(savepath, "~")
        savepath = expanduser(savepath)
    else
        savepath = if savepath == root savepath else joinpath(root, savepath) end
    end

    readpath = get("read path", savepath)
    if startswith(readpath, "~")
        readpath = expanduser(readpath)
    else
        readpath = if readpath == savepath readpath else joinpath(root, readpath) end
    end
    isotope  = get("isotope", "")

    # Different calibration coefficients
    calibrationpaths = Dict{Symbol, Union{String, Vector{String}}}(
        :e  => get("coefficients e",     "coefficients_e.csv"),
        :Δe => get("coefficients de",    "coefficients_de.csv"),
        :γ  => get("coefficients gamma", "coefficients_g.csv"),
#        :γi => get("coefficients gamma inverse", "coefficients_g_i.csv"),
        :t  => get("coefficients time",  "coefficients_t.csv"),
        :tc => get("coefficients time correction", "coefficients_t_c.csv")
    )
    # Try to read from the coefficient paths, defaulting to identity polynomials
    calibrator = Dict{Symbol, AbstractCalibrator}(
        :e  => make(ParticleCalibrator, calibrationpaths[:e],  :e),
        :Δe => make(ParticleCalibrator, calibrationpaths[:Δe], :Δe),
        :γ  => make(GammaCalibrator,    calibrationpaths[:γ],  :eᵧ),
#        :γi => make(GammaCalibrator,    calibrationpaths[:γi], :eᵧ),
        :t  => make(TimeCalibrator,     calibrationpaths[:t],  :t),
        :tc => make(TimeCorrector,      calibrationpaths[:tc], :t)
    )

    gates = getgates(get("gates", nothing))

    # Read the datafile paths
    # datafiles = get("data files", nothing)
    # if isnothing(datafiles)
    #     @warn "No data files will be processed!"
    # else
    #     map!(x -> joinpath(readpath, x), datafiles, datafiles)
    # end

    percent = get("percent to read", 100)

    # Qkinz
    qkinz = get("qkinz", nothing) |> abspath |> parseqkinz

    # Bad channels
    badchannels = Set{Int}(i for i in get("channels to ignore", Set{Int}()))

    atleast1d(x::AbstractArray) = x
    atleast1d(x) = [x]
    for key in keys(parameters)
        if key ∉ readkeys
            @warn "The parameter $key was supplied but not used."
        end
        if occursin("coefficient", key) || occursin("file", key)
            paths = get(key, "") |> atleast1d
            for path in paths
              if !isfile(path)
                  @warn "Could not open $key"
              end
            end
        end
    end

    Parameters(path, readpath, savepath, isotope,
               gates, calibrator,
               percent, badchannels, qkinz)
end

function make(::Type{T}, path::AbstractString, var)::T where T<:AbstractCalibrator
    if isfile(path)
        read(path, T, var=var)
    else
        T(var=var)
    end
end

function make(::Type{T}, paths::AbstractArray{<:AbstractString}, var)::T where T<:AbstractCalibrator
    base = make(T, paths[1], var)
    if length(paths) > 1
      for path in paths[2:end]
          calib = make(T, path, var)
          combine!(base, calib)
      end
    end
    base
end

### PARSE GATES ###
getgates(::Nothing) = Dict(:ex => Gate(), :t => Gate())
function getgates(unparsedgates::AbstractDict)
    # Gates on the excitation energy
    ukeys = lowercase.(keys(unparsedgates))
    ekey =  ukeys ∩ ["ex", "excitation", "ex_int"]
    promptkey = ukeys ∩ ["prompt"]
    bgkey = ukeys ∩ ["bg", "background"]

    gates = Dict{Symbol,Gate}(:ex => Gate(), :background => Gate(),
                              :promt => Gate())

    if length(ekey) == 1
        for gate in unparsedgates[ekey[1]]
            gates[:ex] = gates[:ex] ∩ parsegate(gate)
        end
    end

    if length(promptkey) == 1
        gates[:prompt] = unparsedgates[promptkey[1]] |> parsegate
    end

    if length(bgkey) == 1
        gates[:background] = unparsedgates[bgkey[1]] |> parsegate
    end

    if length(ukeys) > 0 && length(ekey ∪ promptkey ∪ bgkey) == 0
        throw(IOError("Gates $ukeys not supported."))
    end
    gates
end

function parsegate(gate::AbstractDict)::Gate
    ukeys = gate |> keys .|> lowercase
    low = ukeys ∩ ["low", "l"]
    high = ukeys ∩ ["high", "h"]
    wrong = setdiff(ukeys, low ∪ high)
    if length(wrong) > 0
        throw(ArgumentError("Unsupported gate bounds: $wrong"))
    end

    if length(low) > 0 && length(high) > 0
        return Gate(low=gate[low[1]], high=gate[high[1]])
    elseif length(low) > 0
        return Gate(low=gate[low[1]])
    else
        return Gate(high=gate[high[1]])
    end
end

function parsegate(gate::AbstractArray)::Gate
    if length(gate) ≠ 2
        throw(ArgumentError("Prove an upper and lower limit for gate"))
    end
    return Gate(low=gate[1], high=gate[2])
end

## Parse Qkinz
function parseqkinz(path::AbstractString)::Vector{Poly}
    dir = dirname(path)
    pattern = Regex(basename(path))
    polys = Tuple{Int, Poly{Float64}}[]
    j(x) = joinpath(dir, x)
    for fname in readdir(dir)
        m = match(pattern, fname)
        if m ≢ nothing
            f = parse(Int, m[:f])
            push!(polys, (f, JSort.BetheBloch.Banana(j(fname)).exfromede))
        end
    end
    (f, poly) = polys |> sort |> x -> zip(x...)
    @assert length(f) == 8 "Missing Qkinz strips. Got $f"
    [p for p in poly]
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "jlqml"
version = v"0.1.3"

const is_yggdrasil = haskey(ENV, "BUILD_BUILDNUMBER")
git_repo = is_yggdrasil ? "https://github.com/barche/jlqml.git" : joinpath(ENV["HOME"], "src/julia/jlqml/.git")
unpack_target = is_yggdrasil ? "" : "jlqml"

# Collection of sources required to complete build
sources = [
    GitSource(git_repo, "7cfe4025dff8ecdbdffb1178943c53f073e9d212", unpack_target=unpack_target),
    ArchiveSource("https://github.com/JuliaBinaryWrappers/Julia_jll.jl/releases/download/Julia-v1.4.1+1/Julia.v1.4.1.x86_64-linux-gnu-libgfortran4-cxx11.tar.gz", "378b6a23ce4363eeb7afd5bd8092f902caa512f2f987dfc47fc51ae6bdff0e56"; unpack_target="julia-x86_64-linux-gnu"),
    ArchiveSource("https://github.com/JuliaBinaryWrappers/Julia_jll.jl/releases/download/Julia-v1.4.1+1/Julia.v1.4.1.x86_64-w64-mingw32-libgfortran4-cxx11.tar.gz", "621029838e895bf5f201d0858fdbd31f1bb7f458aa0bc0646b4b30185a7d8e7c"; unpack_target="julia-x86_64-w64-mingw32"),
    ArchiveSource("https://github.com/JuliaBinaryWrappers/Julia_jll.jl/releases/download/Julia-v1.4.1+1/Julia.v1.4.1.armv7l-linux-gnueabihf-libgfortran4-cxx11.tar.gz", "0d733c2e0147d6ffb731b638a8b1bd4225069c5735df22bc3a953dffce663d74"; unpack_target="julia-arm-linux-gnueabihf"),
    ArchiveSource("https://github.com/JuliaBinaryWrappers/Julia_jll.jl/releases/download/Julia-v1.4.1+1/Julia.v1.4.1.x86_64-apple-darwin14-libgfortran4-cxx11.tar.gz", "f6d94a3184b0241f20f78523de581949afa038f3e320fb9fd20a83019968adca"; unpack_target="julia-x86_64-apple-darwin14"),
]

# Bash recipe for building across all platforms
script = raw"""
if test -f "$prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake"; then
    sed -i 's/_qt5gui_find_extra_libs.*AGL.framework.*//' $prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake
fi

# Override compiler ID to silence the horrible "No features found" cmake error
if [[ $target == *"apple-darwin"* ]]; then
  macos_extra_flags="-DCMAKE_CXX_COMPILER_ID=AppleClang -DCMAKE_CXX_COMPILER_VERSION=10.0.0 -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
fi

Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target

mkdir build
cd build
cmake -DJulia_PREFIX=$Julia_PREFIX -DCMAKE_FIND_ROOT_PATH=$prefix -DJlCxx_DIR=$prefix/lib/cmake/JlCxx -DQt5Core_DIR=$prefix/lib/cmake/Qt5Core -DQt5Quick_DIR=$prefix/lib/cmake/Qt5Quick -DQt5Svg_DIR=$prefix/lib/cmake/Qt5Svg -DQt5Widgets_DIR=$prefix/lib/cmake/Qt5Widgets -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} $macos_extra_flags -DCMAKE_BUILD_TYPE=Release ../jlqml/
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
install_license $WORKSPACE/srcdir/jlqml*/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:armv7l; libc=:glibc, compiler_abi=CompilerABI(cxxstring_abi=:cxx11)),
    Linux(:x86_64; libc=:glibc, compiler_abi=CompilerABI(cxxstring_abi=:cxx11)),
    MacOS(:x86_64; compiler_abi=CompilerABI(cxxstring_abi=:cxx11)),
    Windows(:x86_64; compiler_abi=CompilerABI(cxxstring_abi=:cxx11)),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libjlqml", :libjlqml),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"),
    Dependency("Qt_jll"),
    BuildDependency("Libglvnd_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8")


text:
# import stuff
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using Statistics
@everywhere using GRASS
@everywhere using SharedArrays
@everywhere using EchelleCCFs
using CSV
using JLD2
using LsqFit
using FileIO
using DataFrames
using LaTeXStrings
using HypothesisTests

# define rms loop function
include(GRASS.moddir * "figures/fig_functions.jl")

# some global stuff
const N = 132
const Nloop = 1200

# get command line args and output directories
run, plot = parse_args(ARGS)
grassdir, plotdir, datadir = check_plot_dirs()

function main()
    # set observing parameters
    N_obs = range(2, 4, step=1)  # number of observations
    exp_time = 600.0             # ~2 p-mode periods
    snr = 5000.0                 # SNR per res element to get 50 cm/s precision
    new_res = 1.17e5

    # dead times
    short_dead = 30.0            # representative for NEID read-out time
    large_dead = 3600.0          # ~60 mins

    # set up parameters for spectrum
    lines = [5434.5]
    depths = [0.8]
    resolution = 7e5

    # allocate shared arrays for short off-target time
    avg_shortdead = SharedArray{Float64}(length(N_obs), Nloop)
    rms_shortdead = SharedArray{Float64}(length(N_obs), Nloop)

    # allocate shared arrays for large off-target time
    avg_largedead = SharedArray{Float64}(length(N_obs), Nloop)
    rms_largedead = SharedArray{Float64}(length(N_obs), Nloop)

    # create spec params instance
    spec1 = SpecParams(lines=lines, depths=depths, resolution=resolution)

    # loop over number of observations
    for i in eachindex(N_obs)
        # create observation plan instances
        obs_shortdead = GRASS.ObservationPlan(N_obs=N_obs[i],
                                              obs_per_night=N_obs[i],
                                              time_per_obs=exp_time,
                                              dead_time=short_dead)
        obs_largedead = GRASS.ObservationPlan(N_obs=N_obs[i],
                                              obs_per_night=N_obs[i],
                                              time_per_obs=exp_time,
                                              dead_time=large_dead)

        # compute vels repeatedly to get good stats
        @sync @distributed for j in 1:Nloop
            # first simulate short off-target time
            rvs1, sigs1 = GRASS.simulate_observations(obs_shortdead, spec1, N=N, snr=snr, new_res=new_res)
            avg_shortdead[i, j] = mean(rvs1)
            rms_shortdead[i, j] = calc_rms(rvs1)

            # then simulate large off-target time
            rvs2, sigs2 = GRASS.simulate_observations(obs_largedead, spec1, N=N, snr=snr, new_res=new_res)
            avg_largedead[i, j] = mean(rvs2)
            rms_largedead[i, j] = calc_rms(rvs2)
        end
    end

    # convert to plain arrays
    avg_shortdead = convert(Array{Float64}, avg_shortdead)
    rms_shortdead = convert(Array{Float64}, rms_shortdead)
    avg_largedead = convert(Array{Float64}, avg_largedead)
    rms_largedead = convert(Array{Float64}, rms_largedead)

    # save the output
    outfile = datadir * "observe_" * string(N) * "_loop_" * string(Nloop) * ".jld2"
    save(outfile,
         "avg_shortdead", avg_shortdead,
         "rms_shortdead", rms_shortdead,
         "avg_largedead", avg_largedead,
         "rms_largedead", rms_largedead,
         "exp_time", exp_time,
         "depths", depths,
         "N_obs", N_obs,
         "snr", snr)
    return nothing
end

# run the simulation
if run
    main()
end

# plotting code block
if plot
    # plotting imports
    import PyPlot; plt = PyPlot; mpl = plt.matplotlib; plt.ioff()
    using PyCall; animation = pyimport("matplotlib.animation")
    mpl.style.use(GRASS.moddir * "figures/fig.mplstyle")

    # read in the data
    file = datadir * "observe_" * string(N) * "_loop_" * string(Nloop) * ".jld2"
    d = load(file)
    avg_shortdead = d["avg_shortdead"]
    rms_shortdead = d["rms_shortdead"]
    avg_largedead = d["avg_largedead"]
    rms_largedead = d["rms_largedead"]
    exp_time = d["exp_time"]
    depths = d["depths"]
    N_obs = d["N_obs"]
    snr = d["snr"]

    # initialize plotting objects
    fig = plt.figure(figsize=(10.5, 6.75))
    ax1 = plt.subplot2grid((2,4),(0,0), colspan=2)
    ax2 = plt.subplot2grid((2,4),(0,2), colspan=2)
    ax3 = plt.subplot2grid((2,4),(1,1), colspan=2)
    axs = [ax1, ax2, ax3]

    # make arrays for limits and loop over plotting windows
    xlims = []
    ylims = []
    for i in eachindex(axs)
        # histogram for short dead time
        axs[i].hist(avg_shortdead[i,:], density=true, histtype="stepfilled", fc="k", ec="k", lw=2, alpha=0.5, label="Short wait")
        axs[i].hist(avg_shortdead[i,:], density=true, histtype="step", ec="k", lw=1.0, alpha=0.8)

        # histogram for long wait time
        axs[i].hist(avg_largedead[i,:], density=true, histtype="stepfilled", fc="tab:blue", ec="tab:blue", lw=1.5, alpha=0.5, label="Long wait")
        axs[i].hist(avg_largedead[i,:], density=true, histtype="step", ec="tab:blue", lw=1.5, alpha=0.8)

        # get axis limits
        append!(xlims, [axs[i].get_xlim()])
        append!(ylims, [axs[i].get_ylim()])

        # report two-sample KS test
        @show(ApproximateTwoSampleKSTest(avg_shortdead[i,:], avg_largedead[i,:]))
        println()
    end

    # set the axes limits
    bbox = Dict("ec"=>"black", "fc"=>"white")
    xy = (maximum(last.(xlims)) - 0.25, maximum(last.(ylims)) - 0.4)
    for i in eachindex(axs)
        # annotate axes and set labels, etc.
        axs[i].annotate(L"N_{\rm obs} =\ " * latexstring(N_obs[i]), xy=xy, bbox=bbox)
        axs[i].set_xlabel(L"{\rm Mean\ velocity\ (m s}^{-1})")
        axs[i].set_ylabel(L"{\rm Probability\ density}")
        axs[i].set_xlim(round(minimum(first.(xlims)), digits=1),
                        round(maximum(last.(xlims)), digits=1))
        axs[i].set_ylim(minimum(first.(ylims)), maximum(last.(ylims)))
        axs[i].xaxis.set_major_locator(mpl.ticker.MaxNLocator(5))
    end

    # set legend in one panel and write figure
    axs[1].legend(loc="upper left")
    fig.savefig(plotdir * "fig7.pdf")
    plt.clf(); plt.close()
    println(">>> Figure written to: " * plotdir * "fig7.pdf")
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
function reshop_reg_eqns(ctx, m::ReSHOPMathProgBaseModel, equil=false)
    if has_objective(m)
        obj_offset = 1
    else
        obj_offset = 0
    end

    reshop_declare_eqns(ctx, m, equil)

    if has_objective(m)
        reshop_add_obj_nl(ctx, m)
        reshop_add_obj_lin(ctx, m)
        reshop_add_obj_quad(ctx, m)
    end

    if m.ncon > 0 || length(m.quad_equs) > 0
#        write_nl_k_block(f, m)
        if length(m.nonquad_idx) != m.ncon
            error("The number of constraint and the equation index do not match")
        end
        reshop_add_cons_lin(ctx, m, obj_offset)
        reshop_add_cons_nl(ctx, m, obj_offset)
        reshop_add_cons_quad(ctx, m, obj_offset)
        reshop_add_contraint_sense(ctx, m, obj_offset)
    end

    # This is most likely broken, we need another way of getting this data if
    # it is important
#    return m.ncon+obj_offset+length(m.quad_equs)
end

function create_reshop_ctx(m::ReSHOPMathProgBaseModel)
    # Not all model have an objective function
    if has_objective(m)
        obj_offset = 1
    else
        obj_offset = 0
    end

    ctx = ctx_create(m.nvar, length(m.quad_equs)+m.ncon+obj_offset)

    reshop_declare_vars(ctx, m)
    reshop_add_var_guess(ctx, m)

    reshop_reg_eqns(ctx, m)

    return ctx
end

function has_objective(m::ReSHOPMathProgBaseModel)
    CONFIG[:debug] && println("DEBUG: has_objective $(m.lin_obj) $(m.obj) $(m.quad_obj)")
    res = !isempty(m.lin_obj) || (isa(m.obj, Expr)) || length(m.quad_obj) == 3
    if !res return res end
    res = any(v != 0. for v in values(m.lin_obj)) || isa(m.obj, Expr) && (m.obj != 0.) || length(m.quad_obj) == 3
    return res
end

# Nonlinear constraint trees
function reshop_add_cons_nl(ctx, m::ReSHOPMathProgBaseModel, offset)
    for idx in 1:m.ncon
        if isa(m.constrs[idx], Expr) && (m.constrs[idx] != 0.)
            eidx = m.nonquad_idx[idx] + offset - 1 + m.offset
            tree, node = reshop_get_treedata(ctx, eidx)
            CONFIG[:debug] && println("reshop_add_cons_nl: storing expression $(eidx)")
            reshop_add_nlexpr(ctx, tree, node, m, m.constrs[idx])
        end
    end
end

function reshop_quad(ctx, m::ReSHOPMathProgBaseModel, idx, equ, offset, isObj::Bool=false)
    ##########################################################################
    #
    ##########################################################################
    if isObj
        eidx = m.offset
    else
        eidx = m.quad_idx[idx] - 1 + m.offset + offset
    end

    CONFIG[:debug] && println("reshop_quad: $equ")
    lidx, lval, rowidx, colidx, qval = equ[1:5]
    CONFIG[:debug] && println("reshop_quad: $lidx $lval $rowidx $colidx $qval")

    ##########################################################################
    # lin_dict collects all the linear terms
    #
    # we don't construct the dict directly, since it may have duplicate values
    ##########################################################################

    lin_dict = Dict(zip(lidx, Iterators.repeated(0.)))
    for i in 1:length(lidx)
        lin_dict[lidx[i]] += lval[i]
    end

    ###########################################################################
    # MPB is a bit insane here: via setquadobj! it provides a triangular matrix
    # representing .5 x^TQx ... But this makes no sense since we assume that
    # the multiplication is commutative here. Therefore the values of the off-
    # diagonal element have to be multiplied by 2
    ###########################################################################

    quad_dict = Dict(zip(zip(rowidx, colidx), Iterators.repeated(0.)))

    if isObj
        for i in 1:length(rowidx)
            if rowidx[i] != colidx[i]
                quad_dict[(rowidx[i], colidx[i])] += 2*qval[i]
            else
                quad_dict[(rowidx[i], colidx[i])] += qval[i]
            end
        end
    else
        for i in 1:length(rowidx)
            quad_dict[(rowidx[i], colidx[i])] += qval[i]
        end
    end

    qidxC = keys(quad_dict)
    rowidxU = [elt[1]-1 for elt in qidxC]
    colidxU = [elt[2]-1 for elt in qidxC]
    qvalU = collect(values(quad_dict))

    ##########################################################################
    # Add the linear part <c,x> for the purely linear variables
    ##########################################################################

    lidxS = Set(lidx)
    lin_vars = setdiff(lidxS, union(Set(rowidx), Set(colidx)))
    qidxS = setdiff(lidxS, lin_vars)
    for lindx in lin_vars
        # add as linear variable
        v = lin_dict[lindx]
        if abs(v) > 0.
            CONFIG[:debug] && println("DEBUG: in reshop_quad, lin var $lindx with value $v")
            ctx_add_lin_var(ctx, eidx, m.v_index_map[lindx], v)
        end
    end

    ##########################################################################
    # Add the term <c,x> for the quadratic variables
    ##########################################################################

    if length(lin_dict) > 0 && length(qidxS) > 0
        # TODO(xhub) filter for qvals[i] = 0.
        qvals = collect(lin_dict[qidx] for qidx in qidxS)
        CONFIG[:debug] && println("DEBUG: in reshop_quad, qidxS = $qidxS; qval= $qvals")
        abs_var = reshop_avar(length(qidxS), collect(qidxS) .- 1)
        CONFIG[:debug] && println("DEBUG: in reshop_quad, quad var $(collect(qidxS)) with value $qvals")
        rhp_equ_add_lin_tree(ctx, eidx, qvals, abs_var, 1.)
        reshop_avar_free(abs_var)
    end

    ##########################################################################
    # Add the quadratic terms <x, Mx>
    ##########################################################################

    if length(qvalU) > 0
        CONFIG[:debug] && println("DEBUG: in reshop_quad, quad term $rowidxU $colidxU $qvalU")
        @assert length(qvalU) == length(rowidxU)
        @assert length(qvalU) == length(colidxU)
        mat = reshop_mat_coo(rowidxU, colidxU, qvalU)
        midxS = union(BitSet(rowidx), BitSet(colidx))
        avar = reshop_avar(length(midxS), collect(i-1 for i in midxS))
        if isObj
            c = 1.
        else
            c = 2.
        end
        # the coeff is 2 since this function adds .5 x^TMx
        rhp_equ_add_quadratic(ctx, eidx, mat, avar, c)
        reshop_avar_free(avar)
        reshop_mat_free(mat)
    end
end

function reshop_add_obj_quad(ctx, m::ReSHOPMathProgBaseModel)
    CONFIG[:debug] && println("DEBUG: quad_obj = $(m.quad_obj)")
    if length(m.quad_obj) == 3
        rowidx = m.quad_obj[1]
        colidx = m.quad_obj[2]
        quadval = m.quad_obj[3]
        reshop_quad(ctx, m, 0, (Int[], Float64[], rowidx, colidx, quadval), 0, true)
    elseif length(m.quad_obj) == 0
        #doing noting here
    else
        error("reshop_add_obj_quad :: invalid quad_obj object $(m.quad_obj)")
    end
end

function reshop_add_cons_quad(ctx, m::ReSHOPMathProgBaseModel, offset)
    for (idx, equ) in enumerate(m.quad_equs)
        reshop_quad(ctx, m, idx, equ, offset)
    end
end

# Nonlinear objective tree
function reshop_add_obj_nl(ctx, m::ReSHOPMathProgBaseModel)
    # Get tree, node, ...
    if m.obj != 0.
        tree, node = reshop_get_treedata(ctx, m.offset)
        CONFIG[:debug] && println("reshop_add_obj_nl: storing expression $(m.offset)")
        reshop_add_nlexpr(ctx, tree, node, m, m.obj)
    end
end

# Initial primal guesses
function reshop_add_var_guess(ctx, m::ReSHOPMathProgBaseModel)
    for idx in 0:(m.nvar - 1)
        i = m.v_index_map_rev[idx]
        CONFIG[:debug] && println("Setting level value variable $i to $(m.x_0[i])")
        ctx_setvarval(ctx, idx, m.x_0[i])
    end
end

# Constraint bounds
function reshop_add_contraint_sense(ctx, m::ReSHOPMathProgBaseModel, offset)
    for idx in 1:m.ncon
        lower = m.g_l[idx]
        upper = m.g_u[idx]
        rel = m.r_codes[idx]
        # TODO(xhub) document
        if rel == -1
            error("Doubly constraint equation is not yet supported")
        elseif rel == 4         # ==
            value = lower
        elseif rel == 1         # >=
            value = lower
        elseif rel == 2         # <=
            value = upper
        elseif rel == 0         # ???
            error("unsupported rel = $(rel) for equation index $idx")
            value == upper
        else
            error("unsupported rel = $(rel) for equation index $idx")
        end
        eidx = m.nonquad_idx[idx] + offset - 1 + m.offset
        CONFIG[:debug] && println("Setting sense and rhs for equation $eidx: $rel $value")
        reshop_set_rhs(ctx, eidx, value)
        reshop_set_equtype(ctx, eidx, rel)
    end

    for (idx, equ) in enumerate(m.quad_equs)
        rel, value = equ[end-1:end]
        eidx = m.quad_idx[idx] + offset - 1 + m.offset
        reshop_set_rhs(ctx, eidx, value)
        reshop_set_equtype(ctx, eidx, relation_to_reshop[rel])
    end
end

function reshop_declare_var(ctx, vtype, lower, upper)
    if lower == -Inf
        if upper == Inf
            reshop_add_free_var(ctx, 1)
        elseif upper == 0
            reshop_add_neg_var(ctx, 1)
        else
            reshop_add_box_var(ctx, lower, upper)
        end
    else
        if lower == upper
            reshop_add_box_var(ctx, lower, upper)
        elseif upper == Inf
            if lower == 0
                reshop_add_pos_var(ctx, 1)
            else
                reshop_add_box_var(ctx, lower, upper)
            end
        else
            reshop_add_box_var(ctx, lower, upper)
        end
    end
    idx = hack_last_vidx(ctx)
    CONFIG[:debug] && println("DEBUG: var $(idx) has type $(vtype)")
    # We have no consistency check
    # TODO SOS, semicont
    if vtype == :Bin
        reshop_set_vartype(ctx, idx, 1)
    elseif vtype == :Int
        reshop_set_vartype(ctx, idx, 2)
    end

end

# Variable bounds
function reshop_declare_vars(ctx, m::ReSHOPMathProgBaseModel)
    for idx in 0:(m.nvar - 1)
        i = m.v_index_map_rev[idx]
        lower = m.x_l[i]
        upper = m.x_u[i]
        reshop_declare_var(ctx, m.vartypes[i], lower, upper)
    end

    # Set variable names
    if !isnothing(m.d)
        ctx_setvarnames(ctx, m.d.m.colNames)
    end
end

# Jacobian counts
# TODO(xhub) resuse this to prealloc the model_repr?
#function write_nl_k_block(f, m::ReSHOPMathProgBaseModel)
#    println(f, "k$(m.nvar - 1)")
#    total = 0
#    for index = 0:(m.nvar - 2)
#        i = m.v_index_map_rev[index]
#        total += m.j_counts[i]
#        println(f, total)
#    end
#end

function reshop_declare_eqns(ctx, m::ReSHOPMathProgBaseModel, equil=false)
    if has_objective(m)
        cor = 0
    else
        cor = -1
    end

    for idx in 0:m.ncon+cor+length(m.quad_equs)
        reshop_decl_eqn(ctx)
    end
    if equil || !has_objective(m)
        rhp_set_objeqn(ctx, -1)
    else
        rhp_set_objeqn(ctx, 0)
    end
#    return m.ncon+cor+length(m.quad_equs)
end

# Linear constraint expressions
function reshop_add_cons_lin(ctx, m::ReSHOPMathProgBaseModel, offset)
    for idx in 1:m.ncon
        num_vars = length(m.lin_constrs[idx])
        if num_vars > 0
#TODO(Xhub)            reshop_reserve_lequ(ctx, idx, num_vars)
            eidx = m.nonquad_idx[idx] + offset - 1 + m.offset
            for (k,v) in m.lin_constrs[idx]
                if abs(v) > 0.
                    CONFIG[:debug] && println("DEBUG: add var $k with value $v in eqn $(eidx)")
                    ctx_add_lin_var(ctx, eidx, m.v_index_map[k], v)
                end
            end
        end
    end
end

# Linear objective expression
function reshop_add_obj_lin(ctx, m::ReSHOPMathProgBaseModel)
    for (k,v) in m.lin_obj
        if abs(v) > 0.
            ctx_add_lin_var(ctx, m.offset, m.v_index_map[k], v)
        end
    end
end

function reshop_add_quad_lin(ctx, m::ReSHOPMathProgBaseModel)
end

function reshop_set_modeltype(m::ReSHOPMathProgBaseModel)
    discrete = any((m.vartypes .== :Int) + (m.vartypes .== :Bin) .> 0)
    if discrete
        if m.model_type == qcp
            m.model_type = miqcp
        elseif m.model_type == nlp
             m.model_type = minlp
        end
    end
    reshop_set_objsense(m.reshop_ctx, m.sense)
    reshop_set_modeltype(m.reshop_ctx, m.model_type)
end


text:
"""
    const PointType = UInt8

Stores certain information about a grid point via bit-flags.

Right now, there are:

    const update_bit      = 0x01
    const undepleted_bit  = 0x02
    const pn_junction_bit = 0x04
    const bulk_bit     = 0x08

## Examples

How to get information out of a `PointType` variable `point_type`:
1. `point_type & update_bit == 0` -> do not update this point (for fixed points)
2. `point_type & update_bit >  0` -> do update this point
3. `point_type & undepleted_bit > 0` -> this point is undepleted
4. `point_type & pn_junction_bit > 0` -> this point belongs to the solid state detector, meaning that it is in the volume of the n-type or p-type material.
5. `point_type & bulk_bit > 0` -> this point is only surrounded by points marked as `pn_junction_bit`
"""
const PointType       = UInt8

# Point types for electric potential calculation
const update_bit      = 0x01 # parse(UInt8, "00000001", base=2) # 1 -> do update; 0 -> do not update
const undepleted_bit  = 0x02 # parse(UInt8, "00000010", base=2) # 0 -> depleted point; 1 -> undepleted point
const pn_junction_bit = 0x04 # parse(UInt8, "00000100", base=2) # 0 -> point is not part of pn-junction; 1 -> point is part of the pn-junction
const bulk_bit     = 0x08 # parse(UInt8, "00001000", base=2) # 0 -> point is surrounded by points that do not belong to the pn-junction; 1 -> point is only surrounded by points in the pn-junction

is_pn_junction_point_type(p::PointType) = p & pn_junction_bit > 0
is_undepleted_point_type(p::PointType) = p & undepleted_bit > 0
is_fixed_point_type(p::PointType) = p & update_bit == 0
"""
    struct PointTypes{T, N, S, AT} <: AbstractArray{T, N}

Information about the grid points used to calculate the [`ElectricPotential`](@ref)
stored via bit-flags. Data is stored as [`PointType`](@ref) which is an `UInt8`.

## Parametric types
* `T`: Element type of `grid.axes`.
* `N`: Dimension of the `grid` and `data` array.
* `S`: Coordinate system (`Cartesian` or `Cylindrical`).
* `AT`: Axes type.

## Fields
* `data::Array{PointType, N}`: Array containing the point type values at the discrete points of the `grid`.
* `grid::Grid{T, N, S, AT}`: [`Grid`](@ref) defining the discrete points for which the point types are determined.

See also [`PointType`](@ref).
"""
struct PointTypes{T, N, S, AT} <: AbstractArray{T, N}
    data::Array{PointType, N}
    grid::Grid{T, N, S, AT}
end

@inline size(point_types::PointTypes{T, N, S}) where {T, N, S} = size(point_types.data)
@inline length(point_types::PointTypes{T, N, S}) where {T, N, S} = length(point_types.data)
@inline getindex(point_types::PointTypes{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(point_types.data, I...)
@inline getindex(point_types::PointTypes{T, N, S}, i::Int) where {T, N, S} = getindex(point_types.data, i)
@inline getindex(point_types::PointTypes{T, N, S}, s::Symbol) where {T, N, S} = getindex(point_types.grid, s)
@inline getindex(point_types::PointTypes{T, N, S}, pt::AbstractCoordinatePoint{T}) where {T, N, S} = point_types.data[find_closest_gridpoint(pt, point_types.grid)...]

function in(pt::AbstractCoordinatePoint{T}, point_types::PointTypes{T, 3, S})::Bool where {T <: SSDFloat, S}
    cpt = _convert_point(pt, S)
    point_type::PointType = point_types[cpt]
    return point_type & bulk_bit > 0
end

function mark_bulk_bits!(point_types::Array{PointType, 3})
    i1max, i2max, i3max = size(point_types)
    for i1 in Base.OneTo(i1max), i2 in Base.OneTo(i2max), i3 in Base.OneTo(i3max)
        (point_types[i1,i2,i3] & pn_junction_bit == 0) && continue
        point_types[i1,i2,i3] += bulk_bit * begin
            in_bulk = true
            for j1 in max(i1-1,1):min(i1+1,i1max), j2 in max(i2-1,1):min(i2+1,i2max), j3 in max(i3-1,1):min(i3+1,i3max)
                in_bulk &= point_types[j1,j2,j3] & pn_junction_bit > 0
                !in_bulk && break
            end
            in_bulk
        end
    end
    point_types
end

"""
    is_depleted(point_types::PointTypes)::Bool

Returns `true` if all [`PointType`](@ref) values of
the [`PointTypes`](@ref) of a [`Simulation`](@ref) are marked as depleted
and `false` if any point in the [`PointTypes`](@ref) is marked as undepleted.

It can be used to determine whether the [`SolidStateDetector`](@ref) is
depleted at the provided bias voltage.

## Arguments
* `point_types::PointTypes`: [`PointTypes`](@ref) of a [`Simulation`](@ref).

## Examples
```julia
is_depleted(sim.point_types)
```
"""
is_depleted(point_types::PointTypes)::Bool =
    !any(b -> undepleted_bit & b > 0, point_types.data)

"""
    get_active_volume(point_types::PointTypes{T}) where {T}

Returns an approximation of the active volume of the detector by summing up the cell volumes of
all cells marked as depleted.

## Arguments
* `point_types::PointTypes{T}`: Point types of a [`Simulation`](@ref).

## Examples
    get_active_volume(sim.point_types)

!!! note
    Only `φ`-symmetries are taken into account.
"""
function get_active_volume(point_types::PointTypes{T, 3, Cylindrical}) where {T}
    active_volume::T = 0

    only_2d::Bool = length(point_types.grid.axes[2]) == 1
    cyclic::T = point_types.grid.axes[2].interval.right

    r_ext::Vector{T} = get_extended_ticks(point_types.grid.axes[1])
    φ_ext::Vector{T} = get_extended_ticks(point_types.grid.axes[2])
    z_ext::Vector{T} = get_extended_ticks(point_types.grid.axes[3])

    mpr::Vector{T} = midpoints(r_ext)
    mpφ::Vector{T} = midpoints(φ_ext)
    mpz::Vector{T} = midpoints(z_ext)
    Δmpφ::Vector{T} = diff(mpφ)
    Δmpz::Vector{T} = diff(mpz)
    Δmpr_squared::Vector{T} = T(0.5) .* ((mpr[2:end].^2) .- (mpr[1:end-1].^2))
    if r_ext[2] == 0
        Δmpr_squared[1] = T(0.5) * (mpr[2]^2)
    end

    isclosed::Bool = typeof(point_types.grid.axes[2].interval).parameters[2] == :closed
    for iz in eachindex(point_types.grid.axes[3])
        if !isclosed || only_2d
            for iφ in eachindex(point_types.grid.axes[2])
                for ir in eachindex(point_types.grid.axes[1])
                    point_type::PointType = point_types[ir, iφ, iz]
                    if (point_type & pn_junction_bit > 0) && (point_type & undepleted_bit == 0) && (point_type & update_bit > 0)
                        dV::T = Δmpz[iz] * Δmpφ[iφ] * Δmpr_squared[ir]
                        active_volume += dV
                    end
                end
            end
        elseif isclosed && !only_2d
            for iφ in eachindex(point_types.grid.axes[2])
                for ir in eachindex(point_types.grid.axes[1])
                    point_type::PointType = point_types[ir, iφ, iz]
                    if (point_type & pn_junction_bit > 0) && (point_type & undepleted_bit == 0) && (point_type & update_bit > 0)
                        dV = Δmpz[iz] * Δmpφ[iφ] * Δmpr_squared[ir]
                        active_volume += if iφ == length(point_types.φ) || iφ == 1
                            dV / 2
                        else
                            dV
                        end
                    end
                end
            end
        end
    end
    if cyclic > 0
        active_volume *= 2π / cyclic
    end

    f::T = 10^6
    return active_volume * f * Unitful.cm * Unitful.cm * Unitful.cm
end

function get_active_volume(point_types::PointTypes{T, 3, Cartesian}) where {T}
    active_volume::T = 0

    x_ext::Vector{T} = get_extended_ticks(point_types.grid.axes[1])
    y_ext::Vector{T} = get_extended_ticks(point_types.grid.axes[2])
    z_ext::Vector{T} = get_extended_ticks(point_types.grid.axes[3])

    mpx::Vector{T} = midpoints(x_ext)
    mpy::Vector{T} = midpoints(y_ext)
    mpz::Vector{T} = midpoints(z_ext)
    Δmpx::Vector{T} = diff(mpx)
    Δmpy::Vector{T} = diff(mpy)
    Δmpz::Vector{T} = diff(mpz)

    for iz in eachindex(point_types.grid.axes[3])
        for iy in eachindex(point_types.grid.axes[2])
            for ix in eachindex(point_types.grid.axes[1])
                point_type::PointType = point_types[ix, iy, iz]
                if (point_type & pn_junction_bit > 0) && (point_type & undepleted_bit == 0) && (point_type & update_bit > 0)
                    dV = Δmpx[ix] * Δmpy[iy] * Δmpz[iz]
                    active_volume += dV
                end
            end
        end
    end

    f::T = 10^6
    return active_volume * f * Unitful.cm * Unitful.cm * Unitful.cm
end

function PointTypes(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(grid.axes[1].ticks[1])
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    PointTypes{T, N, S, typeof(grid.axes)}( nt.values, grid )
end

Base.convert(T::Type{PointTypes}, x::NamedTuple) = T(x)

function NamedTuple(point_types::PointTypes{T, 3}) where {T}
    return (
        grid = NamedTuple(point_types.grid),
        values = point_types.data,
    )
end

Base.convert(T::Type{NamedTuple}, x::PointTypes) = T(x)


text:
module CSEM1DEr
include("DigFilters.jl")
using DataInterpolations, LinearAlgebra

abstract type RadialEr end

mutable struct RadialErLagged <: RadialEr
    thickness :: Array{Float64, 1}
    pz        :: Array{Complex{Float64}, 1}
    epsc      :: Array{ComplexF64, 1}
    zintfc    :: Array{Float64, 1}
    rTE       :: Array{ComplexF64, 1}
    rTM       :: Array{ComplexF64, 1}
    zRx       :: Array{Float64, 1}
    zTx       :: Array{Float64, 1}
    rRx       :: Array{Float64, 1}
    freqs     :: Array{Float64, 1}
    Er        :: Array{ComplexF64, 2}
    ErhJ0     :: Array{ComplexF64, 2}
    ErhJ1     :: Array{ComplexF64, 2}
    ErvJ1     :: Array{ComplexF64, 2}
    ErBase    :: Array{ComplexF64, 2}
    lambdaR   :: Array{Float64, 1}
    rR        :: Array{Float64, 1}
    cst       :: Float64
    sit       :: Float64
    csb       :: Float64
    sib       :: Float64
end

const mu      = 4*pi*1e-7
const eps0    = 8.854e-12

function RadialErLagged(;
      nmax      = 200,
      zTx       = [975.],
      rRx       = collect(LinRange(500, 5000, 20)),
      freqs     = [0.1, 0.3, 0.7],
      zRx       = [1000.],
      TxDip     = 0.0,
      RxAzim    = 0.0
  )
    @assert all(freqs .> 0.)
    thickness = zeros(nmax)
    zintfc  = zeros(nmax)
    pz  = zeros(Complex{Float64}, nmax)
    epsc = similar(pz)
    rTE = zeros(length(pz)-1)
    rTM = similar(rTE)
    Er  = zeros(ComplexF64, length(rRx),length(freqs)) # space domain fields
    #lagged convolution, good for all Rx-Tx depths same
    rMin, rMax          = extrema(rRx)
    lamMin              = Filter_base[1]/rMax
    lamMax              = Filter_base[end]/rMin
    # the filters are spaced by a multiplicative factor
    filterSpacing       = Filter_base[2] / Filter_base[1]
    nLambda             = ceil(log(lamMax/lamMin)/log(filterSpacing))+1
    lambdaR             = collect(exp.(log(lamMin):log(filterSpacing):log(lamMin)+log(filterSpacing)*(nLambda-1)))
    nExtra              = nLambda - length(Filter_base)
    # ranges corresponding to lambdaR plus nExtra 1 lags in convolution:
    rR                  = rMax*exp.(-(0:nExtra)*log(filterSpacing))
    # define kr domain fields here for lagged convolution
    ErhJ0, ErhJ1, ErvJ1 = map(x->zeros(ComplexF64, length(lambdaR), length(freqs)), 1:3)
    # needed for spline computations in space domain
    ErBase              = zeros(ComplexF64, length(rR), length(freqs))
    # Dip and azimuth
    # Rx Azimuths and Tx Dips
    cst = cosd(TxDip)  # needed for all HED terms
    sit = sind(TxDip)  # needed for all VED terms
    csb = cosd(RxAzim) # needs to be multiplied into all HED/VED terms inside the loop
    sib = sind(RxAzim) # can be multiplied into all non-VED terms outside the loop
    RadialErLagged(thickness, pz, epsc, zintfc, rTE, rTM, zRx, zTx, rRx, freqs, Er,
            ErhJ0, ErhJ1, ErvJ1, ErBase, lambdaR, rR,
            cst, sit, csb, sib)
end

mutable struct RadialErDHT <: RadialEr
    thickness :: Array{Float64, 1}
    pz        :: Array{Complex{Float64}, 1}
    epsc      :: Array{ComplexF64, 1}
    zintfc    :: Array{Float64, 1}
    rTE       :: Array{ComplexF64, 1}
    rTM       :: Array{ComplexF64, 1}
    zRx       :: Array{Float64, 1}
    zTx       :: Array{Float64, 1}
    rRx       :: Array{Float64, 1}
    freqs     :: Array{Float64, 1}
    Er        :: Array{ComplexF64, 2}
    ErhJ0     :: Array{ComplexF64, 2}
    ErhJ1     :: Array{ComplexF64, 2}
    ErvJ1     :: Array{ComplexF64, 2}
    lambda    :: Array{Float64, 1}
    cst       :: Array{Float64, 1}
    sit       :: Array{Float64, 1}
    csb       :: Array{Float64, 1}
    sib       :: Array{Float64, 1}
end

function RadialErDHT(;
      nmax      = 200,
      rRx       = collect(LinRange(500, 5000, 20)),
      freqs     = [0.1, 0.3, 0.7],
      zTx       = 975*ones(length(rRx)),
      zRx       = 1000*ones(length(rRx)),
      TxDip     = zeros(length(zTx)),
      RxAzim    = zeros(length(rRx))
  )
    @assert all(freqs .> 0.)
    @assert length(zTx) == length(zRx)
    thickness = zeros(nmax)
    zintfc  = zeros(nmax)
    pz  = zeros(Complex{Float64}, nmax)
    epsc = similar(pz)
    rTE = zeros(length(pz)-1)
    rTM = similar(rTE)
    Er  = zeros(ComplexF64, length(rRx),length(freqs)) # space domain fields
    lambda = zeros(length(Filter_base))
    # define kr domain fields here for lagged convolution
    ErhJ0, ErhJ1, ErvJ1 = map(x->zeros(ComplexF64, length(lambda), length(freqs)), 1:3)
    # Dip and azimuth
    # Rx Azimuths and Tx Dips
    cst = cosd.(TxDip)  # needed for all HED terms
    sit = sind.(TxDip)  # needed for all VED terms
    csb = cosd.(RxAzim) # needs to be multiplied into all HED/VED terms inside the loop
    sib = sind.(RxAzim) # can be multiplied into all non-VED terms outside the loop
    RadialErDHT(thickness, pz, epsc, zintfc, rTE, rTM, zRx, zTx, rRx, freqs, Er,
            ErhJ0, ErhJ1, ErvJ1, lambda, cst, sit, csb, sib)
end

function stacks!(F::RadialEr, iTxLayer::Int, nlayers::Int, omega::Float64)

    rTE              = F.rTE
    rTM              = F.rTM
    pz               = view(F.pz, 1:nlayers)
    #The last and first layer thicknesses are infinite
    d                = view(F.thickness, 1:nlayers)
    d[1]             = 1e60
    d[nlayers]       = 1e60

    # Capital R is for a stack
    # Starting from the bottom up, for Rs_down
    Rlowerstack_TE, Rlowerstack_TM = Complex(0., 0.), Complex(0., 0.)
    for k = (nlayers-1):-1:iTxLayer
      Rlowerstack_TE = lowerstack(Rlowerstack_TE, pz, rTE, d, k, omega)
      Rlowerstack_TM = lowerstack(Rlowerstack_TM, pz, rTM, d, k, omega)
    end

    Rupperstack_TE, Rupperstack_TM = Complex(0., 0.), Complex(0., 0.)
    # Starting from the top down for Rs_up
    for k   = 2:iTxLayer
        Rupperstack_TE = upperstack(Rupperstack_TE, pz, rTE, d, k, omega)
        Rupperstack_TM = upperstack(Rupperstack_TM, pz, rTM, d, k, omega)
    end
return Rupperstack_TE, Rlowerstack_TE, Rupperstack_TM, Rlowerstack_TM
end

function lowerstack(Rlowerstack::ComplexF64, pz::SubArray{ComplexF64, 1},
                    r::Array{ComplexF64, 1}, d::SubArray{Float64, 1}, k::Int, omega::Float64)
    a  = Rlowerstack*exp(2im*omega*pz[k+1]*d[k+1])
    Rs_d = (r[k] + a) / (1. + r[k]*a)
end

function upperstack(Rupperstack::ComplexF64, pz::SubArray{ComplexF64, 1},
                    r::Array{ComplexF64, 1}, d::SubArray{Float64, 1}, k::Int, omega::Float64)
    a = Rupperstack*exp(2im*omega*pz[k-1]*d[k-1])
    Rs_u = (-r[k-1] + a) / (1. - r[k-1]*a)
end

function getCurlyR(Rs_u::ComplexF64, Rs_d::ComplexF64, pz::ComplexF64,
                  zR::Float64, z::SubArray{Float64, 1}, iTxLayer::Int, omega::Float64)
    d=z[iTxLayer+1]-z[iTxLayer]
    denom = (1. - Rs_u*Rs_d*exp(2im*omega*pz*(d)))
    if (zR>=0)
      a = exp(-im*omega*pz*2*z[iTxLayer])
      b = exp(im*omega*pz*zR)
      c = exp(im*omega*pz*(2*z[iTxLayer+1] - zR))
      finRA = (1. + a*Rs_u) * (b + c*Rs_d) / denom
      finRB = (1. - a*Rs_u) * (b + c*Rs_d) / denom
    else
      a = exp(im*omega*pz*2*z[iTxLayer+1])
      b = exp(-im*omega*pz*zR)
      c = exp(-im*omega*pz*(2*z[iTxLayer] - zR))
      finRA = (1. + a*Rs_d) * (b + c*Rs_u) /  denom
      finRB = (-1. + a*Rs_d) * (b + c*Rs_u) / denom
    end

    return finRA, finRB
end

getepsc(rho, omega)      = eps0 + 1im/(rho*omega)
getpz(epsc, krho, omega) = unsafesqrt(mu*epsc - (krho/omega)^2)
checkz(z, zcheck)        = round(Int, z<zcheck)
ztxorignify(z, zTx)      = z - zTx
makesane(pz::Complex)    = imag(pz)  < 0.0 ? ( pz*=-1.) : pz
function unsafesqrt(z::ComplexF64)
    # x = real(z)
    # d = sqrt(x^2+imag(z)^2)
    # ComplexF64(sqrt(0.5(d+x)), sqrt(0.5(d-x)))
    x, y = reim(z)
    ρ = x*x + y*y
    ρ=abs(x)+sqrt(ρ)
    k = -1
    ρ += ρ
    ρ = ldexp(sqrt(ρ),k) #sqrt((abs(z)+abs(x))/2) without over/underflow
    ξ = ρ
    η = y
    η=(η/ρ)/2
    if x<0
        ξ = abs(η)
        η = copysign(ρ,y)
    end
    Complex(ξ,η)
end
function getCSEM1DKernelsEr!(F::RadialEr, krho::Float64, f::Float64, zz::Array{Float64, 1}, rho::Array{Float64, 1},
                            rxno::Int, txno::Int)
    nlayers = length(rho)
    z       = view(F.zintfc, 1:nlayers)
    zRx     = ztxorignify(F.zRx[rxno], F.zTx[txno])
    omega   = 2. *pi*f
    epsc    = F.epsc
    pz      = F.pz
    # reflection coefficients (downward) for an intfc: pz vertical slowness
    rTE, rTM = F.rTE, F.rTM

    l = 1
    z[l]       = ztxorignify(zz[l], F.zTx[txno])
    epsc[l]    = getepsc(rho[l], omega)
    pz[l]      = getpz(epsc[l], krho, omega)
    iTxLayer   =  checkz(z[l], 0)
    iRxLayer   =  checkz(z[l], zRx)
    @inbounds @fastmath for intfc in 1:nlayers-1
        l = intfc+1
        z[l]       = ztxorignify(zz[l], F.zTx[txno])
        epsc[l]    = getepsc(rho[l], omega)
        pz[l]      = getpz(epsc[l], krho, omega)
        iTxLayer  += checkz(z[l], 0)
        iRxLayer  += checkz(z[l], zRx)
        rTE[intfc] = (pz[intfc] - pz[intfc+1])/(pz[intfc] + pz[intfc+1])
        rTM[intfc] = (epsc[intfc]*pz[intfc+1] - epsc[intfc+1]*pz[intfc]) /
                   (epsc[intfc]*pz[intfc+1] + epsc[intfc+1]*pz[intfc])
        F.thickness[intfc] = z[intfc+1] - z[intfc]
    end
    # Find the layer containing the transmitter:
    @assert iTxLayer == iRxLayer "Error, transmitter and receiver need to be located in the same layer! Stopping"
    if iTxLayer ==  0
        iTxLayer = length(z) # in the last layer
        iRxLayer = length(z) # in the last layer
    end
    # TE and TM modes
    Rs_uTE, Rs_dTE, Rs_uTM, Rs_dTM  = stacks!(F, iTxLayer, nlayers, omega)
    curlyRA,         = getCurlyR(Rs_uTE, Rs_dTE, pz[iTxLayer], zRx, z, iTxLayer, omega)
    gA_TE            = mu/pz[iTxLayer]*curlyRA

    curlyRA, curlyRB = getCurlyR(Rs_uTM, Rs_dTM, pz[iTxLayer], zRx, z, iTxLayer, omega)
    gA_TM            = pz[iTxLayer]/epsc[iTxLayer]*curlyRA
    gB_TM            = curlyRB

    # Kernels according to Loseth, without the bessel functions multiplied
    # Er from HED and VED
    ErhJ0              = -krho*gA_TM/4/pi
    ErhJ1              = -(gA_TE - gA_TM)/4/pi
    ErvJ1              = krho^2*gB_TM*1im/(4*pi*omega*epsc[iTxLayer])

    return ErhJ0, ErhJ1, ErvJ1
end

function getfield!(F::RadialErLagged, z::Array{Float64, 1}, ρ::Array{Float64, 1})
    ErhJ0, ErhJ1, ErvJ1, ErBase = F.ErhJ0, F.ErhJ1, F.ErvJ1, F.ErBase
    for (ifreq, freq) in enumerate(F.freqs), (ikr, kr) in enumerate(F.lambdaR)
        ErhJ0[ikr,ifreq],
        ErhJ1[ikr,ifreq],
        ErvJ1[ikr,ifreq] = getCSEM1DKernelsEr!(F, kr, freq, z, ρ, 1, 1)
    end # first loop over freq for evaluations at digital filter specified kr

    for (ifreq, freq) in enumerate(F.freqs)
        for (ir, r) in enumerate(F.rR)
        # Extract results for all rR ranges:
            iShift = (1:length(Filter_base)) .+ ir .-1
            ErBase[ir,ifreq] = F.cst*F.csb*dot(ErhJ0[iShift,ifreq], Filter_J0)/r +
                               F.cst*F.csb*dot(ErhJ1[iShift,ifreq], Filter_J1)/r^2 +
                               F.sit*      dot(ErvJ1[iShift,ifreq], Filter_J1)/r
        end

        # fit spline at computed receiver locations
        # Dierckx needs x to be increasing so end:-1:1
        # Er
        splReal = CubicSpline(real(ErBase[end:-1:1,ifreq]), log10.(F.rR[end:-1:1]))
        splImag = CubicSpline(imag(ErBase[end:-1:1,ifreq]), log10.(F.rR[end:-1:1]))
        # evaluate spline at required rRx
        F.Er[:,ifreq] = splReal.(log10.(F.rRx)) - 1im*splImag.(log10.(F.rRx))
    end # loop over frequencies
    # done with lagged convolution
end

function getfield!(F::RadialErDHT, z::Array{Float64, 1}, ρ::Array{Float64, 1})
    ErhJ0, ErhJ1, ErvJ1 = F.ErhJ0, F.ErhJ1, F.ErvJ1
    for (ifreq, freq) in enumerate(F.freqs)
        for (ir, r) in enumerate(F.rRx)
            copy!(F.lambda, Filter_base/r)
            for (ikr, kr) in enumerate(F.lambda)
                ErhJ0[ikr,ifreq],
                ErhJ1[ikr,ifreq],
                ErvJ1[ikr,ifreq] = getCSEM1DKernelsEr!(F, kr, freq, z, ρ, 1, ir)
            end # kr loop
            F.Er[ir,ifreq] = F.cst[ir]*F.csb[ir]*dot(ErhJ0[:,ifreq], Filter_J0)/r +
                             F.cst[ir]*F.csb[ir]*dot(ErhJ1[:,ifreq], Filter_J1)/r^2 +
                             F.sit[ir]*          dot(ErvJ1[:,ifreq], Filter_J1)/r
        end # rRx loop
    end # freq loop
end

end # module

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module ZygoteLSTM
include("../../shared/load.jl")
using ADPerfTest
using LSTMData
using Zygote

export get_lstm_test

# objective
function sigmoid(x::Float64)
    1. / (1. + exp(-x))
end

function lstmmodel(weight::T1, bias::T2, hidden::T3, cell::T4, input::T5)::Tuple{Vector{Float64}, Vector{Float64}} where {T1<:AbstractVector{Float64}, T2<:AbstractVector{Float64}, T3<:AbstractVector{Float64}, T4<:AbstractVector{Float64}, T5<:AbstractVector{Float64}}
    hsize = size(hidden, 1)
    forget = sigmoid.(input .* view(weight, 1:hsize) .+ view(bias, 1:hsize))
    ingate = sigmoid.(hidden .* view(weight, hsize+1:2hsize) .+ view(bias, hsize+1:2hsize))
    outgate = sigmoid.(input .* view(weight, 2hsize+1:3hsize) .+ view(bias, 2hsize+1:3hsize))
    change = tanh.(hidden .* view(weight, 3hsize+1:4hsize) .+ view(bias, 3hsize+1:4hsize))

    cell2 = cell .* forget .+ ingate .* change
    hidden2 = outgate .* tanh.(cell2)
    hidden2, cell2
end

function lstmpredict(main_params::Matrix{Float64}, extra_params::Matrix{Float64}, state::Matrix{Float64}, input::T)::Tuple{Vector{Float64}, Matrix{Float64}} where {T<:AbstractVector{Float64}}
    x1 = input .* view(extra_params, :, 1)
    x = view(x1, :)
    lenstate = size(state, 2)
    s2 = Matrix{Float64}(undef, size(input, 1), 0)
    for i ∈ 1:2:lenstate
        h, c = lstmmodel(view(main_params, :, i), view(main_params, :, i + 1), view(state, :, i), view(state, :, i + 1), x)
        x = h
        # Zygote does not support mutating arrays
        # TODO: rewrite as array comprehension - initial attempts
        # were not successfully differentiated by Zygote
        s2 = hcat(s2, h, c)
    end
    (x .* view(extra_params, :, 2) .+ view(extra_params, :, 3), s2)
end

function lstmobjective(main_params::Matrix{Float64}, extra_params::Matrix{Float64}, state::Matrix{Float64}, sequence::Matrix{Float64})::Float64
    total = 0.
    count = 0
    input = view(sequence, :, 1)
    b, c = size(sequence)
    curstate = state
    for t ∈ 1:c-1
        ypred, curstate = lstmpredict(main_params, extra_params, curstate, input)
        ynorm = ypred .- log(sum(exp, ypred) + 2.)
        ygold = view(sequence, :, t + 1)
        total += sum(ygold .* ynorm)
        count += b
        input = ygold
    end
    -total / count
end

function pack_zygote_gradient(main_params_grad::Matrix{Float64}, extra_params_grad::Matrix{Float64})
    vcat(reshape(main_params_grad, :), reshape(extra_params_grad, :))
end

mutable struct ZygoteLSTMContext
    input::Union{LSTMInput, Nothing}
    wrapper_lstm_objective::Union{Function, Nothing}
    objective::Float64
    zygote_gradient::Tuple{Matrix{Float64}, Matrix{Float64}}
end

function zygote_lstm_prepare!(ctx::ZygoteLSTMContext, input::LSTMInput)
    # Running Zygote.gradient with test input, because the first invocation
    # of gradient on a given function is very long.
    # Using test input ensures that all computations related to the actual input
    # are done in calculate_jacobian!
    testmainparams = [1. 1.; 1. 1.; 1. 1.; 1. 1.]
    testextraparams = [1. 1. 1.]
    teststate = [1. 1.]
    testsequence = [1. 1.]
    testwrapper_lstm_objective = (main_params, extra_params) -> lstmobjective(main_params, extra_params, teststate, testsequence)
    Zygote.gradient(testwrapper_lstm_objective, testmainparams, testextraparams)

    ctx.input = input
    ctx.wrapper_lstm_objective = (main_params, extra_params) -> lstmobjective(main_params, extra_params, input.state, input.sequence)
end

function zygote_lstm_calculate_objective!(ctx::ZygoteLSTMContext, times)
    for i in 1:times
        ctx.objective = lstmobjective(ctx.input.main_params, ctx.input.extra_params, ctx.input.state, ctx.input.sequence)
    end
end

function zygote_lstm_calculate_jacobian!(ctx::ZygoteLSTMContext, times)
    for i in 1:times
        ctx.zygote_gradient = Zygote.gradient(ctx.wrapper_lstm_objective, ctx.input.main_params, ctx.input.extra_params)
    end
end

function zygote_lstm_output!(out::LSTMOutput, ctx::ZygoteLSTMContext)
    out.objective = ctx.objective
    out.gradient = pack_zygote_gradient(ctx.zygote_gradient...)
end

get_lstm_test() = Test{LSTMInput, LSTMOutput}(
    ZygoteLSTMContext(nothing, nothing, 0, (Matrix{Float64}(undef, 0, 0), Matrix{Float64}(undef, 0, 0))),
    zygote_lstm_prepare!,
    zygote_lstm_calculate_objective!,
    zygote_lstm_calculate_jacobian!,
    zygote_lstm_output!
)

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, Langevin, LinearAlgebra, Plots
using ProgressMeter: @showprogress

begin
    x0 = 0
    x1 = 1
    ncell = 100
    nface = ncell + 1
    dx = (x1 - x0) / ncell
    deg = 2 # polynomial degree
    nsp = deg + 1
    γ = 5 / 3
    cfl = 0.05
    dt = cfl * dx
    t = 0.0
end

ps = FRPSpace1D(x0, x1, ncell, deg)

begin
    uqMethod = "galerkin"
    nr = 5
    nRec = 10
    opType = "uniform"
    parameter1 = 0.95
    parameter2 = 1.05
end

uq = UQ1D(nr, nRec, parameter1, parameter2, opType, uqMethod)

u = zeros(ncell, nsp, 3, uq.nm+1)
for i = 1:ncell, j = 1:nsp
    prim = zeros(3, uq.nm+1)
    if ps.x[i] <= 0.5
        prim[1, :] .= uq.pce
        prim[2, 1] = 0.0
        prim[3, 1] = 0.5
    else
        prim[:, 1] .= [0.3, 0.0, 0.625]
    end

    u[i, j, :, :] .= uq_prim_conserve(prim, γ, uq)
end

function dudt!(du, u, p, t)
    du .= 0.0
    J, ll, lr, lpdm, dgl, dgr, γ, uq = p

    nm = uq.nm
    nq = uq.nq

    ncell = size(u, 1)
    nsp = size(u, 2)

    u_ran = zeros(ncell, nsp, 3, nq)
    for i = 1:ncell, j = 1:nsp, k = 1:3
        u_ran[i, j, k, :] .= chaos_ran(u[i, j, k, :], uq)
    end

    f = zeros(ncell, nsp, 3, nm+1)
    for i = 1:ncell, j = 1:nsp
        _f = zeros(3, nq)
        for k = 1:nq
            _f[:, k] .= euler_flux(u_ran[i, j, :, k], γ)[1] ./ J[i]
        end

        for k = 1:3
            f[i, j, k, :] .= ran_chaos(_f[k, :], uq)
        end
    end

    f_face = zeros(ncell, 3, nm+1, 2)
    for i = 1:ncell, j = 1:3, k = 1:nm+1
        f_face[i, j, k, 1] = dot(f[i, :, j, k], lr)
        f_face[i, j, k, 2] = dot(f[i, :, j, k], ll)
    end

    u_face = zeros(ncell, 3, nq, 2)
    for i = 1:ncell, j = 1:3, k = 1:nq
        u_face[i, j, k, 1] = dot(u_ran[i, :, j, k], lr)
        u_face[i, j, k, 2] = dot(u_ran[i, :, j, k], ll)
    end

    fq_interaction = zeros(ncell + 1, 3, nq)
    for i = 2:ncell, j = 1:nq
        fw = @view fq_interaction[i, :, j]
        flux_hll!(fw, u_face[i-1, :, j, 1], u_face[i, :, j, 2], γ, 1.0)
    end

    f_interaction = zeros(ncell + 1, 3, nm+1)
    for i = 2:ncell, j = 1:3
        f_interaction[i, j, :] .= ran_chaos(fq_interaction[i, j, :], uq)
    end

    rhs1 = zero(u)
    for i = 1:ncell, ppp1 = 1:nsp, k = 1:3, l = 1:nm+1
        rhs1[i, ppp1, k, l] = dot(f[i, :, k, l], lpdm[ppp1, :])
    end

    idx = 2:ncell-1 # ending points are Dirichlet
    for i in idx, ppp1 = 1:nsp, k = 1:3, l = 1:nm+1
        du[i, ppp1, k, l] =
            -(
                rhs1[i, ppp1, k, l] +
                (f_interaction[i, k, l] / J[i] - f_face[i, k, l, 2]) * dgl[ppp1] +
                (f_interaction[i+1, k, l] / J[i] - f_face[i, k, l, 1]) * dgr[ppp1]
            )
    end
end

tspan = (0.0, 0.15)
p = (ps.J, ps.ll, ps.lr, ps.dl, ps.dhl, ps.dhr, γ, uq)
prob = ODEProblem(dudt!, u, tspan, p)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

@showprogress for iter = 1:nt
    step!(itg)
end

sol = zeros(ncell, nsp, 3, 2)
for i in axes(sol, 1), j in axes(sol, 2)
    p1 = zeros(3, uq.nm+1)
    p1 = uq_conserve_prim(itg.u[i, j, :, :], γ, uq)
    p1[end, :] .= lambda_tchaos(p1[end, :], 1.0, uq)

    for k = 1:3
        sol[i, j, k, 1] = mean(p1[k, :], uq.op)
        sol[i, j, k, 2] = std(p1[k, :], uq.op)
    end
end

plot(ps.x, sol[:, 2, 1, 1], label="mean", xlabel="x", ylabel="ρ")
plot!(ps.x, sol[:, 2, 1, 1] .+ sol[:, 2, 1, 2], label="mean+std")
plot!(ps.x, sol[:, 2, 1, 1] .- sol[:, 2, 1, 2], label="mean-std")

cd(@__DIR__)
savefig("rho.png")

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG