text:
struct ResidualBlock
  conv_layers
  norm_layers
  shortcut
end

@treelike ResidualBlock

function ResidualBlock(filters, kernels::Array{Tuple{Int,Int}}, pads::Array{Tuple{Int,Int}}, strides::Array{Tuple{Int,Int}}, shortcut = identity)
  local conv_layers = []
  local norm_layers = []
  for i in 2:length(filters)
    push!(conv_layers, Conv(kernels[i-1], filters[i-1]=>filters[i], pad = pads[i-1], stride = strides[i-1]))
    push!(norm_layers, BatchNorm(filters[i]))
  end
  ResidualBlock(Tuple(conv_layers),Tuple(norm_layers),shortcut)
end

function ResidualBlock(filters, kernels::Array{Int}, pads::Array{Int}, strides::Array{Int}, shortcut = identity)
  ResidualBlock(filters, [(i,i) for i in kernels], [(i,i) for i in pads], [(i,i) for i in strides], shortcut)
end

function (block::ResidualBlock)(input)
  local value = copy.(input)
  for i in 1:length(block.conv_layers)-1
    value = relu.((block.norm_layers[i])((block.conv_layers[i])(value)))
  end
  relu.(((block.norm_layers[end])((block.conv_layers[end])(value))) + block.shortcut(input))
end

function Bottleneck(filters::Int, downsample::Bool = false, res_top::Bool = false)
  if(!downsample && !res_top)
    return ResidualBlock([4 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1])
  elseif(downsample && res_top)
    return ResidualBlock([filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1], Chain(Conv((1,1), filters=>4 * filters, pad = (0,0), stride = (1,1)), BatchNorm(4 * filters)))
  else
    shortcut = Chain(Conv((1,1), 2 * filters=>4 * filters, pad = (0,0), stride = (2,2)), BatchNorm(4 * filters))
    return ResidualBlock([2 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,2], shortcut)
  end
end

function resnet50()
  local layers = [3, 4, 6, 3]
  local layer_arr = []

  push!(layer_arr, Conv((7,7), 3=>64, pad = (3,3), stride = (2,2)))
  push!(layer_arr, MaxPool((3,3), pad = (1,1), stride = (2,2)))

  initial_filters = 64
  for i in 1:length(layers)
    push!(layer_arr, Bottleneck(initial_filters, true, i==1))
    for j in 2:layers[i]
      push!(layer_arr, Bottleneck(initial_filters))
    end
    initial_filters *= 2
  end

  push!(layer_arr, MeanPool((7,7)))
  push!(layer_arr, x -> reshape(x, :, size(x,4)))
  push!(layer_arr, (Dense(2048, 1000)))
  push!(layer_arr, softmax)

  Chain(layer_arr...)
end

function resnet_layers()
  weight = Metalhead.weights("resnet.bson")
  weights = Dict{Any ,Any}()
  for ele in keys(weight)
    weights[string(ele)] = convert(Array{Float64, N} where N, weight[ele])
  end
  ls = resnet50()
  ls[1].weight.data .= weights["gpu_0/conv1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]
  count = 2
  for j in [3:5, 6:9, 10:15, 16:18]
    for p in j
      ls[p].conv_layers[1].weight.data .= weights["gpu_0/res$(count)_$(p-j[1])_branch2a_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]
      ls[p].conv_layers[2].weight.data .= weights["gpu_0/res$(count)_$(p-j[1])_branch2b_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]
      ls[p].conv_layers[3].weight.data .= weights["gpu_0/res$(count)_$(p-j[1])_branch2c_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]
    end
    count += 1
  end
  ls[21].W.data .= transpose(weights["gpu_0/pred_w_0"]); ls[21].b.data .= weights["gpu_0/pred_b_0"]
  Flux.testmode!(ls)
  return ls
end

struct ResNet <: ClassificationModel{ImageNet.ImageNet1k}
  layers::Chain
end

ResNet() = ResNet(resnet_layers())

Base.show(io::IO, ::ResNet) = print(io, "ResNet()")

@treelike ResNet

(m::ResNet)(x) = m.layers(x)


text:
immutable ERBFilterbank{C,G,T<:Real,U<:Real,V<:Real} <: Filterbank
	filters::Vector{SecondOrderSections{C,G}}
	ERB::Vector{T}
  center_frequencies::Vector{U}
  fs::V
end

function make_erb_filterbank(fs, num_channels, low_freq, EarQ = 9.26449, minBW = 24.7, order = 1)
    T = 1/fs
    if length(num_channels) == 1
	cf = erb_space(low_freq, fs/2, num_channels)
    else
	cf = num_channels
	if size(cf,2) > size(cf,1)
	    cf = cf'
	end
    end
    ERB = ((cf/EarQ).^order .+ minBW^order).^(1/order)
    B = 1.019*2*pi*ERB
    B0 = T
    B2 = 0.0
    A0 = 1.0
    A1 = -2*cos(2*cf*pi*T)./exp(B*T)
    A2 = exp(-2*B*T)

    B11 = -(2*T*cos(2*cf*pi*T)./exp(B*T) .+ 2*sqrt(3+2^1.5)*T*sin(2*cf*pi*T)./exp(B*T))/2
    B12 = -(2*T*cos(2*cf*pi*T)./exp(B*T) .- 2*sqrt(3+2^1.5)*T*sin(2*cf*pi*T)./exp(B*T))/2
    B13 = -(2*T*cos(2*cf*pi*T)./exp(B*T) .+ 2*sqrt(3-2^1.5)*T*sin(2*cf*pi*T)./exp(B*T))/2
    B14 = -(2*T*cos(2*cf*pi*T)./exp(B*T) .- 2*sqrt(3-2^1.5)*T*sin(2*cf*pi*T)./exp(B*T))/2

    gain = abs((-2*exp(4*im*cf*pi*T)*T + 2*exp(-(B*T) +
      2*im*cf*pi*T).*T.*(cos(2*cf*pi*T) - sqrt(3 - 2^(3/2))*
      sin(2*cf*pi*T))) .* (-2*exp(4*im*cf*pi*T)*T + 2*exp(-(B*T) +
      2*im*cf*pi*T).*T.* (cos(2*cf*pi*T) + sqrt(3 - 2^(3/2)) *
      sin(2*cf*pi*T))).* (-2*exp(4*im*cf*pi*T)*T + 2*exp(-(B*T) +
      2*im*cf*pi*T).*T.* (cos(2*cf*pi*T) - sqrt(3 +
      2^(3/2))*sin(2*cf*pi*T))) .* (-2*exp(4*im*cf*pi*T)*T +
      2*exp(-(B*T) + 2*im*cf*pi*T).*T.* (cos(2*cf*pi*T) + sqrt(3 +
      2^(3/2))*sin(2*cf*pi*T))) ./ (-2 ./ exp(2*B*T) -
      2*exp(4*im*cf*pi*T) + 2*(1 + exp(4*im*cf*pi*T))./exp(B*T)).^4)

  C = typeof(B0)
	filters = Array(SOSFilter{C,C}, num_channels)
	for ch=1:num_channels
		biquads = Array(BiquadFilter{C}, 4)
		biquads[1] = BiquadFilter(B0, B11[ch], B2, A0, A1[ch], A2[ch])
		biquads[2] = BiquadFilter(B0, B12[ch], B2, A0, A1[ch], A2[ch])
		biquads[3] = BiquadFilter(B0, B13[ch], B2, A0, A1[ch], A2[ch])
		biquads[4] = BiquadFilter(B0, B14[ch], B2, A0, A1[ch], A2[ch])
		filters[ch] = SOSFilter(biquads, 1/gain[ch])
	end
	ERBFilterbank(filters, ERB, cf, fs)
end

function erb_space(low_freq, high_freq, num_channels, EarQ = 9.26449, minBW = 24.7, order = 1)
    # All of the following expressions are derived in Apple TR #35, "An
    # Efficient Implementation of the Patterson-Holdsworth Cochlear
    # Filter Bank."  See pages 33-34.
    cfArray = -(EarQ*minBW) .+ exp([1:num_channels]*(-log(high_freq + EarQ*minBW) + log(low_freq + EarQ*minBW))/num_channels) * (high_freq + EarQ*minBW)
end

ext: jl


text:
# This file is part of Kpax3. License is MIT.

function test_likelihoods_marginal()
  for (α, β) in ([0.1, 0.1], [0.5, 0.5], [1.0, 1.0], [10.0, 10.0], [100.0, 100.0], [0.2, 0.1], [1.0, 0.5], [2.0, 1.0], [20.0, 10.0], [200.0, 100.0], [0.1, 0.2], [0.5, 1.0], [1.0, 2.0], [10.0, 20.0], [100.0, 200.0])
    for (n, y) in ([1.0, 0.0], [1.0, 1.0], [5.0, 0.0], [5.0, 1.0], [5.0,  3.0], [5.0,  5.0], [10.0, 0.0], [10.0, 1.0], [10.0,  5.0], [10.0, 10.0], [100.0, 0.0], [100.0, 1.0], [100.0, 10.0], [100.0, 100.0])
      logp = SpecialFunctions.loggamma(α + y) + SpecialFunctions.loggamma(β + n - y) - SpecialFunctions.loggamma(α + β + n) + SpecialFunctions.loggamma(α + β) - SpecialFunctions.loggamma(α) - SpecialFunctions.loggamma(β)

      logcp0 = log(β + n - y) - log(α + β + n)
      logcp1 = log(α + y) - log(α + β + n)

      @test isapprox(Kpax3.marglik(y, n, α, β), exp(logp), atol=ε)
      @test isapprox(Kpax3.logmarglik(y, n, α, β), logp, atol=ε)

      @test isapprox(Kpax3.condmarglik(0x00, y, n, α, β), exp(logcp0), atol=ε)
      @test isapprox(Kpax3.condmarglik(0x01, y, n, α, β), exp(logcp1), atol=ε)

      @test isapprox(Kpax3.logcondmarglik(0x00, y, n, α, β), logcp0, atol=ε)
      @test isapprox(Kpax3.logcondmarglik(0x01, y, n, α, β), logcp1, atol=ε)
    end
  end

  for (α, β) in ([0.1, 0.1], [0.5, 0.5], [1.0, 1.0], [10.0, 10.0], [100.0, 100.0], [0.2, 0.1], [1.0, 0.5], [2.0, 1.0], [20.0, 10.0], [200.0, 100.0], [0.1, 0.2], [0.5, 1.0], [1.0, 2.0], [10.0, 20.0], [100.0, 200.0])
    for (n, y) in ((1.0, [0.0, 0.0, 1.0, 0.0, 1.0]), (10.0, [2.0, 8.0, 5.0, 10.0, 0.0]), (100.0, [100.0, 2.0, 72.0, 34.0, 0.0]))
      m = length(y)
      n1s = cld(m, 2)

      a = fill(α, m)
      b = fill(β, m)

      q = Kpax3.marglik(y, n, a, b)
      logq = Kpax3.logmarglik(y, n, a, b)

      q0 = Kpax3.condmarglik(fill(0x00, m), y, n, a, b)
      logq0 = Kpax3.logcondmarglik(fill(0x00, m), y, n, a, b)

      q1 = Kpax3.condmarglik(fill(0x01, m), y, n, a, b)
      logq1 = Kpax3.logcondmarglik(fill(0x01, m), y, n, a, b)

      qx = Kpax3.condmarglik([fill(0x01, n1s); fill(0x00, m - n1s)], y, n, a, b)
      logqx = Kpax3.logcondmarglik([fill(0x01, n1s); fill(0x00, m - n1s)], y, n, a, b)

      logp = SpecialFunctions.loggamma.(a .+ y) .+ SpecialFunctions.loggamma.(b .+ n .- y) .- SpecialFunctions.loggamma.(a .+ b .+ n) .+ SpecialFunctions.loggamma.(a .+ b) .- SpecialFunctions.loggamma.(a) .- SpecialFunctions.loggamma.(b)

      logcp0 = log.(b .+ n .- y) .- log.(a .+ b .+ n)
      logcp1 = log.(a .+ y) .- log.(a .+ b .+ n)

      logcpx = log.([a[1:n1s] .+ y[1:n1s]; b[(n1s + 1):m] .+ n .- y[(n1s + 1):m]]) - log.(a .+ b .+ n)

      for i in 1:m
        @test isapprox(q[i], exp(logp[i]), atol=ε)
        @test isapprox(logq[i], logp[i], atol=ε)

        @test isapprox(q0[i], exp(logcp0[i]), atol=ε)
        @test isapprox(logq0[i], logcp0[i], atol=ε)

        @test isapprox(q1[i], exp(logcp1[i]), atol=ε)
        @test isapprox(logq1[i], logcp1[i], atol=ε)

        @test isapprox(qx[i], exp(logcpx[i]), atol=ε)
        @test isapprox(logqx[i], logcpx[i], atol=ε)
      end
    end
  end

  nothing
end

test_likelihoods_marginal()

function test_likelihoods_loglik()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [1; 1; 1; 4; 6; 6]

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  copyto!(state.R, R)

  # fill C with harmless wrong values, just in case
  state.C = UInt8[1 1 1 4 2 1 3 3 1 1 1 1 1 1 1 1 1 2;
                  0 0 2 0 2 0 2 0 0 0 0 2 1 1 0 3 0 4;
                  1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 0 0 0;
                  1 1 1 3 2 1 3 4 1 1 1 1 1 1 1 1 1 2;
                  0 2 1 2 0 0 1 0 4 0 0 0 2 3 0 0 0 0;
                  1 1 1 3 2 1 4 3 1 1 1 1 1 1 1 1 1 2]

  state.emptycluster = [false; true; true; false; true; false]
  state.cl = [1; 4; 6; 0; 0; 0]
  state.k = 3

  state.v = [3; 0; 0; 1; 0; 2]

  state.n1s = zeros(Float64, n, m)
  state.n1s[1, :] .= vec(sum(float(data[:, R .== 1]), dims=2))
  state.n1s[4, :] .= vec(sum(float(data[:, R .== 4]), dims=2))
  state.n1s[6, :] .= vec(sum(float(data[:, R .== 6]), dims=2))

  state.unit = Vector{Int}[zeros(Int, n) for g in 1:n]
  state.unit[1] = [1; 2; 3; 0; 0; 0]
  state.unit[4] = [4; 0; 0; 0; 0; 0]
  state.unit[6] = [5; 6; 0; 0; 0; 0]

  support = Kpax3.MCMCSupport(state, priorC);

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  loglik = zeros(Float64, 3)

  linearidx = Int[state.C[1, b] + 4 * (b - 1) for b in 1:m]
  loglik[1] = sum(Kpax3.logmarglik(vec(state.n1s[1, :]), state.v[1], priorC.A[linearidx], priorC.B[linearidx]))

  linearidx = Int[state.C[4, b] + 4 * (b - 1) for b in 1:m]
  loglik[2] = sum(Kpax3.logmarglik(vec(state.n1s[4, :]), state.v[4], priorC.A[linearidx], priorC.B[linearidx]))

  linearidx = Int[state.C[6, b] + 4 * (b - 1) for b in 1:m]
  loglik[3] = sum(Kpax3.logmarglik(vec(state.n1s[6, :]), state.v[6], priorC.A[linearidx], priorC.B[linearidx]))

  ll = loglik[1] + loglik[2] + loglik[3]

  @test isapprox(Kpax3.loglikelihood(state.C, state.cl, state.k, state.v, state.n1s, priorC), ll, atol=ε)
  @test isapprox(Kpax3.loglikelihood(state.C, state.cl, state.k, support), ll, atol=ε)

  nothing
end

test_likelihoods_loglik()

function test_likelihoods_logmarginal()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [3; 3; 1; 1; 5; 5]
  g = sort(unique(R))
  k = length(g)
  len = maximum(g)

  C = ones(UInt8, len, m)

  emptycluster = trues(len)
  emptycluster[g] .= false

  cl = zeros(Int, len)
  cl[1:k] .= g

  v = zeros(Int, len)
  v[g] .= [2; 2; 2]

  n1s = zeros(Float64, len, m)
  n1s[1, :] .= vec(sum(float(data[:, R .== 1]), dims=2))
  n1s[3, :] .= vec(sum(float(data[:, R .== 3]), dims=2))
  n1s[5, :] .= vec(sum(float(data[:, R .== 5]), dims=2))

  unit = Vector{Int}[sum(R .== g) > 0 ? findall(R .== g) : [0] for g in 1:n]

  logpR = -6.5792512120101012129680384532548487186431884765625
  logpC = [-9.1948612277878307708078864379785954952239990234375;
           -7.2468962917275181467857692041434347629547119140625]
  loglik = -66.7559125873377894322402426041662693023681640625
  logpp = logpR + logpC[1] + loglik

  state = Kpax3.AminoAcidState(copy(R), copy(C), copy(emptycluster), copy(cl), k, copy(v), copy(n1s), deepcopy(unit), logpR, copy(logpC), loglik, logpp)

  support = Kpax3.MCMCSupport(state, priorC)

  logmlik = Kpax3.logmarglikelihood(state.cl, state.k, support.lp, priorC)

  lml = 0.0
  tmp = zeros(Float64, 3)
  for b in 1:m
    tmp[1] = priorC.logγ[1]
    tmp[2] = priorC.logγ[2]
    tmp[3] = priorC.logγ[3]

    for l in 1:k
      g = cl[l]

      tmp[1] += support.lp[1, g, b]
      tmp[2] += support.lp[2, g, b]
      tmp[3] += log(exp(priorC.logω[k][1] + support.lp[3, g, b]) + exp(priorC.logω[k][2] + support.lp[4, g, b]))
    end

    lml += log(exp(tmp[1]) + exp(tmp[2]) + exp(tmp[3]))
  end

  @test isapprox(logmlik, lml, atol=ε)

  nothing
end

test_likelihoods_logmarginal()


text:
using SimulatedNeuralMoments, Flux, MCMCChains, StatsPlots, DelimitedFiles
using BSON:@save
using BSON:@load

# get the things to define the structure for the model
include("JDlib.jl") # contains the functions for the DSGE model
function main()
lb, ub = PriorSupport()

# fill in the structure that defines the model
model = SNMmodel("JD example", lb, ub, InSupport, Prior, PriorDraw, auxstat)

# Here, you can train the net from scratch, or use a previous run
# train the net, and save it and the transformation info
#nnmodel, nninfo = MakeNeuralMoments(model)
#@save "neuralmodel.bson" nnmodel nninfo  # use this line to save the trained neural net
@load "neuralmodel.bson" nnmodel nninfo # use this to load a trained net

# draw a sample at the design parameters, from the prior, or use the official "real" data
data = JDmodel(TrueParameters(), 1, 10)
@show size(data)
#=
# define the neural moments using the data
m = NeuralMoments(auxstat(data), model, nnmodel, nninfo)
# Here, you can create a new chain, or use the results from a previous run
# draw a chain of length 10000 plus 500 burnin
chain, junk, junk = MCMC(m, 10500, model, nnmodel, nninfo, covreps = 100, verbosity=true, do_cue = false)
chain = chain[501:end,:]
writedlm("chain.txt", chain)
#chain = readdlm("chain.txt")

# visualize results
chn = Chains(chain, ["ϕ", "ρ", "σ"])
plot(chn)
savefig("chain.png")
display(chn)
=#
end
main()

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
using DifferentialEquations, ReservoirComputing, Surrogates, ModelingToolkit, Plots, MLJLinearModels, SparseArrays, ProgressBars, CellMLToolkit, Sundials

model_root = joinpath(splitdir(pathof(CellMLToolkit))[1], "..", "models")
ml = CellModel(joinpath(model_root, "ohara_rudy_cipa_v1_2017.cellml.xml"))

tspan = (0, 2000.0)
prob = ODEProblem(ml, tspan);
sol = solve(prob, Rosenbrock23(autodiff=false), dtmax=0.5, save_idxs=49)
plot(sol, vars=1)

#p = list_params(ml)
#update_list!(p, :membrane₊i_Stim_Period, 160.0)
#prob = ODEProblem(ml, tspan; p=p)
p = prob.p; p[201] = 300.0; prob = remake(prob, p=p)
tt = collect(tspan[1]:1.0:tspan[2])
sol = solve(prob, Rosenbrock23(autodiff=false), dtmax=0.5, save_idxs=49:49, saveat=tt)
plot!(sol, vars=1)

modelODESize = 1 #length(prob.u0) # 1 => just fit one component; length(prob.u0) => fit em all

reservoirSize = 300
iterative = reservoirSize > 300 # this is kind of machine-specific?

W   = ReservoirComputing.init_reservoir_givendeg(reservoirSize, 0.8, 6)
Win = ReservoirComputing.init_dense_input_layer(reservoirSize, modelODESize, 0.1)
if iterative
	W = sparse(W)
end

function CTESN!(du, u, p, t)
	du .= p[1].*tanh.(p[2]*u .+ p[3]*p[4](t; idxs=1:1)) .- (1.0 - p[1]).*u # p[3](t; idxs=49:49)
	return nothing
end
esnprob = ODEProblem(CTESN!, randn(Float64, reservoirSize), tspan, (0.8, W, Win, sol))
r = solve(esnprob; saveat = sol.t)

@assert isapprox(sol.t, r.t)

# compute first Wout
Wout = zeros(Float64, modelODESize, reservoirSize)
function fitData!(Wout, xData, rData; beta = 1e-3, iterative = iterative)
	fitModel = RidgeRegression(beta)
	#=
	Note:
		rData = hcat(r(sol.t)...)
		xData = hcat(sol(sol.t)...)
	But we avoid re-interpolating them every time by requiring the r solution use `saveat=sol.t`
	=#
	for n=1:modelODESize
		tmp = transpose(fit(fitModel, transpose(rData[:,:]), xData[n,:]; solver=Analytical(iterative=iterative)))
		Wout[n,:] .= tmp[1:reservoirSize]
	end
	return nothing
end
fitData!(Wout, sol[1:1,:], r[:,:]; beta=1e-3)

# define lower and upper bounds for each variable
#u_lower = minimum(sol, dims=2)[49:49]
#u_upper = maximum(sol, dims=2)[49:49]
p_lower = (40.0)
p_upper = (1000.0)

# select set of Sobol-sampled state vectors
P	= sample(200, p_lower, p_upper, SobolSample()) # this generates tuples; convert to arrays in prob_func
Wout 	= [zeros(Float64, modelODESize, reservoirSize) for _ in 1:length(P)]

# evaluate the ESN at each p, returning WOut(p)
function fitSamples!(Wout, P; p0=p, beta=1e-3)

	function prob_func(prob, i, repeat)
		p0[201] = P[i][1]
		println(i)
		return remake(prob, p=p0)
	end
	ensprob = EnsembleProblem(prob, prob_func = prob_func)
	sim = solve(ensprob, Rosenbrock23(autodiff=false), EnsembleThreads(); trajectories=length(P), dtmax=0.5, save_idxs=49:49, saveat=tt)
	for n in ProgressBar(1:length(P))
		#print("Fitting u[$(n)]...   \t")
		fitData!(Wout[n], sim[n][1:1,:], r; beta = beta)
		#print("Done.\n")
	end
	return nothing
end
fitSamples!(Wout, P)

# form interpolant over sampled parameter vectors
WoutInterpolant = RadialBasis(P, Wout, p_lower, p_upper)

# define unseen p̂ ensure it is not in p
p̂ = (125.0)
# form Wout(p̂) and approximate x(t;p̂) from WoutInterpolant at unseen p̂
Ŵout = reshape(WoutInterpolant(p̂), size(Wout[1]))
x = Ŵout*r
#xInterpolant = RadialBasis(collect(transpose(sol.t)), x, [tspan[begin]], [tspan[end]])
p0 = p
p0[201] = p̂[1]
sol_hat = solve(remake(prob, p=p0), Rosenbrock23(autodiff=false); dtmax=0.5, save_idxs=49:49, saveat=sol.t)

plt = plot(sol, vars=(0,1))
plot!(plt, sol_hat, vars=(0,1), linestyle=:dash, linewidth=2)
plot!(plt, sol.t, x[1,:], linestyle=:dot, linewidth=3)
plot!(legend=false)
#if n==2; plot!(plt, yscale=:log); end
plot(plt)
savefig("./ohara_rudy_CTESN.svg")

max_stars_repo_name: cmarcotte/CTESN

max_issues_repo_name: cmarcotte/CTESN

max_forks_repo_name: cmarcotte/CTESN


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineAgentInstanceView <: SwaggerModel
    vmAgentVersion::Any # spec type: Union{ Nothing, String } # spec name: vmAgentVersion
    extensionHandlers::Any # spec type: Union{ Nothing, Vector{VirtualMachineExtensionHandlerInstanceView} } # spec name: extensionHandlers
    statuses::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: statuses

    function VirtualMachineAgentInstanceView(;vmAgentVersion=nothing, extensionHandlers=nothing, statuses=nothing)
        o = new()
        validate_property(VirtualMachineAgentInstanceView, Symbol("vmAgentVersion"), vmAgentVersion)
        setfield!(o, Symbol("vmAgentVersion"), vmAgentVersion)
        validate_property(VirtualMachineAgentInstanceView, Symbol("extensionHandlers"), extensionHandlers)
        setfield!(o, Symbol("extensionHandlers"), extensionHandlers)
        validate_property(VirtualMachineAgentInstanceView, Symbol("statuses"), statuses)
        setfield!(o, Symbol("statuses"), statuses)
        o
    end
end # type VirtualMachineAgentInstanceView

const _property_map_VirtualMachineAgentInstanceView = Dict{Symbol,Symbol}(Symbol("vmAgentVersion")=>Symbol("vmAgentVersion"), Symbol("extensionHandlers")=>Symbol("extensionHandlers"), Symbol("statuses")=>Symbol("statuses"))
const _property_types_VirtualMachineAgentInstanceView = Dict{Symbol,String}(Symbol("vmAgentVersion")=>"String", Symbol("extensionHandlers")=>"Vector{VirtualMachineExtensionHandlerInstanceView}", Symbol("statuses")=>"Vector{InstanceViewStatus}")
Base.propertynames(::Type{ VirtualMachineAgentInstanceView }) = collect(keys(_property_map_VirtualMachineAgentInstanceView))
Swagger.property_type(::Type{ VirtualMachineAgentInstanceView }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineAgentInstanceView[name]))}
Swagger.field_name(::Type{ VirtualMachineAgentInstanceView }, property_name::Symbol) =  _property_map_VirtualMachineAgentInstanceView[property_name]

function check_required(o::VirtualMachineAgentInstanceView)
    true
end

function validate_property(::Type{ VirtualMachineAgentInstanceView }, name::Symbol, val)
end


text:
struct Data
    data::Vector
end

Base.length(data::Data) = length(data.data)

function Merlin.todevice(data::Data, dev::Int)
    dev < 0 && return data
    data = map(data.data) do (w,c,dims_c,t)
        w = todevice(w, dev)
        c = todevice(c, dev)
        t = todevice(t, dev)
        (w, c, dims_c, t)
    end
    Data(data)
end

function Base.getindex(data::Data, indexes::Vector{Int})
    data = data.data[indexes]
    data = sort(data, by=x->length(x[1]), rev=true)
    ws = map(x -> x[1], data)
    dims_w = length.(ws)
    w = cat(ws..., dims=1)
    w = reshape(w, 1, length(w))
    c = cat(map(x -> x[2], data)..., dims=1)
    c = reshape(c, 1, length(c))
    dims_c = cat(map(x -> x[3], data)..., dims=1)
    t = cat(map(x -> x[4], data)..., dims=1)
    (w=Var(w), c=Var(c), dims_w=dims_w, dims_c=dims_c, t=Var(t))
end

function readconll(path::String, dicts)
    data = []
    words = String[]
    tagids = Int[]
    unkword = dicts.w["UNKNOWN"]
    unkchar = dicts.c["UNKNOWN"]

    lines = open(readlines, path)
    push!(lines, "")
    for i = 1:length(lines)
        line = lines[i]
        if isempty(line)
            isempty(words) && continue
            wordids = Int[]
            charids = Int[]
            chardims = Int[]
            for w in words
                # w0 = replace(lowercase(w), r"[0-9]", '0')
                id = get(dicts.w, lowercase(w), unkword)
                push!(wordids, id)

                chars = Vector{Char}(w)
                push!(chardims, length(chars))
                cids = map(chars) do c
                    get(dicts.c, string(c), unkchar)
                end
                append!(charids, cids)
            end
            push!(data, (wordids,charids,chardims,tagids))
            words = String[]
            tagids = Int[]
        else
            items = Vector{String}(split(line,"\t"))
            word = strip(items[1])
            @assert !isempty(word)
            push!(words, word)
            if length(items) >= 2
                tag = strip(items[2])
                push!(tagids, dicts.t[tag])
            else
                throw("")
            end
        end
    end
    data = Vector{typeof(data[1])}(data)
    Data(data)
end

function initvocab(path::String)
    chardict = Dict{String,Int}()
    tagdict = Dict{String,Int}()
    lines = open(readlines, path)
    for line in lines
        isempty(line) && continue
        items = Vector{String}(split(line,"\t"))
        word = strip(items[1])
        chars = Vector{Char}(word)
        for c in chars
            c = string(c)
            if haskey(chardict, c)
                chardict[c] += 1
            else
                chardict[c] = 1
            end
        end

        tag = strip(items[2])
        haskey(tagdict,tag) || (tagdict[tag] = length(tagdict)+1)
    end

    chars = String[]
    for (k,v) in chardict
        v >= 3 && push!(chars,k)
    end
    chardict = Dict(chars[i] => i for i=1:length(chars))
    chardict["UNKNOWN"] = length(chardict) + 1
    chardict, tagdict
end


text:
function trimap(i::Integer, j::Integer)
    if i < j
        trimap(j, i)
    else
        div((i - 1) * i, 2) + j
    end
end

"""
    extract_eigenvalues(model, f::MOI.VectorAffineFunction{T}, d::Int, offset::Int) where T

The vector `f` contains `t` (if `offset = 1`) or `(t, u)` (if `offset = 2`)
followed by the matrix `X` of dimension `d`.
This functions extracts the eigenvalues of `X` and returns a vector containing `t` or `(t, u)`,
a vector `MOI.VariableIndex` containing the eigenvalues of `X`,
the variables created and the index of the constraint created to extract the eigenvalues.
"""
function extract_eigenvalues(
    model,
    f::MOI.VectorAffineFunction{T},
    d::Int,
    offset::Int,
) where {T}
    f_scalars = MOIU.eachscalar(f)
    tu = [f_scalars[i] for i in 1:offset]

    n = trimap(d, d)
    X = f_scalars[offset.+(1:n)]
    m = length(X.terms)
    M = m + n + d

    terms = Vector{MOI.VectorAffineTerm{T}}(undef, M)
    terms[1:m] = X.terms
    N = trimap(2d, 2d)
    constant = zeros(T, N)
    constant[1:n] = X.constants

    Δ = MOI.add_variables(model, n)

    cur = m
    for j in 1:d
        for i in j:d
            cur += 1
            terms[cur] = MOI.VectorAffineTerm(
                trimap(i, d + j),
                MOI.ScalarAffineTerm(one(T), Δ[trimap(i, j)]),
            )
        end
        cur += 1
        terms[cur] = MOI.VectorAffineTerm(
            trimap(d + j, d + j),
            MOI.ScalarAffineTerm(one(T), Δ[trimap(j, j)]),
        )
    end
    @assert cur == M
    Y = MOI.VectorAffineFunction(terms, constant)
    sdindex =
        MOI.add_constraint(model, Y, MOI.PositiveSemidefiniteConeTriangle(2d))

    D = Δ[trimap.(1:d, 1:d)]

    return tu, D, Δ, sdindex
end

"""
    LogDetBridge{T}

The `LogDetConeTriangle` is representable by a `PositiveSemidefiniteConeTriangle` and `ExponentialCone` constraints.
Indeed, ``\\log\\det(X) = \\log(\\delta_1) + \\cdots + \\log(\\delta_n)`` where ``\\delta_1``, ..., ``\\delta_n`` are the eigenvalues of ``X``.
Adapting the method from [1, p. 149], we see that ``t \\le u \\log(\\det(X/u))`` for ``u > 0`` if and only if there exists a lower triangular matrix ``Δ`` such that
```math
\\begin{align*}
  \\begin{pmatrix}
    X & Δ\\\\
    Δ^\\top & \\mathrm{Diag}(Δ)
  \\end{pmatrix} & \\succeq 0\\\\
  t & \\le u \\log(Δ_{11}/u) + u \\log(Δ_{22}/u) + \\cdots + u \\log(Δ_{nn}/u)
\\end{align*}
```

[1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex optimization: analysis, algorithms, and engineering applications*. Society for Industrial and Applied Mathematics, 2001.
```
"""
struct LogDetBridge{T} <: AbstractBridge
    Δ::Vector{MOI.VariableIndex}
    l::Vector{MOI.VariableIndex}
    sdindex::CI{
        MOI.VectorAffineFunction{T},
        MOI.PositiveSemidefiniteConeTriangle,
    }
    lcindex::Vector{CI{MOI.VectorAffineFunction{T},MOI.ExponentialCone}}
    tlindex::CI{MOI.ScalarAffineFunction{T},MOI.LessThan{T}}
end
function bridge_constraint(
    ::Type{LogDetBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::MOI.LogDetConeTriangle,
) where {T}
    return bridge_constraint(
        LogDetBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end
function bridge_constraint(
    ::Type{LogDetBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::MOI.LogDetConeTriangle,
) where {T}
    d = s.side_dimension
    tu, D, Δ, sdindex = extract_eigenvalues(model, f, d, 2)
    t, u = tu
    l = MOI.add_variables(model, d)
    lcindex = [sublog(model, l[i], u, D[i], T) for i in eachindex(l)]
    tlindex = subsum(model, t, l, T)
    return LogDetBridge(Δ, l, sdindex, lcindex, tlindex)
end

function MOI.supports_constraint(
    ::Type{LogDetBridge{T}},
    ::Type{<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction{T}}},
    ::Type{MOI.LogDetConeTriangle},
) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:LogDetBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{LogDetBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle),
        (MOI.VectorAffineFunction{T}, MOI.ExponentialCone),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
    ]
end

"""
    sublog(model, x::MOI.VariableIndex, y::MOI.VariableIndex, z::MOI.VariableIndex, ::Type{T}) where T

Constrains ``x \\le y \\log(z/y)`` and returns the constraint index.
"""
function sublog(
    model,
    x::MOI.VariableIndex,
    y::MOI.ScalarAffineFunction{T},
    z::MOI.VariableIndex,
    ::Type{T},
) where {T}
    return MOI.add_constraint(
        model,
        MOIU.operate(vcat, T, MOI.SingleVariable(x), y, MOI.SingleVariable(z)),
        MOI.ExponentialCone(),
    )
end

"""
    subsum(model, t::MOI.ScalarAffineFunction, l::Vector{MOI.VariableIndex}, ::Type{T}) where T

Constrains ``t \\le l_1 + \\cdots + l_n`` where `n` is the length of `l` and returns the constraint index.
"""
function subsum(
    model,
    t::MOI.ScalarAffineFunction,
    l::Vector{MOI.VariableIndex},
    ::Type{T},
) where {T}
    n = length(l)
    f = MOIU.operate!(-, T, t, MOIU.operate(sum, T, l))
    return MOIU.normalize_and_add_constraint(
        model,
        f,
        MOI.LessThan(zero(T)),
        allow_modify_function = true,
    )
end

# Attributes, Bridge acting as a model
MOI.get(b::LogDetBridge, ::MOI.NumberOfVariables) = length(b.Δ) + length(b.l)
MOI.get(b::LogDetBridge, ::MOI.ListOfVariableIndices) = [b.Δ; b.l]
function MOI.get(
    b::LogDetBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T},
        MOI.PositiveSemidefiniteConeTriangle,
    },
) where {T}
    return 1
end
function MOI.get(
    b::LogDetBridge{T},
    ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{T},MOI.ExponentialCone},
) where {T}
    return length(b.lcindex)
end
function MOI.get(
    b::LogDetBridge{T},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.LessThan{T}},
) where {T}
    return 1
end
function MOI.get(
    b::LogDetBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T},
        MOI.PositiveSemidefiniteConeTriangle,
    },
) where {T}
    return [b.sdindex]
end
function MOI.get(
    b::LogDetBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T},
        MOI.ExponentialCone,
    },
) where {T}
    return b.lcindex
end
function MOI.get(
    b::LogDetBridge{T},
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},MOI.LessThan{T}},
) where {T}
    return [b.tlindex]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::LogDetBridge)
    MOI.delete(model, bridge.tlindex)
    MOI.delete(model, bridge.lcindex)
    MOI.delete(model, bridge.sdindex)
    MOI.delete(model, bridge.l)
    return MOI.delete(model, bridge.Δ)
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::LogDetBridge,
)
    d = length(bridge.lcindex)
    Δ = MOI.get(model, MOI.VariablePrimal(), bridge.Δ)
    t =
        MOI.get(model, attr, bridge.tlindex) +
        sum(MOI.get(model, attr, ci)[1] for ci in bridge.lcindex)
    u = MOI.get(model, attr, first(bridge.lcindex))[2]
    x = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t, u, x)
end
# [X Δ; Δ' Diag(Δ)] in PSD
# t - sum(l) >= 0
# (l_i, u, Δ_ii) in Exp
# (t, u, x) in LogDet <=> exists Δ, l such that At + Bu + Cx + DΔ + El in (PSD, >=, Exp_i)
# so LogDet* = [A'; B'; C'] (PSD, >=, Exp_i)*
# and 0 = [D'; E'] (PSD, >=, Exp_i)*
# where
# A = [0, 0, 0, 1, 0, 0, 0]
# B = [0, 0, 0, 0, 0, 1, 0]
# C = [I, 0, 0, 0, 0, 0, 0]
# D = [0, I, I(i=j), 0, 0, 0, I(i=j)]
# E = [0, 0, 0, 1, I, 0, 0]
# so given dual q = (a, b, c, d, e, f, g), we get
# t = A' q = d => t = d
# u = B' q = f => u = sum(f)
# x = C' q = a => x = a
# offdiag(b) = 0
# 0 = D' q = diag(b) + diag(c) + g => g = -diag(b) - diag(c)
# let b = 0, so g_i = -c_i
# 0 = E' q = d + e => d = -e
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::LogDetBridge,
)
    t_dual = MOI.get(model, attr, bridge.tlindex)
    u_dual =
        sum(MOI.get(model, attr, lcindex_i)[2] for lcindex_i in bridge.lcindex)
    x_dual = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t_dual, u_dual, x_dual)
end

"""
    RootDetBridge{T}

The `RootDetConeTriangle` is representable by a `PositiveSemidefiniteConeTriangle` and an `GeometricMeanCone` constraints; see [1, p. 149].
Indeed, ``t \\le \\det(X)^{1/n}`` if and only if there exists a lower triangular matrix ``Δ`` such that
```math
\\begin{align*}
  \\begin{pmatrix}
    X & Δ\\\\
    Δ^\\top & \\mathrm{Diag}(Δ)
  \\end{pmatrix} & \\succeq 0\\\\
  t & \\le (Δ_{11} Δ_{22} \\cdots Δ_{nn})^{1/n}
\\end{align*}
```

[1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex optimization: analysis, algorithms, and engineering applications*. Society for Industrial and Applied Mathematics, 2001.
"""
struct RootDetBridge{T} <: AbstractBridge
    Δ::Vector{MOI.VariableIndex}
    sdindex::CI{
        MOI.VectorAffineFunction{T},
        MOI.PositiveSemidefiniteConeTriangle,
    }
    gmindex::CI{MOI.VectorAffineFunction{T},MOI.GeometricMeanCone}
end
function bridge_constraint(
    ::Type{RootDetBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::MOI.RootDetConeTriangle,
) where {T}
    return bridge_constraint(
        RootDetBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end
function bridge_constraint(
    ::Type{RootDetBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::MOI.RootDetConeTriangle,
) where {T}
    d = s.side_dimension
    tu, D, Δ, sdindex = extract_eigenvalues(model, f, d, 1)
    t = tu[1]
    DF = MOI.VectorAffineFunction{T}(MOI.VectorOfVariables(D))
    gmindex = MOI.add_constraint(
        model,
        MOIU.operate(vcat, T, t, DF),
        MOI.GeometricMeanCone(d + 1),
    )

    return RootDetBridge(Δ, sdindex, gmindex)
end

function MOI.supports_constraint(
    ::Type{RootDetBridge{T}},
    ::Type{<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction{T}}},
    ::Type{MOI.RootDetConeTriangle},
) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:RootDetBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{RootDetBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle),
        (MOI.VectorAffineFunction{T}, MOI.GeometricMeanCone),
    ]
end

# Attributes, Bridge acting as a model
MOI.get(b::RootDetBridge, ::MOI.NumberOfVariables) = length(b.Δ)
MOI.get(b::RootDetBridge, ::MOI.ListOfVariableIndices) = b.Δ
function MOI.get(
    b::RootDetBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T},
        MOI.PositiveSemidefiniteConeTriangle,
    },
) where {T}
    return 1
end
function MOI.get(
    b::RootDetBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T},
        MOI.GeometricMeanCone,
    },
) where {T}
    return 1
end
function MOI.get(
    b::RootDetBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T},
        MOI.PositiveSemidefiniteConeTriangle,
    },
) where {T}
    return [b.sdindex]
end
function MOI.get(
    b::RootDetBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T},
        MOI.GeometricMeanCone,
    },
) where {T}
    return [b.gmindex]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::RootDetBridge)
    MOI.delete(model, bridge.gmindex)
    MOI.delete(model, bridge.sdindex)
    return MOI.delete(model, bridge.Δ)
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::RootDetBridge,
)
    t = MOI.get(model, attr, bridge.gmindex)[1]
    x = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t, x)
end
# (t, x) in RootDet <=> exists Δ such that At + Bx + CΔ in (PSD, GeoMean)
# so RootDet* = [A'; B'] (PSD, GeoMean)*
# and 0 = [C'] (PSD, GeoMean)*
# where
# A = [0, 0, 0, 1, 0]
# B = [I, 0, 0, 0, 0]
# C = [0, I, I(i=j), 0, I(i=j)]
# so given dual q = (a, b, c, d, e), we get
# t = A' q = d => t = d
# x = B' q = a => x = a
# offdiag(b) = 0
# 0 = C' q = diag(b) + diag(c) + e => e = -diag(b) - diag(c)
# let b = 0, so e_i = -c_i
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::RootDetBridge,
)
    t_dual = MOI.get(model, attr, bridge.gmindex)[1]
    x_dual = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t_dual, x_dual)
end


text:
include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_utils.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_curves.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_roadways.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_vehicles.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_trajdatas.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_scenes.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_scene_records.jl"))


text:
#!/usr/bin/julia
using ExperimentalDesign, StatsModels, GLM, DataFrames, Distributions, Random, CSV, StatsBase, Plots, LinearAlgebra

function mse(y, ŷ)
	n = length(y)
	diff = y - ŷ
	acc = dot(diff, diff)

	return acc/n
end

out_file = open("validation.csv", "w")

function output_csv_header()
	header = "Exectution; Training_size; Correct; Correct_parallel; MSE"
	println(out_file, header)
end

function output_csv_line(execution, training_size, correct, correct_parallel, MSE)
	str = string(execution) * "; " * string(training_size) * "; " * string(correct) * "; " * string(correct_parallel) * "; " * string(MSE)

	println(out_file, str)
end

results = CSV.read("output_corei7.csv", delim="; ", DataFrame)

filter!(:Expected_insns => x -> x > 0, results)

# Rename labels to lowercase, as @formula seems to not like uppercase strings.
for name in names(results)
	h = name => lowercase(name)
	rename!(results, replace.(names(results), h))
end

output_csv_header()

all_corrects = Int64[]
all_corrects_parallel = Int64[]
all_mse_mean = Float64[]

for j in 1:100
	train_ratio = 0.2
	train_size = round(Int, train_ratio * nrow(results))
	train_rows = shuffle(1:nrow(results))[1:train_size]
	D = results[train_rows, :]

	D = stack(D, Not([:filename, :functions, :expected_insns, :parallel, :inlined_percentage, :num_partitions]))
	sort!(D, [:expected_insns, :value])

	regression_full = @formula(value ~ ((expected_insns + expected_insns ^ 2) +
												(functions + functions ^ 2 + 1 / functions)) *
												parallel * inlined_percentage)

	screening_model = regression_full
	screening_fit = lm(screening_model, D)

	test_rows = [x for x in 1:nrow(results) if !(x in train_rows)]
	Dtest = results[test_rows, :]

	Dtest[!, :predicted] = predict(screening_fit, Dtest)
	sort!(Dtest, [:filename, :parallel])

	corrects = 0
	parallel_correct = 0
	c = 0
	real_best = Float64[]
	pred_best = Float64[]
	len = Int64(floor(nrow(Dtest)/4)) - 1

	for i in 1:len
		c += 1

		real_1 = Dtest[i*4    , :execution_1]
		real_2 = Dtest[i*4 + 1, :execution_1]
		real_4 = Dtest[i*4 + 2, :execution_1]
		real_8 = Dtest[i*4 + 3, :execution_1]

		predicted_1 = Dtest[i*4    , :predicted]
		predicted_2 = Dtest[i*4 + 1, :predicted]
		predicted_4 = Dtest[i*4 + 2, :predicted]
		predicted_8 = Dtest[i*4 + 3, :predicted]

		real_t      = [real_1, real_2, real_4, real_8]
		predicted_t = [predicted_1, predicted_2, predicted_4, predicted_8]

		argmin_real = argmin(real_t)
		argmin_predicted = argmin(predicted_t)

		if argmin_real == argmin_predicted
			corrects += 1
			parallel_correct += 1
		elseif argmin_real > 1 && argmin_predicted > 1
			parallel_correct += 1
		end

		append!(real_best, real_t[argmin_real])
		append!(pred_best, predicted_t[argmin_predicted])
	end
	best_mse = mse(real_best, pred_best)

	println("Execution " * string(j))
	println("Corrects: " * string(corrects) * " from " * string(c) * ". Ratio = " * string(corrects/c))
	println("Parallel correct: " * string(parallel_correct))
	println("MSE: " * string(best_mse))

	append!(all_corrects, corrects)
	append!(all_corrects_parallel, parallel_correct)
	append!(all_mse_mean, best_mse)

	output_csv_line(j, c, corrects, parallel_correct, best_mse)
end

mean_corrects = mean(all_corrects)
mean_corrects_parallel = mean(all_corrects_parallel)
mean_mse = mean(all_mse_mean)

println("------------------------")
println("Mean corrects: " * string(mean_corrects))
println("Mean corrects parallel: " * string(mean_corrects_parallel))
println("Mean MSE: " * string(mean_mse))

max_stars_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_issues_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_forks_repo_name: giulianobelinassi/LTO-Timing-Analysis


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	using SpecialFunctions
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.24"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	a,b = 0.2, 0.7
	x = 0.75
end

# ╔═╡ cdacaea2-2167-11eb-3ad7-7521db4ca126
begin
	betaAB1 = beta(a,b)
	betaAB2 = (gamma(a)gamma(b))/gamma(a+b)
	betaAB3 = (factorial(a-1)factorial(b-1))/factorial(a+b-1)
	betaPDFAB1 = pdf(Beta(a,b),x)
	betaPDFAB2 = (1/beta(a,b))*x^(a-1) * (1-x)^(b-1)
end

# ╔═╡ cdaceaca-2167-11eb-00b5-5b987d87a6c0
Text("beta($a,$b)    = $betaAB1,\t$betaAB2,\t$betaAB3 ")

# ╔═╡ cdad87e4-2167-11eb-28b7-8b9600d8049d
Text("betaPDF($a,$b) at $x = $betaPDFAB1,\t$betaPDFAB2")

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.24"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═cdacaea2-2167-11eb-3ad7-7521db4ca126
# ╠═cdaceaca-2167-11eb-00b5-5b987d87a6c0
# ╠═cdad87e4-2167-11eb-28b7-8b9600d8049d
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
include("post_process_defs.jl")

OUTDIR = length(ARGS) > 0 ? ARGS[1] : "results/K_MCMC20_L_MCMC5_scale0.1_SEED0/"

post_process(OUTDIR)


text:
module GeneralizedGenerated
using MLStyle
using JuliaVariables
using DataStructures
List = LinkedList

export NGG
export gg, @gg, UnderGlobal, @under_global
export RuntimeFn, closure_conv, mk_function, mkngg, mk_expr
export to_type, from_type, runtime_eval
include("closure_conv.jl")

function mk_function(ex)
    mk_function(@__MODULE__, ex)
end

function mk_function(mod::Module, ex)
    ex = macroexpand(mod, ex)
    ex = simplify_ex(ex)
    ex = solve!(ex)
    fn = closure_conv(mod, ex)
    if !(fn isa RuntimeFn)
        error("Expect an unnamed function expression. ")
    end
    fn
end

"""
process an expression and perform closure conversions for all nested expressions.
"""
function mk_expr(mod::Module, ex)
    ex = macroexpand(mod, ex)
    ex = simplify_ex(ex)
    ex = solve!(ex)
    closure_conv(mod, ex)
end

function mk_function(mod::Module, args, kwargs, body)
    mk_function(mod, Expr(:function, :($(args...), ; $(kwargs...)), body))
end

function mk_function(args, kwargs, body)
    mk_function(Main, args, kwargs, body)
end

function runtime_eval(mod::Module, ex)
    fn_ast = :(function ()
        $ex
    end)
    mk_function(mod, fn_ast)()
end

function runtime_eval(ex)
    runtime_eval(Main, ex)
end

end # module


text:
using BenchmarkTools
using LinearAlgebra

using RootsAndPoles

# Define a parent BenchmarkGroup to contain our suite
const suite = BenchmarkGroup()

# Add some child groups to our benchmark suite.
suite["grpf"] = BenchmarkGroup(["functions"])  # BenchmarkGroup "tags"

simplefcn(z) = (z - 1)*(z - im)^2*(z + 1)^3/(z + im)

function defaultfcn(z)
    f = 1e9
    ϵᵣ = 5 - 2im
    μᵣ = 1 - 2im
    d = 1e-2
    c = 3e8
    ω = 2π*f
    k₀ = ω/c
    cc = ϵᵣ^2*(k₀*d)^2*(ϵᵣ*μᵣ - 1)
    w = ϵᵣ^2*z^2 + z^2*tan(z)^2 - cc
end

function wvgd(z)
    ns = 0.065-4im
    n1 = 1.5835
    nc = 1.0
    d1 = 1.81e-6
    λ₀ = 0.6328e-6
    k₀ = 2π/λ₀
    k₀d1 = k₀*d1
    κ1 = sqrt(n1^2 - z^2)
    γs = sqrt(z^2 - ns^2)
    γc = sqrt(z^2 - nc^2)
    m11 = cos(κ1*k₀d1)
    m12 = im/κ1*sin(κ1*k₀d1)
    m21 = im*κ1*sin(κ1*k₀d1)
    m22 = cos(κ1*k₀d1)
    w = det([1.0    -m11+im*γc*m12
           im*γs  -m21+im*γc*m22])
end

function graphenefunction(z)
    f = 1e12
    c = 299792458
    μ₀ = 4π*1e-7
    ϵ₀ = 1/(μ₀*c^2)

    e = 1.602176565e-19
    kB = 1.3806488e-23
    hk = 1.05457168e-34
    vFe = 1e6
    muc = 0.05*e
    t = 0.135e-12
    T = 300
    ϵᵣ₁ = 1.0
    ϵᵣ₂ = 11.9

    ω = 2π*f
    k₀ = ω/c
    kᵣ₀ = -im*z*k₀

    Slo=-im*e^2*kB*T*log(2+2*cosh(muc/kB/T)) / (π*hk^2*(ω-im/t))

    a = -3*vFe^2*Slo/(4*(ω-im/t)^2)
    b = a/3

    Y1TM = ω*ϵᵣ₁*ϵ₀/sqrt(ϵᵣ₁*k₀^2 - kᵣ₀^2);
    Y2TM = ω*ϵᵣ₂*ϵ₀/sqrt(ϵᵣ₂*k₀^2 - kᵣ₀^2);
    YSTM = Slo + 1*a*kᵣ₀^2 + 1*b*kᵣ₀^2;

    w = (Y1TM + Y2TM + YSTM)*(-Y1TM + Y2TM + YSTM)*(Y1TM - Y2TM + YSTM)*(-Y1TM - Y2TM + YSTM) # four Riemann sheets
end

testfcns = (simplefcn, defaultfcn, wvgd, graphenefunction)

# Add some benchmarks to the "grpf" group
for f in testfcns
    sf = string(f)
    if sf == "simplefcn"
        origcoords = rectangulardomain(complex(-2, -2), complex(2, 2), 0.1)
        suite["grpf"][sf] = @benchmarkable grpf($f, $origcoords)
    elseif sf == "defaultfcn"
        origcoords = rectangulardomain(complex(-2, -2), complex(2, 2), 0.2)
        suite["grpf"][sf] = @benchmarkable grpf($f, $origcoords)
    elseif sf == "wvgd"
        origcoords = rectangulardomain(complex(1.0, -1.0), complex(2.5, 2.5), 0.5)
        suite["grpf"][sf] = @benchmarkable grpf($f, $origcoords)
    elseif sf == "graphenefunction"
        origcoords = rectangulardomain(complex(-100, -100), complex(400, 400), 18)
        suite["grpf"][sf] = @benchmarkable grpf($f, $origcoords)
    end
end

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `suite` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end

# run with a time limit of ~5 second per benchmark
results = run(suite, verbose = true, seconds = 5)

# Appends results into `paramspath`
BenchmarkTools.save(paramspath, results)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(hist::EmpiricalHistogram)
  seriestype --> :step
  seriescolor --> :black
  legend --> false
  framestyle --> :box
  grid --> false

  normalize(hist.hist)
end


text:
#ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля без маркеров.
#РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.

#Доп. условие:  на поле могут находиться внутрение перегородки прямоугольной формы, среди которых
#могут быть и вырожденые прямоугольники (отрезки),эти внутренние перегородки изолированы друг от друга и от внешней рамки.

function putmarker(r,side)
    k = 1
    while k != 0
        if isborder(r,side)
            k = obhod(r,side)
            if k == 0
                break
            end
            putmarker!(r)
        end
        if !isborder(r,side)
            move!(r,side)
            putmarker!(r)
        end
    end
end

function moverobot(r,side)
    while ismarker(r) == true
        if isborder(r,side)
            obhod(r,side)
        end
        if !isborder(r,side)
            move!(r,side)
        end
    end
end

function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end

function next_side(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+3,4))
        end
    end
end

function obhod(r,side)
    k = 0
    while (isborder(r,side) && !isborder(r,next_side(side)))
        move!(r,next_side(side))
        k += 1
    end
    if isborder(r,next_side(side))&&isborder(r,side)
        for i in 1:k
            move!(r,reverse(next_side(side)))
        end
        return 0
    end
    if !isborder(r,side)
        move!(r,side)
    end
    if k != 0
        while isborder(r,reverse(next_side(side)))
            move!(r,side)
        end
        for i in 1:k
            move!(r,reverse(next_side(side)))
        end
    end
    return 1
end

function krest(r)
    for side in (HorizonSide(i) for i=0:3)
        putmarker(r,side) #функция передвигает и ставит марки
        moverobot(r,reverse(side)) #функция перемещает в первоначальное положение
    end
    putmarker!(r)
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
# Issue #6573
srand(0); rand(); x = rand(384);
@test find(x .== rand()) == []

@test rand() != rand()
@test 0.0 <= rand() < 1.0
@test rand(Uint32) >= 0
@test -10 <= rand(-10:-5) <= -5
@test -10 <= rand(-10:5) <= 5
@test minimum([rand(int32(1):int32(7^7)) for i = 1:100000]) > 0
@test(typeof(rand(false:true)) == Bool)

@test length(randn(4, 5)) == 20
@test length(randbool(4, 5)) == 20

@test rand(MersenneTwister()) == 0.8236475079774124
@test rand(MersenneTwister(0)) == 0.8236475079774124
@test rand(MersenneTwister(42)) == 0.5331830160438613
# Try a seed larger than 2^32
@test rand(MersenneTwister(5294967296)) == 0.3498809918210497

# randn
@test randn(MersenneTwister(42)) == -0.5560268761463861
A = zeros(2, 2)
randn!(MersenneTwister(42), A)
@test A == [-0.5560268761463861  0.027155338009193845;
            -0.444383357109696  -0.29948409035891055]

for T in (Int8, Uint8, Int16, Uint16, Int32, Uint32, Int64, Uint64, Int128, Uint128,
          Char, Float16, Float32, Float64, Rational{Int})
    r = rand(convert(T, 97):convert(T, 122))
    @test typeof(r) == T
    @test 97 <= r <= 122
    r = rand(convert(T, 97):convert(T,2):convert(T, 122),2)[1]
    @test typeof(r) == T
    @test 97 <= r <= 122
    @test mod(r,2)==1
end

if sizeof(Int32) < sizeof(Int)
    r = rand(int32(-1):typemax(Int32))
    @test typeof(r) == Int32
    @test -1 <= r <= typemax(Int32)
    @test all([div(0x00010000000000000000,k)*k - 1 == Base.Random.RandIntGen(uint64(1:k)).u for k in 13 .+ int64(2).^(32:62)])
    @test all([div(0x00010000000000000000,k)*k - 1 == Base.Random.RandIntGen(int64(1:k)).u for k in 13 .+ int64(2).^(32:61)])

end

# Test ziggurat tables
ziggurat_table_size = 256
nmantissa           = int64(2)^51 # one bit for the sign
ziggurat_nor_r      = BigFloat("3.65415288536100879635194725185604664812733315920964488827246397029393565706474")
nor_section_area    = ziggurat_nor_r*exp(-ziggurat_nor_r^2/2) + erfc(ziggurat_nor_r/sqrt(BigFloat(2)))*sqrt(big(π)/2)
emantissa           = int64(2)^52
ziggurat_exp_r      = BigFloat("7.69711747013104971404462804811408952334296818528283253278834867283241051210533")
exp_section_area    = (ziggurat_exp_r + 1)*exp(-ziggurat_exp_r)

const ki = Array(Uint64, ziggurat_table_size)
const wi = Array(Float64, ziggurat_table_size)
const fi = Array(Float64, ziggurat_table_size)
# Tables for exponential variates
const ke = Array(Uint64, ziggurat_table_size)
const we = Array(Float64, ziggurat_table_size)
const fe = Array(Float64, ziggurat_table_size)
function randmtzig_fill_ziggurat_tables() # Operates on the global arrays
    wib = big(wi)
    fib = big(fi)
    web = big(we)
    feb = big(fe)
    # Ziggurat tables for the normal distribution
    x1 = ziggurat_nor_r
    wib[256] = x1/nmantissa
    fib[256] = exp(-0.5*x1*x1)
    # Index zero is special for tail strip, where Marsaglia and Tsang
    # defines this as
    # k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
    # where v is the area of each strip of the ziggurat.
    ki[1] = uint64(itrunc(x1*fib[256]/nor_section_area*nmantissa))
    wib[1] = nor_section_area/fib[256]/nmantissa
    fib[1] = one(BigFloat)

    for i = 255:-1:2
        # New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
        # need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
        x = sqrt(-2.0*log(nor_section_area/x1 + fib[i+1]))
        ki[i+1] = uint64(itrunc(x/x1*nmantissa))
        wib[i] = x/nmantissa
        fib[i] = exp(-0.5*x*x)
        x1 = x
    end

    ki[2] = uint64(0)

    # Zigurrat tables for the exponential distribution
    x1 = ziggurat_exp_r
    web[256] = x1/emantissa
    feb[256] = exp(-x1)

    # Index zero is special for tail strip, where Marsaglia and Tsang
    # defines this as
    # k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
    # where v is the area of each strip of the ziggurat.
    ke[1] = uint64(itrunc(x1*feb[256]/exp_section_area*emantissa))
    web[1] = exp_section_area/feb[256]/emantissa
    feb[1] = one(BigFloat)

    for i = 255:-1:2
        # New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
        # need inverse operator of y = exp(-x) -> x = -ln(y)
        x = -log(exp_section_area/x1 + feb[i+1])
        ke[i+1] = uint64(itrunc(x/x1*emantissa))
        web[i] = x/emantissa
        feb[i] = exp(-x)
        x1 = x
    end
    ke[2] = zero(Uint64)

    wi[:] = float64(wib)
    fi[:] = float64(fib)
    we[:] = float64(web)
    fe[:] = float64(feb)
    return nothing
end
randmtzig_fill_ziggurat_tables()
@test all(ki == Base.dSFMT.ki)
@test all(wi == Base.dSFMT.wi)
@test all(fi == Base.dSFMT.fi)
@test all(ke == Base.dSFMT.ke)
@test all(we == Base.dSFMT.we)
@test all(fe == Base.dSFMT.fe)

#same random numbers on for small ranges on all systems

seed = rand(Uint) #leave state nondeterministic as above
srand(seed)
r = int64(rand(int32(97:122)))
srand(seed)
@test r == rand(int64(97:122))

srand(seed)
r = uint64(rand(uint32(97:122)))
srand(seed)
@test r == rand(uint64(97:122))

@test all([div(0x000100000000,k)*k - 1 == Base.Random.RandIntGen(uint64(1:k)).u for k in 13 .+ int64(2).^(1:30)])
@test all([div(0x000100000000,k)*k - 1 == Base.Random.RandIntGen(int64(1:k)).u for k in 13 .+ int64(2).^(1:30)])

import Base.Random: uuid4, UUID

# UUID
a = uuid4()
@test a == UUID(string(a)) == UUID(utf16(string(a))) == UUID(utf32(string(a)))
@test_throws ArgumentError UUID("550e8400e29b-41d4-a716-446655440000")
@test_throws ArgumentError UUID("550e8400e29b-41d4-a716-44665544000098")
@test_throws ArgumentError UUID("z50e8400-e29b-41d4-a716-446655440000")

max_stars_repo_name: alanedelman/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
module NonlinearSolvers

#pde_pumi_interface_path = joinpath(Pkg.dir("PumiInterface"), "src")
#push!(LOAD_PATH, pde_pumi_interface_path)
using ArrayViews
import ArrayViews.view
using PDESolver
using PdePumiInterface  # needed to write vtk files
using ODLCommonTools
using Utils
import ODLCommonTools.sview
import MPI
using PETSc2  #TODO: see if this can be removed
using SummationByParts
using Jacobian
using LinearSolvers

import LinearSolvers: calcPC, applyPC, applyPCTranspose, calcLinearOperator,
                      applyLinearOperator, applyLinearOperatorTranspose,
                      needParallelData

include(joinpath(Pkg.dir("PDESolver"), "src/solver/debug.jl"))  # debug macro
include("rk4.jl")
include("lserk.jl")
export lserk54

include("jac_recalc.jl")
include("newton.jl")
export getNewtonPCandLO
include("preconditioning_types.jl")

include("crank_nicolson.jl")
include("preconditioning.jl")
include("globalization.jl")
include("predictor_corrector.jl")
include("p_homotopy.jl")
#include("newton_fd_old.jl")

#-----------------------------------------------------------
# NON WORKING AND EXPERIMENTAL!!!
# unsteady adjoint crank nicolson code:
include("crank_nicolson_uadj/crank_nicolson_uadj.jl")
include("crank_nicolson_uadj/crank_nicolson_uadj_jacandrhs.jl")
include("crank_nicolson_uadj/crank_nicolson_uadj_objective.jl")
include("crank_nicolson_uadj/crank_nicolson_uadj_checkpointstraight.jl")
#-----------------------------------------------------------

# predictor_corrector.jl
export predictorCorrectorHomotopy

# p_homotopy.jl
export pHomotopy

end


text:
# Import IPhreeqc before running this function
include("../src/IPhreeqc.jl")
pqc_input = """
SOLUTION 0
  pH   7   charge
  Na   0.2
SOLUTION 1-10
  pH   7   charge
  Na   0.1
SELECTED_OUTPUT
  -reset false
  -totals Na
END
"""
id = IPhreeqc.CreateIPhreeqc()
res = IPhreeqc.LoadDatabase(id, "phreeqc.dat")
println("load the database $res")
# IPhreeqc.RunFile(id, "ex2")
res = IPhreeqc.RunString(id, pqc_input)
println("Runnin the string $res")
n_comp = IPhreeqc.GetComponentCount(id)
println("number of components $n_comp")
IPhreeqc.DestroyIPhreeqc(id)


text:
"""
```
mutation!(Φ, Ψ, QQ, det_HH, inv_HH, φ_new, y_t, s_t, s_t1, ϵ_t, c, n_mh_steps;
    parallel = false, poolmodel = false)
```

Mutate particles by taking Metropolis-Hastings steps in the `ϵ_t` space. This
function modifies `s_t` and `ϵ_t` in place and returns `accept_rate`.
"""
function mutation!(Φ::Function, Ψ::Function, QQ::Matrix{Float64},
                   det_HH::Float64, inv_HH::Matrix{Float64}, φ_new::Float64,
                   y_t::Vector{Float64}, s_t::M, s_t1::M,
                   ϵ_t::M, c::Float64, n_mh_steps::Int;
                   parallel::Bool = false,
                   poolmodel::Bool = false) where M<:AbstractMatrix{Float64}
    # Sizes
    n_obs = size(y_t, 1)
    n_particles = size(ϵ_t, 2)

    # Initialize vector of acceptances
    accept_vec = Vector{Int}(undef, n_particles)

    # Used to generate new draws of ϵ
    dist_ϵ = MvNormal(c^2 * diag(QQ))

    # Used to calculate posteriors
    scaled_det_HH = det_HH/(φ_new^n_obs)
    scaled_inv_HH = inv_HH*φ_new

    # Take Metropolis-Hastings steps
    if parallel
        #sendto(workers(), s_t1 = s_t1)
        #sendto(workers(), s_t = s_t)
        #sendto(workers(), ϵ_t = ϵ_t)
        sendto(workers(), dist_ϵ = dist_ϵ)
        # sendto(workers(), y_t = y_t)
        # sendto(workers(), Φ = Φ) ## Should be sent to all workers before calling the function
        # sendto(workers(), Ψ = Ψ) ## Should be sent to all workers before calling the function
        sendto(workers(), scaled_inv_HH = scaled_inv_HH)

        mh_steps_closure(i::Int) = mh_steps(Φ, Ψ, dist_ϵ, y_t, s_t1[:,i], s_t[:,i], ϵ_t[:,i],
                         scaled_det_HH, scaled_inv_HH, n_mh_steps;
                         poolmodel = poolmodel)
        @everywhere mh_steps_closure(i::Int) = mh_steps(Φ, Ψ, dist_ϵ, y_t, s_t1[:,i], s_t[:,i], ϵ_t[:,i],
                         scaled_det_HH, scaled_inv_HH, n_mh_steps;
                         poolmodel = poolmodel)

        s_t2, ϵ_t2, accept_vec2 = @sync @distributed (vector_reduce) for i in 1:n_particles
            vector_reshape(mh_steps_closure(i)...)
        end
        s_t .= s_t2
        ϵ_t .= ϵ_t2
        accept_vec .= vec(accept_vec2)
    else
         for i in 1:n_particles
            s_t[:,i], ϵ_t[:,i], accept_vec[i] =
                mh_steps(Φ, Ψ, dist_ϵ, y_t, s_t1[:,i], s_t[:,i], ϵ_t[:,i],
                         scaled_det_HH, scaled_inv_HH, n_mh_steps;
                         poolmodel = poolmodel)
        end
    end

    # Calculate and return acceptance rate
    accept_rate = sum(accept_vec) / (n_mh_steps*n_particles)
    return accept_rate
end

"""
```
mh_steps(Φ, Ψ, dist_ϵ, y_t, s_t1, s_t, ϵ_t, scaled_det_HH,
    scaled_inv_HH, n_mh_steps; poolmodel = false)
```

Take `n_mh_steps` many steps in the `ϵ_t` space for a single particle. Returns
the new `s_t`, `ϵ_t`, and the number of acceptances `accept`.
"""
function mh_steps(Φ::Function, Ψ::Function, dist_ϵ::MvNormal, y_t::Vector{Float64},
                  s_t1::Vector{Float64}, s_t::Vector{Float64}, ϵ_t::Vector{Float64},
                  scaled_det_HH::Float64, scaled_inv_HH::Matrix{Float64}, n_mh_steps::Int;
                  poolmodel::Bool = false)
    accept = 0

    # Compute posterior at initial ϵ_t
    if poolmodel
        post_1 = Ψ(s_t) # PoolModel measurement eq is the likelihood
    else
        post_1 = fast_mvnormal_pdf(y_t - Ψ(s_t), scaled_det_HH, scaled_inv_HH)
    end
    post_2 = fast_mvnormal_pdf(ϵ_t)
    post   = post_1 * post_2

    for j = 1:n_mh_steps
        # Draw ϵ_new and s_new
        ϵ_new = ϵ_t + rand(dist_ϵ)
        s_new = Φ(s_t1, ϵ_new)

        # Calculate posterior
        post_new_1 = poolmodel ? Ψ(s_new) : fast_mvnormal_pdf(y_t - Ψ(s_new),
                                                              scaled_det_HH, scaled_inv_HH)
        post_new_2 = fast_mvnormal_pdf(ϵ_new)
        post_new   = post_new_1 * post_new_2

        # Calculate α, probability of accepting the new particle
        α = post_new / post
        if rand() < α
            s_t  = s_new
            ϵ_t  = ϵ_new
            post = post_new
            accept += 1
        end
    end
    return s_t, ϵ_t, accept
end

"""
```
update_c(c, accept_rate, target_rate)
```

Return the new proposal covariance matrix scaling `c`, adaptively chosen given
`accept_rate` to match `target_rate`.
"""
@inline function update_c(c::Float64, accept_rate::Float64, target_rate::Float64)
    c*(0.95 + 0.1*exp(20*(accept_rate - target_rate))/(1 + exp(20*(accept_rate - target_rate))))
end


text:
"""
Return default condition based on `PNTD`. Has meaning of true or always.
"""
function default_condition end
default_condition(::PNTD) where {PNTD <: PnmlType} = Condition(true)
default_condition(::Type{PNTD}) where {PNTD <: PnmlType} = Condition(true)
default_condition(::PNTD) where {PNTD <: AbstractContinuousCore} = Condition(true)
default_condition(::Type{PNTD}) where {PNTD <: AbstractContinuousCore} = Condition(true)
default_condition(pntd::PNTD) where {PNTD <: AbstractHLCore} = Condition(true)#default_term(pntd)) #!
default_condition(::Type{PNTD}) where {PNTD <: AbstractHLCore} = Condition(true)

"""
Label of a Transition that determines when the transition fires.

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Condition{T} <: HLAnnotation
    text::Maybe{String}
    term::T # structure
    com::ObjectCommon

end
Condition() = Condition(nothing, true, ObjectCommon())
Condition(value) = Condition(nothing, value, ObjectCommon())
Condition(text, value) = Condition(text, value, ObjectCommon())
Condition(text, value, oc::ObjectCommon) = Condition{typeof(value)}(text, value, oc)
Condition(::Nothing) =  error("Condition must have a `value`, ")
Condition(::Maybe{String}, ::Nothing) = error("Condition must have a `value`, ")
Condition(::Maybe{String}, ::Nothing, ::ObjectCommon) = error("Condition must have a `value`, ")

"""
Evaluate a [`Condition`](@ref) instance.
"""
(mark::Condition)() = _evaluate(mark.term)


text:
@deprecate front(x) first(x)
@deprecate back(x) last(x)
@deprecate top(x) first(x)
#@deprecate find_root find_root! # 2020-03-31 - deprecate in v0.18, or when Julia 1.5 is released.
export find_root
const find_root = find_root!

@deprecate deque(::Type{T}) where {T} Deque{T}()


text:
_package = ARGS[1]
trace_out = ARGS[2]
outfile = ARGS[3]

_package = Symbol(_package)

package = @eval (import $(_package); $(_package))
modules = Set(keys(Base.loaded_modules))

package.include("precomp-utils.jl")

open(outfile, "w") do io
    write(io, read(joinpath(@__DIR__, "precomp-utils.jl")))
    println(io)
    for p in collect(modules)
        import_str = "@import_pkg $(p.name) $(repr(p.uuid))"
        psym = Symbol(p.name)
        if !isdefined(package, psym)
            package.eval(:(const $(psym) = Base.require($(p))))
            println(io, import_str)
        else
            println("$psym already defined in package")
        end
    end
    println(io)
    for line in readlines(trace_out)
        isempty(line) && continue
        try
            include_string(package, line)
            println(io, "@" * line)
        catch e
            @warn "Could not precompile $(line)" exception=e
        end
    end
end


text:
###
# LLVMBuilder -- reliable LLVM builds all the time.
#
# --llvm-asserts: Build the Release+Asserts version
# --llvm-check: Build a RelWithDebInfo+Asserts version on x86-64-linux-musl
#               and run the testsuite. This will build for all targets.
###

using BinaryBuilder

# Collection of sources required to build LLVM
llvm_ver = "8.0.1"
sources = [
    "https://github.com/llvm/llvm-project/releases/download/llvmorg-$(llvm_ver)/llvm-$(llvm_ver).src.tar.xz" =>
    "44787a6d02f7140f145e2250d56c9f849334e11f9ae379827510ed72f12b75e7",
    "https://github.com/llvm/llvm-project/releases/download/llvmorg-$(llvm_ver)/cfe-$(llvm_ver).src.tar.xz" =>
    "70effd69f7a8ab249f66b0a68aba8b08af52aa2ab710dfb8a0fba102685b1646",
    "https://github.com/llvm/llvm-project/releases/download/llvmorg-$(llvm_ver)/compiler-rt-$(llvm_ver).src.tar.xz" =>
    "11828fb4823387d820c6715b25f6b2405e60837d12a7469e7a8882911c721837",
    "https://github.com/llvm/llvm-project/releases/download/llvmorg-$(llvm_ver)/polly-$(llvm_ver).src.tar.xz" =>
    "e8a1f7e8af238b32ce39ab5de1f3317a2e3f7d71a8b1b8bbacbd481ac76fd2d1",
    "https://github.com/llvm/llvm-project/releases/download/llvmorg-$(llvm_ver)/libunwind-$(llvm_ver).src.tar.xz" =>
    "1870161dda3172c63e632c1f60624564e1eb0f9233cfa8f040748ca5ff630f6e",
    "https://github.com/llvm/llvm-project/releases/download/llvmorg-$(llvm_ver)/lld-$(llvm_ver).src.tar.xz" =>
    "9fba1e94249bd7913e8a6c3aadcb308b76c8c3d83c5ce36c99c3f34d73873d88",

    # Include our LLVM patches
    "patches",
]

llvm_ver = VersionNumber(llvm_ver)

# Since we kind of do this LLVM setup twice, this is the shared setup start:
script_setup = raw"""
# We want to exit the program if errors occur.
set -o errexit

cd $WORKSPACE/srcdir/

# First, move our other projects into llvm/projects
for f in *.src; do
    # Don't symlink llvm itself into llvm/projects...
    if [[ ${f} == llvm-*.src ]]; then
        continue
    fi

    # clang lives in tools/clang and not projects/cfe
    if [[ ${f} == cfe-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/clang
    elif [[ ${f} == polly-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/polly
    elif [[ ${f} == lld-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/lld
    else
        mv $(pwd)/${f} $(echo llvm-*.src)/projects/${f%-*}
    fi
done

# Next, boogie on down to llvm town
cd llvm-*.src

# Update config.guess/config.sub stuff
# update_configure_scripts

# Apply all our patches
for f in $WORKSPACE/srcdir/llvm_patches/*.patch; do
    echo "Applying patch ${f}"
    patch -p1 < ${f}
done
"""

# The very first thing we need to do is to build llvm-tblgen for x86_64-linux-gnu
# This is because LLVM's cross-compile setup is kind of borked, so we just
# build the tools natively ourselves, directly.  :/
script = script_setup * raw"""
# Build llvm-tblgen, clang-tblgen, and llvm-config
mkdir build && cd build
CMAKE_FLAGS="-DLLVM_TARGETS_TO_BUILD:STRING=host -DLLVM_TEMPORARILY_ALLOW_OLD_TOOLCHAIN:BOOL=ON"
cmake .. ${CMAKE_FLAGS}
make -j${nproc} llvm-tblgen clang-tblgen llvm-config

# Copy the tblgens and llvm-config into our destination `bin` folder:
mkdir -p $prefix/bin
mv bin/llvm-tblgen $prefix/bin/
mv bin/clang-tblgen $prefix/bin/
mv bin/llvm-config $prefix/bin/
"""

# We'll do this build for x86_64-linux-musl only, as that's the arch we're building on
platforms = [
    Linux(:x86_64, :musl),
]

# We only care about llvm-tblgen and clang-tblgen
products(prefix) = [
    ExecutableProduct(prefix, "llvm-tblgen", :llvm_tblgen)
    ExecutableProduct(prefix, "clang-tblgen", :clang_tblgen)
    ExecutableProduct(prefix, "llvm-config", :llvm_config)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

llvm_ARGS = filter(s->startswith(s, "--llvm"), ARGS)
filter!(s->!startswith(s, "--llvm"), ARGS)

# Build the tarball, overriding ARGS so that the user doesn't shoot themselves in the foot,
# but only do this if we don't already have a Tblgen tarball available:
tblgen_tarball = joinpath("products", "tblgen.x86_64-linux-musl.tar.gz")
if !isfile(tblgen_tarball)
    tblgen_ARGS = ["x86_64-linux-musl"]
    if "--verbose" in ARGS
        push!(tblgen_ARGS, "--verbose")
    end
    if "--debug" in ARGS
        push!(tblgen_ARGS, "--debug")
    end
    product_hashes = build_tarballs(tblgen_ARGS, "tblgen", llvm_ver, sources, script, platforms, products, dependencies)

    # Extract path information to the built tblgen tarball and its hash
    tblgen_tarball, tblgen_hash = product_hashes["x86_64-linux-musl"]
    tblgen_tarball = joinpath("products", tblgen_tarball)
else
    @info("Using pre-built tblgen tarball at $(tblgen_tarball)")
    using SHA: sha256
    tblgen_hash = open(tblgen_tarball) do f
        bytes2hex(sha256(f))
    end
end

# Take that tarball, feed it into our next build as another "source".
push!(sources, tblgen_tarball => tblgen_hash)

# Next, we will Bash recipe for building across all platforms
script = script_setup * raw"""
# This value is really useful later
LLVM_DIR=$(pwd)

# Let's do the actual build within the `build` subdirectory
mkdir build && cd build

# Accumulate these flags outside CMAKE_FLAGS,
# they will be added at the end.
CMAKE_CPP_FLAGS=""
CMAKE_CXX_FLAGS=""
CMAKE_C_FLAGS=""

CMAKE_FLAGS=""
if [[ "${CHECK}" == "0" ]]; then
    # Start the cmake flags off with building for both our host arch, NVidia and AMD
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TARGETS_TO_BUILD:STRING=\"host\;NVPTX\;AMDGPU\""
fi

# Allow us to build with GCC 4.X for at least a little while longer
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TEMPORARILY_ALLOW_OLD_TOOLCHAIN:BOOL=ON"

# Also target Wasm because Javascript is the Platform Of The Future (TM)
# Actually, turn this off because we need it to work with The Julia Of Today (TM)
#CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_EXPERIMENTAL_TARGETS_TO_BUILD:STRING=\"WebAssembly\""

if [[ "${CHECK}" == "0" ]]; then
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_BUILD_TYPE=Release"
else
    # We want a more debugable version
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_BUILD_TYPE=RelWithDebInfo"
    # make sure that asserts are on as well
    ASSERTS=1
fi

if [[ "${ASSERTS}" == "1" ]]; then
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_ENABLE_ASSERTIONS:BOOL=ON"
fi

# We want a build with no bindings
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_BINDINGS_LIST=\"\" "

# Turn off ZLIB and XML2
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_ENABLE_ZLIB=OFF"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_ENABLE_LIBXML2=OFF"

# Disable useless things like docs, terminfo, etc....
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_INCLUDE_DOCS=Off"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_ENABLE_TERMINFO=Off"
CMAKE_FLAGS="${CMAKE_FLAGS} -DHAVE_HISTEDIT_H=Off"
CMAKE_FLAGS="${CMAKE_FLAGS} -DHAVE_LIBEDIT=Off"

# We want a shared library
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_BUILD_LLVM_DYLIB:BOOL=ON"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_LINK_LLVM_DYLIB:BOOL=ON"

# Install things into $prefix, and make sure it knows we're cross-compiling
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_INSTALL_PREFIX=${prefix}"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_CROSSCOMPILING=True"

# Julia expects the produced LLVM tools to be installed into tools and not bin
# We can't simply move bin to tools since on MingW64 it will also contain the shlib.
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TOOLS_INSTALL_DIR=${prefix}/tools"

# Also build and install utils, since we want FileCheck, and lit
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_UTILS_INSTALL_DIR=${prefix}/tools"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_INCLUDE_UTILS=True -DLLVM_INSTALL_UTILS=True"

# Include perf/oprofile/vtune markers
if [[ ${target} == *linux* ]]; then
    CMAKE_FLAGS="${CMAKE_FLAGS} -DUSE_PERF=1"
    CMAKE_FLAGS="${CMAKE_FLAGS} -DUSE_OPROFILE=1"
fi
if [[ ${target} == *linux* ]] || [[ ${target} == *mingw32* ]]; then
    CMAKE_FLAGS="${CMAKE_FLAGS} -DUSE_INTEL_JITEVENTS=1"
fi

# Tell LLVM where our pre-built tblgen tools are
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TABLEGEN=${WORKSPACE}/srcdir/bin/llvm-tblgen"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCLANG_TABLEGEN=${WORKSPACE}/srcdir/bin/clang-tblgen"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_CONFIG_PATH=${WORKSPACE}/srcdir/bin/llvm-config"

# Explicitly use our cmake toolchain file
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_TOOLCHAIN_FILE=/opt/${target}/${target}.toolchain"

# Manually set the host triplet, as otherwise on some platforms it tries to guess using
# `ld -v`, which is hilariously wrong.
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_HOST_TRIPLE=${target}"

# Tell LLVM which compiler target to use, because it loses track for some reason
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_C_COMPILER_TARGET=${target}"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_CXX_COMPILER_TARGET=${target}"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_ASM_COMPILER_TARGET=${target}"

# For now we focus on building llvm, clang, polly, and compiler-rt.
# We would like to build libc++, libc++abi and libunwind eventually
# but we currently don't due to issues on ppc and windows with
# libc++ and the LLVM libunwind conflicts with the one Julia is using.
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TOOL_LIBUNWIND_BUILD=OFF"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TOOL_LIBCXX_BUILD=OFF"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TOOL_LIBCXXABI_BUILD=OFF"

# Patch D50167 broke Polly :(
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_POLLY_BUILD=OFF"

if [[ "${target}" == *apple* ]]; then
    # On OSX, we need to override LLVM's looking around for our SDK
    CMAKE_FLAGS="${CMAKE_FLAGS} -DDARWIN_macosx_CACHED_SYSROOT:STRING=/opt/${target}/${target}/sys-root"

    # LLVM actually won't build against 10.8, so we bump ourselves up slightly to 10.9
    export MACOSX_DEPLOYMENT_TARGET=10.9
    export LDFLAGS=-mmacosx-version-min=10.9

    # We need to link against libc++ on OSX
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_ENABLE_LIBCXX=ON"
fi

if [[ "${target}" == *apple* ]] || [[ "${target}" == *freebsd* ]]; then
    # On clang-based platforms we need to override the check for ffs because it doesn't work with `clang`.
    export ac_cv_have_decl___builtin_ffs=yes

    # We don't use X-ray on BSD systems
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCOMPILER_RT_BUILD_XRAY=OFF"
fi

if [[ "${target}" == *mingw* ]]; then
    CMAKE_CPP_FLAGS="${CMAKE_CPP_FLAGS} -remap -D__USING_SJLJ_EXCEPTIONS__ -D__CRT__NO_INLINE"
    # Windows is case-insensitive and some dependencies take full advantage of that
    echo "BaseTsd.h basetsd.h" >> /opt/${target}/${target}/include/header.gcc

    # We don't Polly on windows
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_POLLY_BUILD=OFF"
fi

if [[ "${target}" == *musl* ]]; then
    # Taken from https://git.alpinelinux.org/cgit/aports/tree/main/compiler-rt/APKBUILD
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCOMPILER_RT_INCLUDE_TESTS=ON"
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCOMPILER_RT_BUILD_SANITIZERS=OFF"
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCOMPILER_RT_BUILD_XRAY=OFF"
fi

if [[ "${target}" == *freebsd* ]]; then
    # On FreeBSD, we must force even statically-linked code to have -fPIC
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_POSITION_INDEPENDENT_CODE=TRUE"
fi

# Build!
cmake .. ${CMAKE_FLAGS} -DCMAKE_C_FLAGS="${CMAKE_CPP_FLAGS} ${CMAKE_C_FLAGS}" -DCMAKE_CXX_FLAGS="${CMAKE_CPP_FLAGS} ${CMAKE_CXX_FLAGS}"
cmake -LA || true
make -j${nproc} VERBOSE=1

# Test
if [[ "${CHECK}" == "1" ]]; then
    make check -j${nproc}
fi

# Install!
make install -j${nproc} VERBOSE=1

# move clang products out of $prefix/bin to $prefix/tools
mv ${prefix}/bin/clang* ${prefix}/tools/
mv ${prefix}/bin/scan-* ${prefix}/tools/
mv ${prefix}/bin/c-index* ${prefix}/tools/
mv ${prefix}/bin/git-clang* ${prefix}/tools/
mv ${prefix}/bin/lld* ${prefix}/tools/

# Life is harsh on Windows and dynamic libraries are
# expected to live alongside the binaries. So we have
# to copy the *.dll from bin/ to tools/ as well...
if [[ "${target}" == *mingw* ]]; then
    cp ${prefix}/bin/*.dll ${prefix}/tools/
fi

# Work around llvm-config bug by creating versioned symlink to libLLVM
# https://github.com/JuliaLang/julia/pull/30033
if [[ "${target}" == *darwin* ]]; then
    LLVM_VER=$(basename $(echo ${prefix}/tools/clang-*.*))
    ln -s libLLVM.dylib ${prefix}/lib/libLLVM-${LLVM_VER##*-}.dylib
fi

# Lit is a python dependency and there is no proper install target
cp -r ../utils/lit ${prefix}/tools/

# Lots of tools don't respect `$DSYMUTIL` and so thus do not find
# our cleverly-named `llvm-dsymutil`.  We create a symlink to help
# Those poor fools along:
#ln -s llvm-dsymutil ${prefix}/tools/dsymutil
"""

if "--llvm-check" in llvm_ARGS
   # BB is using musl as a platform and we don't want to run glibc binaries on it.
   @info("Restricting build to `x86_64-linux-musl`")
   platforms = [
        BinaryProvider.Linux(:x86_64, :musl)
   ]
else
    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = [
        BinaryProvider.Linux(:i686, :glibc),
        BinaryProvider.Linux(:x86_64, :glibc),
        BinaryProvider.Linux(:x86_64, :musl),
        BinaryProvider.Linux(:aarch64, :glibc),
        BinaryProvider.Linux(:armv7l, :glibc),
        BinaryProvider.Linux(:powerpc64le, :glibc),
        BinaryProvider.MacOS(),
        BinaryProvider.Windows(:i686),
        BinaryProvider.Windows(:x86_64),
        BinaryProvider.FreeBSD(:x86_64),
    ]
    platforms = expand_gcc_versions(platforms)
end

# The products that we will ensure are always built
products(prefix) = [
    # libraries
    LibraryProduct(prefix, "libLLVM",  :libLLVM)
    LibraryProduct(prefix, "libLTO",   :libLTO)
    LibraryProduct(prefix, "libclang", :libclang)
    # tools
    ExecutableProduct(prefix, "llvm-tblgen", :llvm_tblgen)
    ExecutableProduct(prefix, "clang-tblgen", :clang_tblgen)
    ExecutableProduct(prefix, "llvm-config", :llvm_config)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
config = ""
name = "LLVM"
if "--llvm-asserts" in llvm_ARGS
    config *= "ASSERTS=1\n"
    name *= ".asserts"
    @warn("Removing PPC64LE from the platform list")
    filter!(p-> p.arch != :powerpc64le, platforms)
else
    config *= "ASSERTS=0\n"
end

if "--llvm-check" in llvm_ARGS
   config *= "CHECK=1\n"
   name *= ".check"
else
   config *= "CHECK=0\n"
end

build_tarballs(ARGS, name, llvm_ver, sources, config * script, platforms, products, dependencies)

if !("--llvm-keep-tblgen" in llvm_ARGS)
    # Remove tblgen tarball as it's no longer useful, and we don't want to upload them.
    rm(tblgen_tarball; force=true)
end


text:
κ(d, r) = exp(- d^2 / (2*r^2))

function symmetricconvolution_cuba(; g = g, κ = κ, tᵢ = tᵢ, tⱼ = tⱼ)

    offset = 50

    targetintegral(t) = g(tⱼ-t[2]) * g(tᵢ-t[1]) * κ(t[1]-t[2])

    trans(x) = -offset + (2*offset)*x

    auxiliary(t) = ((2*offset)^2)*targetintegral(trans.(t))

    resultquad = divonne( (t,out) -> out[1] = auxiliary(t), 2, 1, maxevals = 30_000_000, rtol=1e-6)

    display(resultquad)

    resultquad.integral[1]

end

function test_symmetricconvolution_cuba(tᵢ, tⱼ, r; mass=1e8, accretion=0.5, wavelength = 5000.0)

    # define transfer function

    g = PhysicalTransferFunction(mass=mass, accretion=accretion, wavelength=wavelength)

    symmetricconvolution_cuba(; g = g, κ = d -> κ(d, r), tᵢ = tᵢ, tⱼ = tⱼ)

end

#-------------------------------------------------------------------------------

function nonsymmetricconvolution(; gₖ = gₖ, gₗ = gₗ, κ = κ, tᵢ = tᵢ, tⱼ = tⱼ)

    offset = 50

    targetintegral(t) = gₖ(tᵢ-t[2]) * gₗ(tⱼ-t[1]) * κ(t[1]-t[2])

    trans(x) = -offset + (2*offset)*x

    auxiliary(t) = ((2*offset)^2)*targetintegral(trans.(t))

    resultquad = divonne( (t,out) -> out[1] = auxiliary(t), 2, 1, maxevals = 30_000_000, rtol=1e-6)

    display(resultquad)

    resultquad.integral[1]

end

function test_nonsymmetricconvolution(tᵢ, tⱼ, r; mass=1e8, accretion=0.5, wavelengths = [3000.0; 5000.0])

    # define transfer functions

    gₖ = PhysicalTransferFunction(mass=mass, accretion=accretion, wavelength=wavelengths[1])

    gₗ = PhysicalTransferFunction(mass=mass, accretion=accretion, wavelength=wavelengths[2])

    nonsymmetricconvolution(; gₖ = gₖ, gₗ = gₗ, κ = d -> κ(d, r) , tᵢ = tᵢ, tⱼ = tⱼ)

end


text:
const R = 8.314 #J/(mol*K)
export R
const Na = 6.022e23 #molecules/mol
export Na
const kB = R/Na #J/(molecule*K)
export kB
const F = 96485.33212 #C/mol


text:
#
module SEM

using Reexport
@reexport using DiffEqBase

using LinearAlgebra,SparseArrays
using Zygote,Plots
using UnPack
using StaticArrays

import FastGaussQuadrature
import Zygote,NNlib
#--------------------------------------#
Base.:*(op::Function, x::AbstractArray) = op(x)
#--------------------------------------#
export linspace
linspace(zi::Number,ze::Number,n::Integer) = Array(range(zi,stop=ze,length=n))
#--------------------------------------#
export iscallable
iscallable(op) = !isempty(methods(op))
#--------------------------------------#
sum(A::Array, n::Int) = Base.sum(A, dims=n)
sum(A) = Base.sum(A)
#cumprod(A::Array) = Base.cumprod(A, dims=1)
#cumprod(A::Array, d::Int) = Base.cumprod(A, dims=d)
#flipdim(A, d) = reverse(A, dims=d)
#--------------------------------------#
"""
 DiffEq ecosystem bindings
"""
#abstract type SEMPDEAlgorithm <: AbstractPDEAlgorithm end
#
#struct EllipticPDEAlgorithm <: SEMPDEAlgorithm end
#
#function DiffEqBase.__solve{}(
#                              prob::AbstractSEMPDEProblem
#                             )
#
#    build_solution(prob,alg,ts,timeseries,
#                   du = dures,
#                   dense = dense,
#                   timeseries_errors = timeseries_errors,
#                   retcode = :Success)
#end

#--------------------------------------#

# SEM building blocks
include("jac.jl")
include("interp.jl")
include("derivMat.jl")
include("semmesh.jl")
include("semq.jl")
include("ndgrid.jl")

include("ABu.jl")
include("geom.jl")
include("mask.jl")

include("time.jl")
include("mesh.jl")

include("pcg.jl")
include("adjoint.jl")

include("gatherScatter.jl")
include("grad.jl")
include("mass.jl")
include("lapl.jl")
include("hlmz.jl")
include("advect.jl")
include("diver.jl")

include("diffusion.jl")
include("convectionDiffusion.jl")
#include("stokes.jl")

include("plt.jl")

end # module
#--------------------------------------#


text:
include("production.jl")

"""A dotted production, also known as an *Item*."""
struct Item <: AbstractProduction
    left::Union{AbstractString, AbstractArray}
    right::Union{AbstractString, AbstractArray}
    pos::Int
    function Item(left, right, pos::Int = 1)
        if pos < 1 || pos > length(right)+1 throw(ArgumentError("The dot position is invalid.")) end
        p = Production(left, right)
        return Item(p, pos)
    end
    Item(prod::Production, pos::Int = 1) = new(prod.left, prod.right, pos)
end

### Functions ###

"""
    parseitem(input::AbstractString, iscontextfree::Bool = true)::Array{Item}
Returns an Array of [`Item`](@ref) obtained from the given string.
"""
parseitem(input::AbstractString, iscontextfree::Bool = true)::Array{Item} = parseproduction(input, iscontextfree) |> (x->map(y-> Item(y), x))

"""
    astype0(p::Item)::Item
Returns a new [`Item`](@ref) that is type 0
"""
astype0(i::Item)::Item = isa(i.left, AbstractArray) ? i : Item([i.left], i.right, i.pos)

"""
    afterdotsymbol(item::Item)::Union{AbstractString,Nothing}
Returns the symbol after the dot.
"""
afterdotsymbol(item::Item)::Union{AbstractString,Nothing} = item.pos <= length(item.right) ? item.right[item.pos] : nothing

"""
    advance(item::Item, x::AbstractString)::Item
Returns a new [`Item`](@ref) obtained advancing the dot past the given symbol.
"""
advance(item::Item, x::AbstractString)::Union{Item,Nothing} = if item.pos <= length(item.right) && item.right[item.pos] == x return Item(item.left, item.right, item.pos+1) else nothing end

### Operators ###

Base.:(==)(x::Item, y::Item) = (x.left, x.right, x.pos) == (y.left, y.right, y.pos)

Base.:<(x::Item, y::Item) = (x.left, x.right, x.pos) == (y.left, y.right, y.pos)

Base.:>(x::Item, y::Item) = !(x==y || x<y)

Base.hash(x::Item) = Base.hash((x.left, x.right, x.pos))

Base.show(io::IO, x::Item) = Base.show(io, string(x.left, "->", x.right[begin:x.pos-1],"•", x.right[x.pos:end]))

Base.iterate(x::Item, i...) = Base.iterate((x.left, x.right, x.pos), i...)

Base.isless(x::Item, y::Item) = Base.isless((x.left, x.right, x.pos), (y.left, y.right, y.pos))


text:
mutable struct VolumeMarket{T<:Real,N} <: OptimalBids.Market
    budget::T
    volume_fee::T
    allow_short_selling::Bool
    risk_free_rate::T
    volume_bids::Union{Vector{T},Nothing}
    clearing_prices::Union{Vector{T},Nothing}
end

function VolumeMarket(N::Int, budget::T, volume_fee::T, allow_short_selling::Bool, risk_free_rate::T) where {T<:Real}
    VolumeMarket{T,N}(budget, volume_fee, allow_short_selling, risk_free_rate, nothing, nothing)
end

function VolumeMarket(N::Int; budget::T=1.0, volume_fee::T=0.0, allow_short_selling::Bool=false, risk_free_rate::T=0.0) where {T<:Real}
    VolumeMarket(N, budget, volume_fee, allow_short_selling, risk_free_rate)
end

risk_free_rate(market::VolumeMarket) = market.risk_free_rate

market_budget(market::VolumeMarket) = market.budget
function set_market_budget(market::VolumeMarket{T,N}, val::T) where {T,N}
    market.budget = val
end
market_volume_fee(market::VolumeMarket) = market.volume_fee

function eltype(::VolumeMarket{T,N}) where {T,N}
    return T
end

function change_bids!(market::VolumeMarket{T, N}, new_bids::Vector{T}) where {T<:Real,N}
    @assert length(new_bids) == N
    total_volume = norm(new_bids, 1)
    if  total_volume > market.budget
        market.volume_bids = new_bids / total_volume
    else
        market.volume_bids = new_bids
    end

    return nothing
end

function clear_market!(market::VolumeMarket)
    @warn "Prices have not yet been set for this market.
    This function can only be dispatched passing cleared prices (`clearing_prices`)"
    throw(MethodError(clear_market!, (market)))
end

function clear_market!(market::VolumeMarket{T, N}, clearing_prices::Vector{T}) where {T<:Real,N}
    @assert length(clearing_prices) == N
    market.clearing_prices = clearing_prices
end

function calculate_profit(market::VolumeMarket)
    return (;
        cleared_volumes=market.volume_bids,
        clearing_prices=market.clearing_prices,
        profit=dot(market.volume_bids, market.clearing_prices),
        risk_free_profit=market.risk_free_rate * (market.budget - norm(market.volume_bids, 1))
    )
end

function total_profit(market::VolumeMarket)
    profits = calculate_profit(market)
    return sum(profits.profit) + profits.risk_free_profit
end

function length(::VolumeMarket{T, N}) where {T, N}
    return N
end

"""
    market_model(market::VolumeMarket{T,N}, optimizer_factory::Function) -> JuMP.Model, Vector{VariableRef}

Creates a JuMP model with appropriate PO variable and constraints:
    - Investment vector of variables `w` (portfolio weights if `budget = 1`).
    - Invested monney should be lower budget.

Returns the model and the reference to the vector of decision variables (length N).

Aditional arguments:
 - `optimizer_factory`: callable with zero arguments and return an empty `MathOptInterface.AbstractOptimizer``.
"""
function market_model(market::VolumeMarket{T,N}, optimizer_factory::Any;
    sense::MOI.OptimizationSense, model::JuMP.Model=Model(optimizer_factory)
) where {T,N}
    w = @variable(model, [1:N])
    sum_invested = @variable(model)
    if !market.allow_short_selling
        set_lower_bound.(w, 0.0)
        @constraint(model, sum_invested == sum(w))
    else
        @constraint(model, [sum_invested; w] in MOI.NormOneCone(length(w) + 1))
    end
    @constraint(model, sum_invested <= market_budget(market))
    obj = (1-sum_invested) * risk_free_rate(market) - sum_invested * market_volume_fee(market)
    if sense === MAX_SENSE
        @objective(model, sense, obj)
    else
        @objective(model, sense, - obj)
    end
    return model, w
end

"""
    change_bids!(market::VolumeMarket, model::JuMP.Model, decision_variables)

Solves optimization model and set new bids to optimal values of the `decision_variables`.
"""
function change_bids!(market::VolumeMarket, model::JuMP.Model, decision_variables)
    optimize!(model)
    status = termination_status(model)
    if status !== MOI.OPTIMAL
        @warn "Did not find an optimal solution: status=$status"
        change_bids!(market, zeros(length(decision_variables)))
    else
        new_bids = value.(decision_variables)
        change_bids!(market, new_bids)
    end

    return nothing
end

JuMP.owner_model(decision_variables::Vector{VariableRef}) = owner_model(first(decision_variables))
JuMP.owner_model(decision_variables::AffExpr) = owner_model(first(keys(decision_variables.terms)))
JuMP.owner_model(decision_variables::Vector{AffExpr}) = owner_model(first(decision_variables))

abstract type MarketHistory{T<:Real,N} end

function eltype(::MarketHistory{T,N}) where {T,N}
    return T
end

struct VolumeMarketHistory{T<:Real,N} <: MarketHistory{T,N}
    market::VolumeMarket{T,N}
    history_clearing_prices::Any
    history_risk_free_rates::Any
    timestamp::Any

    function VolumeMarketHistory{T,N}(market::VolumeMarket{T,N}, history_clearing_prices, history_risk_free_rates, timestamp) where {T,N}
        @assert eltype(values(history_clearing_prices)) === T
        @assert eltype(values(history_risk_free_rates)) === T
        @assert issetequal(keys(history_risk_free_rates), timestamp)
        @assert issetequal(keys(history_clearing_prices), timestamp)

        return new(market, history_clearing_prices, history_risk_free_rates, timestamp)
    end
end

function VolumeMarketHistory(market::VolumeMarket{T,N}, history_clearing_prices) where {T<:Real,N}
    rf = Dict(keys(history_clearing_prices) .=> market.risk_free_rate)
    return VolumeMarketHistory{T,N}(market, history_clearing_prices, rf, sort(unique(keys(history_clearing_prices))))
end

function VolumeMarketHistory(history_clearing_prices; kwards...)
    return VolumeMarketHistory(VolumeMarket(length(values(first(history_clearing_prices))); kwards...), history_clearing_prices)
end

number_assets(::VolumeMarketHistory{T,N}) where {T,N} = N

timestamp(hist::VolumeMarketHistory) = hist.timestamp
keys(hist::VolumeMarketHistory) = timestamp(hist)
function past_prices(hist::VolumeMarketHistory, t)
    timestamps = timestamp(hist)
    idx = findfirst(x -> x == t, timestamps)
    hist.history_clearing_prices[timestamps[1:idx-1]]
end
current_prices(hist::VolumeMarketHistory, t) = values(hist.history_clearing_prices[t])[1,:]
risk_free_rate(hist::VolumeMarketHistory, t) = first(values(hist.history_risk_free_rates[t]))
function size(hist::VolumeMarketHistory{T,N}) where {T,N}
    return (length(eachindex(hist.history_risk_free_rates)), N)
end

market_template(hist::VolumeMarketHistory) = hist.market

function market_template(hist, t)
    market = hist.market
    market.risk_free_rate = risk_free_rate(hist, t)
    return market
end

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
# function Mvz(vz_pml::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dz::Float64, dt::Float64, ext::Int64, iflag::Int64)
#     rho = modExpand(rho, ext, iflag)
#     (m,n) = size(rho)
#     a1 = 9/8;   a2 = -1/24;
#     c1 = a1/dz; c2 = a2/dz;
#     C3 = zeros(m*n)
#     denum = zeros(m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             if iz < m
#                a = (rho[iz+1,ix]+rho[iz,ix]) / (2*dt)
#             else
#                a = rho[iz,ix] / dt
#             end
#             b = vz_pml[iz,ix]/2
#             denum[(ix-1)*m+iz] = 1 / (a+b)
#             C3[(ix-1)*m+iz]= (a-b) / (a+b)
#         end
#     end
#     tmp = spzeros(m, m)
#     tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
#     for iz = 2: m-2
#         tmp[iz,iz  ] = -c1; tmp[iz, iz+1] = c1;
#         tmp[iz,iz-1] = -c2; tmp[iz, iz+2] = c2;
#     end
#     tmp[m-1,m-1] = -c1; tmp[m-1, m  ] = c1; tmp[m-1, m-2] = -c2;
#     tmp[m  ,m  ] = -c1; tmp[m  , m-1] =-c2;
#     MvzBvz = spdiagm(C3)
#     MvzBp  = spdiagm(denum) * kron(speye(n), tmp)
#     return MvzBvz, MvzBp
# end

function Mvz(vz_pml::SparseMatrixCSC{Float64,Int64}, dz::Float64, dt::Float64, ext::Int64, iflag::Int64)
    (m,n) = size(vz_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dz; c2 = a2/dz;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt + vz_pml[iz,ix]/2
            b = 1/dt - vz_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = 1 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(m, m)
    tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
    for iz = 2: m-2
        tmp[iz,iz  ] = -c1; tmp[iz, iz+1] = c1;
        tmp[iz,iz-1] = -c2; tmp[iz, iz+2] = c2;
    end
    tmp[m-1,m-1] = -c1; tmp[m-1, m  ] = c1; tmp[m-1, m-2] = -c2;
    tmp[m  ,m  ] = -c1; tmp[m  , m-1] =-c2;
    MvzBvz = spdiagm(C3)
    MvzBp  = spdiagm(denum) * kron(speye(n), tmp)
    return MvzBvz, MvzBp
end

function Mvz_back(vz_pml::SparseMatrixCSC{Float64,Int64}, dz::Float64, dt::Float64, ext::Int64, iflag::Int64)
    (m,n) = size(vz_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dz; c2 = a2/dz;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt - vz_pml[iz,ix]/2
            b = 1/dt + vz_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = 1 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(m, m)
    tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
    for iz = 2: m-2
        tmp[iz,iz  ] = -c1; tmp[iz, iz+1] = c1;
        tmp[iz,iz-1] = -c2; tmp[iz, iz+2] = c2;
    end
    tmp[m-1,m-1] = -c1; tmp[m-1, m  ] = c1; tmp[m-1, m-2] = -c2;
    tmp[m  ,m  ] = -c1; tmp[m  , m-1] =-c2;
    MvzBvz = spdiagm(C3)
    MvzBp  = -spdiagm(denum) * kron(speye(n), tmp)
    return MvzBvz, MvzBp
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
using InteractiveUtils
export mass_coeff, attenuation_coeff, mean_free_path, cross_section, lookup
export XAAMDI
export XCOM
export ESTAR

export EnergyLoss
export TotalAttenuation
export CoherentScattering
export IncoherentScattering
export PhotoelectricAbsorption
export PairProductionNuclearField
export PairProductionElectronField
export TotalAttenuationWithoutCoherentScattering

export Collision
export Radiative
export Total
export CSDA
export RadiationYield
export DensityEffect
export transmission

export Photon, Proton, Electron, Alpha

"""
    lookup(s::DataSource, E::E, m::Element, p::Process)

Must be implemented.
"""
abstract type DataSource end

abstract type Process end

abstract type Particle end
struct Photon <: Particle
    energy::typeof(1.0MeV)
end
struct Electron <: Particle
    energy::typeof(1.0MeV)
end
struct Proton <: Particle
    energy::typeof(1.0MeV)
end
struct Alpha <: Particle
    energy::typeof(1.0MeV)
end

# XAAMDI
struct EnergyLoss <: Process end
struct TotalAttenuation <: Process end

# XCOM
struct CoherentScattering <: Process end
struct IncoherentScattering <: Process end
struct PhotoelectricAbsorption <: Process end
struct PairProductionNuclearField <: Process end
struct PairProductionElectronField <: Process end
struct TotalAttenuationWithoutCoherentScattering <: Process end

# ESTAR
struct Collision      <: Process end
struct Radiative      <: Process end
"""
    Total

Total stopping power.
"""
struct Total          <: Process end # merge with TotalAttenuation?
struct CSDA           <: Process end
"""
    RadiationYield

Fraction of energy converted into Bremsstrahlung photons.
"""
struct RadiationYield <: Process end
struct DensityEffect  <: Process end

function lookup(mat, pt, proc::Process)
    s = default_data_source(mat, pt, proc)
    lookup(s, mat, pt, proc)
end

function mass_coeff(mat, pt, proc::Process;
                    datasource = default_data_source(mat, pt, proc))

    lookup(datasource, mat, pt, proc)
end

function attenuation_coeff(mat, args...; density=density(mat), kw...)
    q = mass_coeff(mat, args...;kw...) * density
    uconvert(inv(cm), q)
end

function cross_section(mat, args...; kw...)
    q = mass_coeff(mat, args...;kw...) * mat.atomic_mass
    uconvert(cm^2, q)
end

function stopping_power(mat, pt, proc::Process; datasource=default_data_source(mat,E,proc))
    q = lookup_stopping_power(datasource,mat,pt,proc)
    uconvert(MeV*cm^2/g, q)
end

mean_free_path(mat, pt, proc::Process;kw...) = uconvert(cm, 1 / attenuation_coeff(mat, pt, proc; kw...))

function transmission(mat, pt, proc::Process=EnergyLoss(); length, kw...)
    charlen = mean_free_path(mat, pt, proc; kw...)
    tau = uconvert(NoUnits, length/charlen)
    exp(-tau)
end


text:
using EMD, Test

run(ys) = (emdrice(ys), emd(ys))

@testset "Sin" begin
    let
        xs = range(0, 2π; length=200)
        ys = sin.(xs)
        expected, got = run(ys)
        @test length(expected[1]) == length(got[1])
        for (e, g) in zip(expected[1], got[1])
            @test all(g .≈ e)
        end
        @test all(got[2] .≈ expected[2])
    end
end

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore