text:
#=
Simulation for ESSE

Ignacio Quintero Mächler

t(-_-t)

January 12 2017
=#

"""
    simesse(λ       ::Array{Float64,1},
            μ       ::Array{Float64,1},
            Q       ::Array{Float64,2},
            x       ::Array{Float64,1},
            y       ::Array{Float64,1},
            esse_mod::String;
            nspp    ::Int64   = 20,
            δt      ::Float64 = 1e-3)

Simulate tree according to ESSE.
"""
function simesse(λ       ::Array{Float64,1},
                 μ       ::Array{Float64,1},
                 Q       ::Array{Float64,2},
                 β       ::Array{Float64},
                 x       ::Array{Float64,1},
                 y       ::Array{Float64};
                 esse_mod::String  = "speciation",
                 nspp    ::Int64   = 20,
                 δt      ::Float64 = 1e-3)

  # make simulation
  ed, el, st, simt, af, r = simedges(λ, μ, Q, β, x, y, nspp, δt, esse_mod)

  # organize in postorder
  ed     = numberedges(ed, nspp)
  ed, el = postorderedges(ed, el, nspp)

  tip_val = Dict(i => st[i] for i = 1:nspp)

  ##
  # organize new x and y for inference
  th = tree_height(el, ed, nspp)

  # is y multidimensional
  md = size(y,2) > 1

  # if simulation time is equal to simulated tree height
  if isapprox(th, simt, atol = 1e-9)

    idx = searchsortedlast(x, simt)

    nx = x[1:idx]

    push!(nx, th)

    if md
      ny = y[1:idx, :]
      af(th, r)
      ny = vcat(ny, reshape(r,(1,:)))
    else
      ny = y[1:idx]
      push!(ny, af(th))
    end

    ny = reverse(ny, dims = 1)

  # if simulation time is not equal to simulated tree height
  else
    idxs = searchsortedlast(x, simt)
    idxt = searchsortedlast(x, th)

    # new start
    nzero = simt - th

    idxz = searchsortedfirst(x, nzero)

    nx = x[idxz:idxs]
    ny = y[idxz:idxs, :]

    # push to time
    pushfirst!(nx, nzero)
    push!(nx, simt)

    # make x start at 0
    nx .-= nzero

    if md
      ny = y[idxz:idxs, :]
      af(nzero, r)
      ny = vcat(ny, reshape(r,(1,:)))
    else
      ny = y[idxz:idxs]
      push!(ny, af(nzero))
    end

    if md
      af(th, r)
      ny = vcat(ny, reshape(r,(1,:)))
    else
      push!(ny, af(th))
    end

    ny = reverse(ny, dims = 1)
  end

  tip_val = states_to_values(tip_val, lastindex(λ))

  return tip_val, ed, el, nx, ny
end

"""
    simedges(λ       ::Array{Float64,1},
             μ       ::Array{Float64,1},
             Q       ::Array{Float64,2},
             β       ::Array{Float64,1},
             af      ::Function,
             nspp    ::Int64,
             δt      ::Float64,
             esse_mod::String)

Simulate edges tree according to ESSE.
"""
function simedges(λ       ::Array{Float64,1},
                  μ       ::Array{Float64,1},
                  Q       ::Array{Float64,2},
                  β       ::Array{Float64,1},
                  x       ::Array{Float64,1},
                  y       ::Array{Float64},
                  nspp    ::Int64,
                  δt      ::Float64,
                  esse_mod::String)

  k = length(λ)

  # if multidimensional
  md = size(y, 2) > 1

  # check number of parameters and data
  ws, we, wq = id_mod(esse_mod, k, md, size(y, 2), λ, μ, β, Q)

  # total simulation time
  simt = 0.0

  # make approximate time function
  af = make_approxf(x, y)

  if md
    if wq
      r = Array{Float64,1}(undef, k*k - k)
    else
      r = Array{Float64,1}(undef, k)
    end
  end

  # get function estimates at `0.0`
  if md
    af(0.0, r)
  else
    r = af(0.0)
  end

  # preallocate vector of transitions probabilities after a transition has
  # been selected
  qpr = Array{Float64,1}(undef,k)
  spr = Array{Float64,1}(undef,k)

  # estimate starting q probabilities for initial states
  if wq
    # non-diagonal Q indices
    Qndi = setdiff(1:(k*k), 1:(k+1):(k*k))

    # Q column indices without diagonals
    Qci = Array{Int64,1}[]
    for i in Base.OneTo(k)
      push!(Qci, setdiff(1+k*(i-1): k+k*(i-1), 1:(k+1):(k*k)))
    end

    # β indices to match the ones for Q
    βi = UnitRange{Int64}[]
    for j in Base.OneTo(k)
      push!(βi, (j-1)k + (2-j):(j*k-j))
    end

    if md
      qsum = sum(l -> expf(Q[Qndi[l]], β[l], r[l]), 1:lastindex(Qndi))
      for j in Base.OneTo(k)
        qpr[j] =
          sum(i -> expf(Q[Qci[j][i]], β[βi[j][i]], r[βi[j][i]]), 1:(k-1))/qsum
      end
    else
      qsum = sum(l -> expf(Q[Qndi[l]], β[l], r), 1:lastindex(Qndi))
      for j in Base.OneTo(k)
        qpr[j] =
          sum(i -> expf(Q[Qci[j][i]], β[βi[j][i]], r), 1:(k-1))/qsum
      end
    end
  else
    qsum = sum(Q)
    sum!(qpr, transpose(Q))./qsum
  end

  # sample initial state
  si = prop_sample(spr, qpr, k)

  # n = 2 starting lineages (assume no stem branch)
  n = 2

  # edges alive
  ea = [1, 2]

  # state vector
  st = fill(si, 2)

  # edge array
  ed = zeros(Int64, 2nspp, 2)
  ed[ea,:] = [1 2;
              1 3]

  # edge lengths
  el = zeros(2nspp)

  # make λ, μ & Q probability functions
  λpr = make_λpr(λ, β, δt, md, ws)
  μpr = make_μpr(μ, β, δt, md, we)
  Qpr = make_Qpr(Q, β, δt, md, wq)

  # make function for estimating transition probabilites
  rqpr! = make_rqpr(Q, β, δt, md, wq)

  # start simulation
  while true

    # keep track of time
    simt += δt

    # estimate `z(simt)` or `z_i(simt)`
    if md
      af(simt, r)
    else
      r = af(simt)
    end

    # one time step
    for i in Base.OneTo(n)

      # update edge length
      el[ea[i]] += δt

      #=
          Q event?
      =#
      if rand() < Qpr(st[i], r)
        # get `out` transition probabilities
        rqpr!(st[i], qpr, r)
        # sample from the transition probabilities
        st[i] = prop_sample(spr, qpr, k)
      end

      #=
          μ event?
      =#
      if rand() < μpr(st[i], r)

        # update time in other extant lineages
        el[ea[(i+1):n]] .+= δt

        # node to remove
        nod = ed[ea[i],1]

        # top or bottom edge of node?
        top = ea[i] == findfirst(isequal(nod), ed[:,1])

        # remove edge (only preserving persisting species)
        ned = findfirst(isequal(0), ed)[1]

        if ned == 3
          printstyled("What would you do if an endangenered animal is eating an endangered plant? Sometimes nature is too harsh... \n",
            color=:red)
          error("tree went extinct... rerun simulation")

        end

        @views ed[ea[i]:(ned-1),:] = ed[(ea[i]+1):ned,:]

        # remove edge length
        @views el[ea[i]:(ned-1)] = el[(ea[i]+1):ned]

        # update alive lineages
        deleteat!(ea,i)
        ea[i:end] .-= 1

        # remove node and extra edge
        @views pr = findfirst(isequal(nod), ed[:,1])
        @views da = findfirst(isequal(nod), ed[:,2])

        if da == nothing
          da = 0
        end

        if da != 0
          ed[da,2] = ed[pr,2]
          el[da]  += el[pr]
        end

        # remove intervening edge
        @views ed[pr:(ned-2),:] = ed[(pr+1):(ned-1),:]

        # remove edge lengths of intervening node
        @views el[pr:(ned-2)] = el[(pr+1):(ned-1)]

        ## update alive lineages
        # is the remaining node terminal?
        if da == 0 || in(ed[da,2], ed[:,1])
          ea[i:end] .-= 1
        else
          ea[findfirst(isequal(pr),ea)] = da
          sort!(ea)
          if top
            ea[(i+1):end] .-= 1
          else
            ea[i:end]     .-= 1
          end
        end

        # update living species states
        deleteat!(st,i)

        # update n species
        n = lastindex(ea)

        # no more events at this time
        break
      end

      #=
          λ event?
      =#
      if rand() < λpr(st[i], r)

        ### add new edges
        # start node
        ed[ea[end] + 1,1] = ed[ea[i],2]
        ed[ea[end] + 2,1] = ed[ea[i],2]

        # end nodes
        ed[ea[end] + 1,2] = maximum(ed)+1
        ed[ea[end] + 2,2] = maximum(ed)+1

        # update time in other extant lineages
        el[ea[(i+1):n]] .+= δt

        # update living edges
        push!(ea, ea[end]+1, ea[end]+2)
        deleteat!(ea, i)

        # update living states
        push!(st,st[i],st[i])
        deleteat!(st,i)

        # update number of species alive
        n = lastindex(ea)

        # no more events at this time
        break
      end

    end

    n == nspp + 1 && break

  end

  ## remove last species and organize
  # identify daughter node and delete last two edges
  @views da = findfirst(isequal(ed[end,1]), ed[:,2])

  ed = ed[1:(2nspp-2),:]
  el = el[1:(2nspp-2)]

  # organize extant species
  pop!(ea)
  pop!(ea)
  push!(ea,da)

  sort!(ea)

  # organize extant states
  stf = st[end]
  pop!(st)
  pop!(st)

  insert!(st,findfirst(isequal(da), ea),stf)

  return ed, el, st, simt, af, r
end

"""
    id_mod(esse_mod::String, k::Int64,  md::Bool, nzt::Int64)

Check if number of parameters and `z(t)` functions
is consistent with specified model.
"""
function id_mod(esse_mod::String,
                k  ::Int64,
                md ::Bool,
                nzt::Int64,
                λ  ::Array{Float64,1},
                μ  ::Array{Float64,1},
                β  ::Array{Float64},
                Q  ::Array{Float64})

  allok = true

  # if speciation, extinction or transition model
  if occursin(r"^[s|S][A-za-z]*", esse_mod)

    mdd = "speciation"
    ws, we, wq = true, false, false

    if lastindex(λ) == lastindex(μ) == lastindex(β) == size(Q,1) == size(Q,2)
      printstyled("Number of parameters are consistent with specified model \n",
        color=:green)
    else
      allok = false
      printstyled("Parameter number for speciation ESSE with $k states:
                   λ = $(k) parameters
                   μ = $(k) parameters
                   Q = $k x $k parameter matrix
                   β = $(k) parameters",
        color=:red)
    end

  elseif occursin(r"^[e|E][A-za-z]*", esse_mod)

    if lastindex(λ) == lastindex(μ) == lastindex(β) == size(Q,1) == size(Q,2)
      printstyled("Number of parameters are consistent with specified model \n",
        color=:green)
    else
      allok = false
      printstyled("Parameter number for extinction ESSE with $k states:
                   λ = $(k) parameters
                   μ = $(k) parameters
                   Q = $k x $k parameter matrix
                   β = $(k) parameters",
        color=:red)
      error("Parameter number not consistent with specified model")
    end

    mdd = "extinction"
    ws, we, wq = false, true, false

  elseif occursin(r"^[t|T|r|R|q|Q][A-za-z]*", esse_mod)

    if lastindex(λ) == lastindex(μ) == size(Q,1) == size(Q,2) &&
       lastindex(β) == (k*k - k)
      printstyled("Number of parameters are consistent with specified model \n",
        color=:green)
    else
      allok = false
      printstyled("Parameter number for extinction ESSE with $k states:
                   λ = $(k) parameters
                   μ = $(k) parameters
                   Q = $k x $k parameter matrix
                   β = $(k*k - k) parameters",
        color=:red)
    end
    mdd = "transition"
    ws, we, wq = false, false, true
  else
    ws, we, wq = false, false, false
    warning("esse_mod does not match any ESSE model, running MuSSE \n")
  end

  if md
    if (ws || we) && nzt != k
      allok = false
      printstyled("For character-specific $mdd ESSE model for $k states:
               z(t) = $(k) functions",
        color=:red)
    end
    if wq && nzt != (k*k - k)
      allok = false
      printstyled("For character-specific $mdd ESSE model for $k states:
               z(t) = $(k*k - k) functions",
        color=:red)
    end
  end

  if allok
    printstyled("Simulating $mdd ESSE model with $k states \n",
      color=:green)
  else
    error("Parameter and z(t) function number not consistent with specified model")
  end

  return ws, we, wq
end

"""
    make_λpr(λ ::Array{Float64,1},
             β ::Array{Float64},
             δt::Float64,
             md::Bool,
             ws::Bool)

Make λ instantaneous probability function
"""
function make_λpr(λ ::Array{Float64,1},
                  β ::Array{Float64,1},
                  δt::Float64,
                  md::Bool,
                  ws::Bool)

  # if dependent on `z(t)` and multidimensional
  function f1(st ::Int64,
              aft::Array{Float64,1})
    return expf(λ[st], β[st], aft[st])*δt
  end

  # if dependent on `z(t)` and unidimensional
  function f2(st ::Int64,
              aft::Float64)
    return expf(λ[st], β[st], aft)*δt
  end

  # if **not** dependent on `z(t)` and multidimensional
  function f3(st ::Int64,
              aft::Array{Float64,1})
    return λ[st]*δt
  end

  # if **not** dependent on `z(t)` and unidimensional
  function f4(st ::Int64,
              aft::Float64)
    return λ[st]*δt
  end

  if ws
    if md
      return f1
    else
      return f2
    end
  else
    if md
      return f3
    else
      return f4
    end
  end
end

"""
    make_μpr(μ ::Array{Float64,1},
             β ::Array{Float64},
             δt::Float64,
             md::Bool,
             we::Bool)

Make `μ` instantaneous probability function
"""
function make_μpr(μ ::Array{Float64,1},
                  β ::Array{Float64,1},
                  δt::Float64,
                  md::Bool,
                  we::Bool)

  # if dependent on `z(t)` and multidimensional
  function f1(st ::Int64,
              aft::Array{Float64,1})
    return expf(μ[st], β[st], aft[st])*δt
  end

  # if dependent on `z(t)` and unidimensional
  function f2(st ::Int64,
              aft::Float64)
    return expf(μ[st], β[st], aft)*δt
  end

  # if **not** dependent on `z(t)` and multidimensional
  function f3(st ::Int64,
              aft::Array{Float64,1})
    return μ[st]*δt
  end

  # if **not** dependent on `z(t)` and unidimensional
  function f4(st ::Int64,
              aft::Float64)
    return μ[st]*δt
  end

  if we
    if md
      return f1
    else
      return f2
    end
  else
    if md
      return f3
    else
      return f4
    end
  end
end

"""
    make_Qpr(Q ::Array{Float64,2},
             β ::Array{Float64},
             δt::Float64,
             md::Bool,
             wq::Bool)

Make `Q` instantaneous probability function
"""
function make_Qpr(Q ::Array{Float64,2},
                  β ::Array{Float64,1},
                  δt::Float64,
                  md::Bool,
                  wq::Bool)
  k = size(Q,1)

  # Q row indices without diagonals
  Qri = Array{Int64,1}[]
  for i in Base.OneTo(k)
    push!(Qri, setdiff(i:k:(k*k), 1:(k+1):(k*k)))
  end

  # β indices to match the ones for Q
  βi = Array{Int64,1}[]
  for j in Base.OneTo(k)
    tv = Int64[]
    for i in Base.OneTo(k-1)
      push!(tv, Qri[j][i] - i)
    end
    push!(βi, tv)
  end

  ## Functions
  # if dependent on `z(t)` and multidimensional
  function f1(st ::Int64,
              aft::Array{Float64,1})
    rsum = 0.0
    for i in Base.OneTo(k-1)
      rsum += expf(Q[Qri[st][i]], β[βi[st][i]], aft[βi[st][i]])
    end
    return rsum*δt
  end

  # if dependent on `z(t)` and unidimensional
  function f2(st ::Int64,
              aft::Float64)
    rsum = 0.0
    for i in Base.OneTo(k-1)
      rsum += expf(Q[Qri[st][i]], β[βi[st][i]], aft)
    end
    return rsum*δt
  end

  # if **not** dependent on `z(t)` and multidimensional
  function f3(st ::Int64,
              aft::Array{Float64,1})
    return @views sum(Q[st,:])*δt
  end

  # if **not** dependent on `z(t)` and unidimensional
  function f4(st ::Int64,
              aft::Float64)
    return @views sum(Q[st,:])*δt
  end

  if wq
    if md
      return f1
    else
      return f2
    end
  else
    if md
      return f3
    else
      return f4
    end
  end
end

"""
    make_rqpr(Q ::Array{Float64,2},
              β ::Array{Float64},
              δt::Float64,
              md::Bool,
              wq::Bool)

Make cumulative probability for instantaneous change in
a given row in `Q` function.
"""
function make_rqpr(Q ::Array{Float64,2},
                   β ::Array{Float64},
                   δt::Float64,
                   md::Bool,
                   wq::Bool)
  k = size(Q,1)

  # Q row indices without diagonals
  Qri = Array{Int64,1}[]
  for i in Base.OneTo(k)
    push!(Qri, setdiff(i:k:(k*k), 1:(k+1):(k*k)))
  end

  # β indices to match the ones for Q
  βi = Array{Int64,1}[]
  for j in Base.OneTo(k)
    tv = Int64[]
    for i in Base.OneTo(k-1)
      push!(tv, Qri[j][i] - i)
    end
    push!(βi, tv)
  end

  ## Functions
  # if dependent on `z(t)` and multidimensional
  function f1(sti ::Int64,
              qpr::Array{Float64,1},
              aft::Array{Float64,1})
    l = 1
    for i in Base.OneTo(k)
      if i == sti
        qpr[i] = 0.0
      else
        qpr[i] = expf(Q[Qri[sti][l]], β[βi[sti][l]], aft[βi[sti][l]])
        l += 1
      end
    end

    return nothing
  end

  # if dependent on `z(t)` and unidimensional
  function f2(sti ::Int64,
              qpr::Array{Float64,1},
              aft::Float64)

    l = 1
    for i in Base.OneTo(k)
      if i == sti
        qpr[i] = 0.0
      else
        qpr[i] = expf(Q[Qri[sti][l]], β[βi[sti][l]], aft)
        l += 1
      end
    end

    return nothing
  end

  # if **not** dependent on `z(t)` and multidimensional
  function f3(sti ::Int64,
              qpr::Array{Float64,1},
              aft::Array{Float64,1})

    for i in Base.OneTo(k)
      qpr[i] = Q[sti,i]
    end

    return nothing
  end

  # if **not** dependent on `z(t)` and unidimensional
  function f4(sti ::Int64,
              qpr::Array{Float64,1},
              aft::Float64)
    for i in Base.OneTo(k)
      qpr[i] = Q[sti,i]
    end

    return nothing
  end

  if wq
    if md
      return f1
    else
      return f2
    end
  else
    if md
      return f3
    else
      return f4
    end
  end
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text: get_inner_vars(dynamic_device::PSY.DynamicInjection) = dynamic_device.ext[INNER_VARS]

ext: jl


text:
# Example with currying and piping

using Accounts, DataFrames

const FILE_ACCOUNTS = "./test_accounts.txt"

# declare Currying, see: https://riptutorial.com/julia-lang/example/20261/implementing-currying

curry(f, x) = (xs...) -> f(x, xs...)

add = curry(add_to_file, FILE_ACCOUNTS)

read = curry(read_from_file, FILE_ACCOUNTS)

# end Currying

email = create(EMAIL, "donald@duckcity.com")

donald = create("Donald Duck", [email])

# using Currying

[donald] |> add # How to use it

r = read() |> DataFrame

println(r)

# Remove data file

cmd = `rm $FILE_ACCOUNTS` # linux remove file statement

run(cmd) # remove file


text:
using ParallelAccelerator

#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)

@acc function opt_At_mul_B!(U, X, W)
	 #At_mul_B!(U, X, W)  # u <- w'x
	 U = X' * W
end

function main(m::Int, k::Int, n::Int)
    U = Array{Float64}(n, k)
    W = Array{Float64}(m, k)
    X = Array{Float64}(m, n)
    fill!(W, 3)
    fill!(X, 5)
    opt_At_mul_B!(U, X, W)
    println("done")

end

main(100, 100, 200)


text:
# Run from command line:
# julia train.jl
# julia train.jl --x 0.2
# julia train.jl --x 0.2 --noise 0.3
#
# Requires ArgParse:
#
# julia> import Pkg
# julia> Pkg.add("ArgParse")

using ArgParse

s = ArgParseSettings()
@add_arg_table! s begin
    "--x"
        arg_type = Float64
        default = 0.1
    "--noise"
        arg_type = Float64
        default = 0.1
end

args = parse_args(ARGS, s)
x = args["x"]
noise = args["noise"]

loss = (sin(5x) * (1 - tanh(x^2)) + randn() * noise)

println("loss:", loss)

max_stars_repo_name: wheatdog/guildai

max_issues_repo_name: wheatdog/guildai

max_forks_repo_name: wheatdog/guildai


text:
using FractalTools
using Debugger

function deneme(x,k)
    # k * x
    1
end

function deneme2(x,w)
    sin(w*x)
end

# An example of one dimensional ifs
A1 = reshape([1/2],1,1)
b1 = [0]
A2 = reshape([1/2],1,1)
b2 = [1/2]

w1 = Transformation(A1,b1)
w2 = Transformation(A2,b2)

w = [w1, w2]
ifs = IFS(w)
generator = ifs.generator
# generator = randalg_sequential_generator(ifs.ws, ifs.probs)

function indicator_measure_1D(x, x0, ϵ)
    return 2*ϵ * is_in_Ball(x,x0,ϵ)
end

elton_integral_1D(deneme, indicator_measure_1D, generator, (1,), (0.5, 0.25), 1e-11, 1e7)

elton_integral_1D(deneme2, indicator_measure_1D, generator, (2,), (0.6,0.2), 1e-11, 1e7)

# elton_integral(dipole_pocklington, is_in_Ball, generator,(0.9,0.005,2π), (0.9, 0.05))


text:
# This file was generated, do not modify it. # hide
p1 = :(model.evo_tree_classifier.η)
p2 = :(model.evo_tree_classifier.max_depth)

r1 = range(iterated_pipe, p1, lower=-2, upper=-0.5, scale=x->10^x)
r2 = range(iterated_pipe, p2, lower=2, upper=6)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# Tsitouras
struct Tsitouras45 <: AbstractAdaptiveExplicitSolver
end

function order(::Tsitouras45)::Int
    return 5 # ??? or just 5
end

function n_intermediates(::Tsitouras45)
    return 7 # to be checked
end

function n_temporaries(::Tsitouras45)
    return 1
end

function stepper!(ode_solver::Tsitouras45, ode_fcn, stepper_options, stepper_temp::StepperTemporaries)

    # renaming for better readability
    t          = stepper_temp.t_old;
	x          = stepper_temp.x_old;
	dt         = stepper_temp.dt;
    n_prob_dim = n_dim(stepper_temp)

    # views into the temporaries
    k1 = view(stepper_temp.der_x, 1:n_prob_dim, 1)
    k2 = view(stepper_temp.der_x, 1:n_prob_dim, 2)
    k3 = view(stepper_temp.der_x, 1:n_prob_dim, 3)
    k4 = view(stepper_temp.der_x, 1:n_prob_dim, 4)
    k5 = view(stepper_temp.der_x, 1:n_prob_dim, 5)
    k6 = view(stepper_temp.der_x, 1:n_prob_dim, 6)
    k7 = view(stepper_temp.der_x, 1:n_prob_dim, 7)

    x_tmp = view(stepper_temp.x_tmp, 1:n_prob_dim, 1)

    # x2 = view(stepper_temp.x_tmp, 1:n_prob_dim, 1)
    # x3 = view(stepper_temp.x_tmp, 1:n_prob_dim, 2)
    # x4 = view(stepper_temp.x_tmp, 1:n_prob_dim, 3)
    # x5 = view(stepper_temp.x_tmp, 1:n_prob_dim, 4)
    # x6 = view(stepper_temp.x_tmp, 1:n_prob_dim, 5)

    # older paper: this works
    # https://www.researchgate.net/publication/251743368_Runge-Kutta_Pairs_of_Orders_54_using_the_Minimal_Set_of_Simplifying_Assumptions
    c2 = 0.231572163526079
    c3 = 0.212252555252816
    c4 = 0.596693497318054
    c5 = 0.797009955708112
    c6 = 1

    b1 = 0.091937670648056
    b2 = 1.156529958312496
    b3 = -0.781330409541651
    b4 = 0.197624776163019
    b5 = 0.271639883438847
    b6 = 0.063598120979232

    bt1 = 0.092167469090589
    bt2 = 1.131750860603267
    bt3 = -0.759749304413104
    bt4 = 0.205573577541223
    bt5 = 0.264767065074229
    bt6 = 0.040490332103796
    bt7 = 1 / 40

    a32 = -0.059103796886580
    a42 = 4.560080615554683
    a52 = -2.443935658802774
    a62 = 9.516251378071800

    a43 = -4.006458683473722
    a53 = 2.631461258707441
    a63 = -8.467630087008555

    a54 = 0.524706566208284
    a64 = -0.987888827522473

    a65 = 0.867009765724064

    a21 = 0.231572163526079 #  c2
    a31 = 0.271356352139396 # c3 - (a32)
    a41 = 0.043071565237093434 # c4 - (a42 + a43)
    a51 = 0.08477778959516069 #  c5 - (a52 + a53 + a54)
    a61 = 0.0722577707351637 # c6 - (a62 + a63 + a64 + a65)

    #=
    c2 = 0.161
    c3 = 0.327
    c4 = 0.9
    c5 = 0.9800255409045097
    c6 = 1.0
    c7 = 1.0

    b1 = 0.09646076681806523
    b2 = 0.01
    b3 = 0.4798896504144996
    b4 = 1.379008574103742
    b5 = -3.290069515436081
    b6 = 2.324710524099774

    # original paper values
    # bt1 = 0.001780011052226
    # bt2 = 0.000816434459657
    # bt3 = -0.007880878010262
    # bt4 = 0.144711007173263
    # bt5 = -0.582357165452555
    # bt6 = 0.458082105929187
    # bt7 = 1 / 66

    # source: https://github.com/SciML/DiffEqDevTools.jl/blob/master/src/ode_tableaus.jl#L866-L920
    # the original paper values do not work
    # but I dont get where these values come from...
    bt1 = 0.09468075576583945
    bt2 = 0.009183565540343254
    bt3 = 0.4877705284247616
    bt4 = 1.2342975669304789
    bt5 = -2.707712349983525
    bt6 = 1.8666284181705870
    bt7 = 1 / 66

    # this still doesnt work
    # https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/tableaus/low_order_rk_tableaus.jl
    # bt1 = -0.001780011052226
    # bt2 = -0.000816434459657
    # bt3 = 0.007880878010262
    # bt4 = -0.144711007173263
    # bt5 = 0.582357165452555
    # bt6 = -0.458082105929187
    # bt7 = 1 / 66

    a32 = 0.3354806554923570
    a42 = -6.359448489975075
    a52 = -11.74888356406283
    a62 = -12.92096931784711

    a43 = 4.362295432869581
    a53 = 7.495539342889836
    a63 = 8.159367898576159

    a54 = -0.09249506636175525
    a64 = -0.07158497328140100

    a65 = -0.02826905039406838

    a21 = 0.161
    a31 = -0.008480655492356992
    a41 = 2.8971530571054944
    a51 = 5.32586482843926
    a61 = 5.86145544294642

    =#

    # intermediate steps
	t2 = t + c2 * dt
    # @. x2 = x + dt * a21 * k1
    @. x_tmp = x + dt * a21 * k1
    ode_fcn(k2, t2, x_tmp)

    t3 = t + c3 * dt
    # @. x3 = x + dt * (a31 * k1 + a32 * k2)
	@. x_tmp = x + dt * (a31 * k1 + a32 * k2)
    ode_fcn(k3, t3, x_tmp)

    t4 = t + c4 * dt
    # @. x4 = x + dt * (a41 * k1 + a42 * k2 + a43 * k3)
	@. x_tmp = x + dt * (a41 * k1 + a42 * k2 + a43 * k3)
    ode_fcn(k4, t4, x_tmp)

    t5 = t + c5 * dt
    # @. x5 = x + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4)
    @. x_tmp = x + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4)
    ode_fcn(k5, t5, x_tmp)

    t6 = t + c6 * dt
    # @. x6 = x + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)
    @. x_tmp = x + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)
    ode_fcn(k6, t6, x_tmp)

    # new step
	stepper_temp.t_new    = t + dt
	@. stepper_temp.x_new = x + dt * (b1 * k1 + b2 * k2 + b3 * k3 + b4 * k4 + b5 * k5 + b6 * k6)
	ode_fcn(k7, stepper_temp.t_new, stepper_temp.x_new)

    # error estimation
	@. stepper_temp.x_est = x + dt * (bt1 * k1 + bt2 * k2 + bt3 * k3 + bt4 * k4 + bt5 * k5 + bt6 * k6 + bt7 * k7)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
type SSVRGsd <: StepData

	t
	mtilde
	wtilde
	phat

	k::Function
	m::Int64
	getStep::Function
	stepString::AbstractString
	function SSVRGsd(k::Function, m::Int64,  numVars::Int64, stepString::AbstractString)
		new(0,0,zeros(numVars), zeros(numVars),k,m, SSVRGComputation, stepString)
	end
end

function SSVRG(k::Function, m::Int64,  numVars::Int64)
	stepString = "ssvrg.m=$m"
	SSVRGsd(k, m,numVars, stepString)
end

function productionOfSandU(s,u)
	k=0
    while true
		produce(s(k,0)+u(k,0))
		k=k+1
	end
end

function kComputation(t::Task,gnum,mtilde)
	thisSum =0
	for i=1:mtilde
		thisSum += consume(t)
	end
	max(thisSum-gnum,1)
end

function mytaskP(myarg)
	produce(myarg)
end

#
# function SSVRG(thisUsd::USD)
# 	thisProducer = @task  productionOfSandU(thisUsd.s,thisUsd.u)
# 	m = 100
#  	stepString = "ssvrgFromUm=100 $(thisEGRsd.stepString)"
# 	SSVRGsd((gnum,mtilde)->kComputation(thisProducer,gnum,mtilde), m,size(thisUsd.A)[1], stepString)
# end

function SSVRGComputation(x, k, gnum, sd::SSVRGsd, problem::Problem; outputLevel = 0)
	outputLevel>0 && println("starting step computation")
	if sd.t ==  sd.mtilde
		outputLevel>0 && println("in outer cycle, recomputing mtilde, wtilde, phat")
		phatsum=zeros(x)
		thisK = sd.k(gnum,sd.mtilde)
		for i = 1:thisK
			(f,sampleG) = (consume(problem.getNextSampleFunction))(x)
			outputLevel>0 && println("sampleG = $(sampleG)")
			phatsum += sampleG
		end
		gnum+=thisK
		sd.phat = phatsum/thisK
		sd.mtilde = rand(1:sd.m)
		sd.wtilde = copy(x)
		outputLevel>0 && println("sd.phat = $(sd.phat)")
		outputLevel>0 && println("sd.mtilde = $(sd.mtilde)")
		outputLevel>0 && println("sd.wtilde = $(sd.wtilde)")
		sd.t=0
	end
	kay = (consume(problem.getNextSampleFunction))
	(f,sampleG) = kay(x)
	(ft,sampleGt) = kay(sd.wtilde)
	gnum+=2
	g=sampleG-sampleGt+sd.phat
	sd.t+=1
	outputLevel>0 && println("finished step computation")
	(g, gnum)
end


text:
using TorsAttempt
using Test

@testset "TorsAttempt.jl" begin
    # Write your tests here.
end


text:
@testset "Prerequisites" begin
    @testset "Tools" begin include("Tools.jl") end
    @testset "TypeTraits" begin include("TypeTraits.jl") end
    @testset "Factories" begin include("Factories.jl") end
    @testset "CompositeStructures" begin include("CompositeStructures.jl") end
    @testset "SimpleTrees" begin include("SimpleTrees.jl") end
    @testset "NamedVectors" begin include("NamedVectors.jl") end
end


text:
# ![valid_cover](../logo.png)

# First valid image link is the cover


text:
using JuMP

m = Model()

macro persec(d,t)
	return quote
		for (k,v) in $d
			$d[k] = v / $t
		end
	end
end

@enum Item IronPlate IronOre Steel Items
@enum Assembler ElecSmelt Assemblers

iNames = Dict{Item, AbstractString}(IronPlate=>"Iron Plate", IronOre=>"Iron Ore", Steel=>"Steel")

n = Int64(Items)-1
f = Int64(Assemblers)-1

type Recipe
	outs::Dict{Item, Float64}
	ins::Dict{Item, Float64}
	function Recipe(time, outs, ins)
		@persec outs time
		@persec ins time
		new(outs, ins)
	end
end

function production(n::AbstractString, r::Recipe, nfacs)
	@printf "%s\n" n
	@printf "\tProduces\n"
	for (k,v) in r.outs
		@printf "\t\t%f %s per min\n" 60v * nfacs iNames[k]
	end
	@printf "\tConsumes\n"
	for (k,v) in r.ins
		@printf "\t\t%f %s per min\n" 60v * nfacs iNames[k]
	end
end

facSpeed = Dict{Assembler, Float64}(ElecSmelt=>2.0)

Steel_Time = 17.5/facSpeed[ElecSmelt]
Iron_Time = 3.5/facSpeed[ElecSmelt]
Ore_Time = 1 / 0.525

recips = Dict{Item, Recipe}()
recips[Steel] = Recipe(Steel_Time, Dict(Steel=>1.), Dict(IronPlate=>5.))
recips[IronPlate] = Recipe(Iron_Time, Dict(IronPlate=>1.), Dict(IronOre=>1.))
recips[IronOre] = Recipe(Ore_Time, Dict(IronOre=>1.), Dict(IronOre=>1.))

@variable(m, facSteel >= 1, Int)
@variable(m, facIron >= 1, Int)
@variable(m, facOre >= 1, Int)

@constraint(m, facSteel * recips[Steel].outs[Steel] >= 0)
@constraint(m, facIron * recips[IronPlate].outs[IronPlate] >= facSteel * recips[Steel].ins[IronPlate])
@constraint(m, facOre * recips[IronOre].outs[IronOre] >= facIron * recips[IronPlate].ins[IronOre])

@objective(m, Min, (facIron * recips[IronPlate].outs[IronPlate] - facSteel * recips[Steel].ins[IronPlate]))

solve(m)

production("Steel", recips[Steel], getvalue(facSteel))
production("Iron Plate", recips[IronPlate], getvalue(facIron))
production("Iron Ore", recips[IronOre], getvalue(facOre))

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
using WhereTraits
using Documenter

# The DocTestSetup does not work for our case, quite unsatisfying
# DocMeta.setdocmeta!(WhereTraits, :DocTestSetup, quote
#     using WhereTraits
#     using WhereTraits.BasicTraits
#     WhereTraits.BasicTraits.@overwrite_Base
# end, recursive=true)

makedocs(;
    modules=[WhereTraits, WhereTraits.BasicTraits],
    authors="Stephan Sahm <stephan.sahm@gmx.de> and contributors",
    repo="https://github.com/schlichtanders/WhereTraits.jl/blob/{commit}{path}#L{line}",
    sitename="WhereTraits.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://schlichtanders.github.io/WhereTraits.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => [
            "`@traits`" => "manual.md",
            "WhereTraits.BasicTraits" => "basictraits.md",
            "Combination with `isdef`" => "isdef.md",
        ],
        "Library" => "library.md",
    ],
    doctest = true,  # we haven't been able to get doctest running
)

deploydocs(;
    repo="github.com/schlichtanders/WhereTraits.jl",
)


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

using Markdown
using InteractiveUtils

# ╔═╡ 35dbb5f2-2ab1-11eb-3614-e7f46f8b16c2
using Images, Statistics, FFTW, Plots, DSP, ImageFIltering, PlutoUI

# ╔═╡ 2a33b6b0-2ab7-11eb-1987-69485426ad01
import Pkg; Pkg.add("Statistics")

# ╔═╡ 585724d0-2ab7-11eb-377d-17e95faba7b7

# ╔═╡ 54fb60ae-2ab1-11eb-04cf-c15e93f0893d

# ╔═╡ Cell order:
# ╠═35dbb5f2-2ab1-11eb-3614-e7f46f8b16c2
# ╠═2a33b6b0-2ab7-11eb-1987-69485426ad01
# ╠═585724d0-2ab7-11eb-377d-17e95faba7b7
# ╠═54fb60ae-2ab1-11eb-04cf-c15e93f0893d

max_stars_repo_name: zuckberj/comp-thinking-lectures

max_issues_repo_name: zuckberj/comp-thinking-lectures

max_forks_repo_name: zuckberj/comp-thinking-lectures


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# # Lagrangian particle tracking
#
# Materials and particles that tend to follow ocean currents can be analyzed in terms of trajectories. These are simply computed by integrating velocities over time within a [Lagrangian framework](https://en.wikipedia.org/wiki/Lagrangian_and_Eulerian_specification_of_the_flow_field).
#
# In `Julia` this is easily done e.g. using the [IndividualDisplacements.jl](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev/) package.

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## 1. Software, Grid, And Background Map

# + {"slideshow": {"slide_type": "-"}, "cell_style": "center"}
include("helper_functions.jl")
using IndividualDisplacements
pth=dirname(pathof(IndividualDisplacements))
include(joinpath(pth,"../examples/helper_functions.jl"));

# +
using OceanStateEstimation
OceanStateEstimation.get_ecco_velocity_if_needed() #download data if needed

function global_ocean_circulation(;k=1,ny=2)
  r_reset = 0.01 #fraction of the particles reset per month (0.05 for k<=10)

  #read grid and set up connections between subdomains
  γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
  Γ=GridLoad(γ;option="full")
  Γ=merge(Γ,NeighborTileIndices_cs(Γ))

  func=(u -> update_location_llc!(u,𝐷))
  Γ=merge(Γ,(; update_location! = func))

  #initialize u0,u1 etc
  𝑃,𝐷=set_up_FlowFields(k,Γ,ECCOclim_path);

  #add parameters for use in reset!
  tmp=(frac=r_reset, Γ=Γ)
  𝐷=merge(𝐷,tmp)

  return 𝑃,𝐷

end

𝑃,𝐷=global_ocean_circulation(k=20,ny=2); #grid etc

# +
lon=[i for i=-179.:2.0:179., j=-89.:2.0:89.]
lat=[j for i=-179.:2.0:179., j=-89.:2.0:89.]

(f,i,j,w)=InterpolationFactors(𝐷.Γ,vec(lon),vec(lat))
λ=(lon=lon,lat=lat,f=f,i=i,j=j,w=w);

OceanDepth=Interpolate(𝐷.Γ.Depth,λ.f,λ.i,λ.j,λ.w)
OceanDepth=reshape(OceanDepth,size(lon));

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## 2. Initialize Individuals

# + {"slideshow": {"slide_type": "-"}}
𝐷.🔄(𝑃,𝐷,0.0) #update velocity fields (here, to Dec and Jan bracketing t=0.0)

np=100
xy = init_global_randn(np,𝐷)
df=DataFrame(x=xy[1,:],y=xy[2,:],f=xy[3,:]) #initial positions

𝐼=Individuals(𝑃,df.x[1:np],df.y[1:np],df.f[1:np]) #Individuals data structure

fieldnames(typeof(𝐼))

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## 3. Compute Trajectories

# +
𝑇=(0.0,𝐼.𝑃.𝑇[2]) #first half of January
∫!(𝐼,𝑇) #mid-Dec to mid-Jan

for m=1:12
    𝐷.🔄(𝑃,𝐷,0.0) #update velocity fields
    ∫!(𝐼) #integrate forward by one more month
end

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## 4. Post-Processing

# + {"slideshow": {"slide_type": "subslide"}}
add_lonlat!(𝐼.🔴,𝐷.XC,𝐷.YC)
𝐼.🔴[end-3:end,:]

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## 5. Plot Trajectories
#
# In this example we simply map out individual positions (red to blue) over a map of ocean depth (log10).

# +
plt=contourf(vec(lon[:,1]),vec(lat[1,:]),permutedims(OceanDepth),clims=(-2500.,5000.),c = :grays, colorbar=false)
xli=extrema(lon)

🔴_by_t = groupby(𝐼.🔴, :t)
lo=deepcopy(🔴_by_t[1].lon); lo[findall(lo.<xli[1])]=lo[findall(lo.<xli[1])].+360
scatter!(lo,🔴_by_t[1].lat,markersize=1.5,c=:red,leg=:none,marker = (:circle, stroke(0)))
lo=deepcopy(🔴_by_t[end].lon); lo[findall(lo.<xli[1])]=lo[findall(lo.<xli[1])].+360
scatter!(lo,🔴_by_t[end].lat,markersize=1.5,c=:blue,leg=:none,marker = (:dot, stroke(0)))
# -

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
using Zygote
using Zygote: @adjoint, literal_getproperty
import Zygote: accum

import Distances: pairwise, colwise
const dtol = 1e-12 # threshold value for precise recalculation of distances

@nograd MersenneTwister, propertynames, broadcast_shape

function accum(D::Diagonal{T}, B::AbstractMatrix) where {T}
    A = Matrix{Union{T, Nothing}}(undef, size(D))
    A[diagind(A)] .= D.diag
    return accum(A, B)
end
accum(A::AbstractMatrix, D::Diagonal) = accum(D, A)
accum(A::Diagonal, B::Diagonal) = Diagonal(accum(diag(A), diag(B)))

@adjoint function colwise(s::SqEuclidean, x::AbstractMatrix, y::AbstractMatrix)
    return colwise(s, x, y), function (Δ::AbstractVector)
        x̄ = 2 .* Δ' .* (x .- y)
        return nothing, x̄, -x̄
    end
end

@adjoint function pairwise(s::SqEuclidean, X::AbstractMatrix; dims=2)
    D = pairwise(s, X; dims=dims)
    return D, function(Δ)
        d1, d2 = Diagonal(vec(sum(Δ; dims=1))), Diagonal(vec(sum(Δ; dims=2)))
        return (nothing, X * (2 .* (d1 .+ d2 .- Δ .- Δ')))
    end
end

@adjoint function colwise(s::Euclidean, x::AbstractMatrix, y::AbstractMatrix)
    d = colwise(s, x, y)
    return d, function (Δ::AbstractVector)
        x̄ = (Δ ./ d)' .* (x .- y)
        return nothing, x̄, -x̄
    end
end

@adjoint function pairwise(::Euclidean, X::AbstractMatrix, Y::AbstractMatrix; dims=2)
    @assert dims == 2
    D, back = Zygote.pullback((X, Y)->pairwise(SqEuclidean(dtol), X, Y; dims=2), X, Y)
    D .= sqrt.(D)
    return D, Δ -> (nothing, back(Δ ./ (2 .* D))...)
end

@adjoint function pairwise(::Euclidean, X::AbstractMatrix; dims=2)
    @assert dims == 2
    D, back = Zygote.pullback(X->pairwise(SqEuclidean(dtol), X; dims=2), X)
    D .= sqrt.(D)
    return D, function(Δ)
        Δ = Δ ./ (2 .* D)
        Δ[diagind(Δ)] .= 0
        return (nothing, first(back(Δ)))
    end
end

@adjoint function literal_getproperty(C::Cholesky, ::Val{:factors})
    error("@adjoint not implemented for :factors as is unsafe.")
    return literal_getproperty(C, Val(:factors)), function(Δ)
        error("@adjoint not implemented for :factors. (I couldn't make it work...)")
    end
end

import LinearAlgebra: HermOrSym, diag, Diagonal

diag(S::Symmetric{T, <:Diagonal{T}} where T) = S.data.diag
Zygote._symmetric_back(Δ::Diagonal) = Δ

# Diagonal matrices are always symmetric...
cholesky(A::HermOrSym{T, <:Diagonal{T}} where T) = cholesky(Diagonal(diag(A)))

#
# Some very specific broadcasting hacks while Zygote has crappy broadcasting.
#

import Base.Broadcast: broadcasted

@adjoint broadcasted(::typeof(-), x::AbstractArray) = .-x, Δ->(nothing, .-Δ,)

@adjoint function broadcasted(::typeof(exp), x::AbstractArray)
    y = exp.(x)
    return y, Δ->(nothing, Δ .* y)
end

@adjoint function \(D::Diagonal{<:Real}, B::AbstractVecOrMat{<:Real})
    Y = D \ B
    function ldiv_Diagonal_pullback(Ȳ::AbstractVecOrMat{<:Real})
        B̄ = D' \ Ȳ
        D̄ = Diagonal(-vec(sum(B̄ .* Y; dims=2)))
        return (D̄, B̄)
    end
    return Y, ldiv_Diagonal_pullback
end


text:
module Interpolate

using Geotherm.Geometry: Point2D, Point3D, Rectangle

export linear_interpolate, bilinear_interpolate

"""
    linear_interpolate(a, b)

Return a function of the linear interpolation between any 2 points `(x₁, f(x₁))` and `(x₂, g(x₂))`.
"""
function linear_interpolate(a::Point2D, b::Point2D)
    x₁, y₁, x₂, y₂ = a.x, a.y, b.x, b.y
    @assert x₁ != x₂ "The x-coordinates of the 2 arguments `a` and `b` cannot be equal!"
    return x -> ((x₂ - x) * y₁ + (x - x₁) * y₂) / (x₂ - x₁)
end

"""
    bilinear_interpolate(q₁₁, q₁₂, q₂₁, q₂₂)

If 4 points consisting of a rectangle and their z-coordinates are known, return a bilinear interpolation of the 4 points.
"""
function bilinear_interpolate(q₁₁::T, q₁₂::T, q₂₁::T, q₂₂::T) where {T<:Point3D}
    x₁, x₂, y₁, y₂ = q₁₁.x, q₂₁.x, q₁₁.y, q₂₂.y
    z₁₁, z₁₂, z₂₁, z₂₂ = map(q -> q.z, (q₁₁, q₁₂, q₂₁, q₂₂))
    rec = Rectangle(x₁, y₁, x₂, y₂)
    return function (x, y)
        if Point2D(x, y) ∉ rec
            error("point ($x, $y) is out of boundary $rec.")
        end
        v₁ = linear_interpolate(Point2D(x₁, z₁₁), Point2D(x₂, z₂₁))(x)
        v₂ = linear_interpolate(Point2D(x₁, z₁₂), Point2D(x₂, z₂₂))(x)
        return linear_interpolate(Point2D(y₁, v₁), Point2D(y₂, v₂))(y)
    end
end

end


text:
"""
    estimate_derivatives(object::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64, Vector{Float64}},x::Matrix{Float64},variable::Int64)

# Arguments

- **`object`** Output object from `BooST` function
- **`x`** Matrix of characteristics
- **`variable`** Column index of the desired variable

Returns the derivative of y with respect to `x[:,variable]`.
"""
function estimate_derivatives(object::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64, Vector{Float64}},x::Matrix{Float64},variable::Int64)::Vector{Float64}
    M::Int64 = object[5][5]
    N::Int64 = size(x,1)
    v::Float64 = object[5][1]
    rho::Vector{Float64} = object[7]
    vrho::Vector{Float64} = v*rho
    peaux::Matrix{Float64} = Matrix{Float64}(N,M)
    derivative::Vector{Float64} = Vector{Float64}(N)
    model::Dict{Float64,Array{Float64,2}} = object[1]
    for i in 1:M
        peaux[:,i] = gradient_st(model[i],x,variable)
    end
    derivative = peaux*vrho
    return derivative
end

"""
    BooSTMore(x::Matrix{Float64},y::Vector{Float64},object::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64, Vector{Float64}}; M::Int64=50, display::Bool=true)

# Arguments

- **`x`** Matrix of characteristics
- **`y`** Response variable vector
- **`object`** Output object from `BooST` function
- **`M`** Number of Boosting iterations to add
- **`display`** if `true`, prints iteration counter

Estimates more trees for a previously estimated BooST and returns the updated BooST.
"""
function BooSTMore(x::Matrix{Float64},y::Vector{Float64},object::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64, Vector{Float64}}; M::Int64=50, display::Bool=true)::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64,Vector{Float64}}

    (v::Float64,p::Float64,d_max::Int64,gamma::Vector{Float64},Mold::Int64,node_obs::Int64) = object[5]
    params::Tuple{Float64,Float64,Int64,Vector{Float64},Int64,Float64} = (v,p,d_max,gamma,M+Mold,node_obs)
    d_max = d_max - 1

    save_rho::Vector{Float64} = Vector{Float64}(Mold+M)
    save_rho[1:Mold] = object[7]
    ybar::Float64 = object[4]
    N::Int64 = length(y)
    phi::Vector{Float64} = predictBooST(object,x)
    brmse::Vector{Float64} = Vector{Float64}(Mold+M)
    save_tree::Dict{Float64,Array{Float64,2}} = object[1]
    brmse[1:Mold] = object[3]
    sizehint!(save_tree,Mold+M)
    for i in (Mold+1):(M+Mold)
        u::Vector{Float64} = y .-phi
        (tr::Matrix{Float64},fitted::Vector{Float64}) = grow_tree(x,u,p,d_max,gamma,node_obs)::Tuple{Array{Float64,2},Array{Float64,1}}
        rho::Float64 = sum([fitted[j]*u[j] for j in 1:N])/sum([fitted[j]*fitted[j] for j in 1:N])
        [phi[j] = phi[j] + v*rho*fitted[j] for j in 1:N]
        save_tree[i] = tr
        save_rho[i] = rho
        sqerror::Vector{Float64} = [(y[j]-phi[j])^2 for j in 1:N]
        rmse::Float64 = sqrt(mean(sqerror))
        brmse[i] = rmse
        if display == true
            println("iteration ",i, " rmse = ", rmse)
        end
    end
    nvar::Int64 = size(x,2)

    return save_tree, phi, brmse, ybar, params, nvar, save_rho

end

"""
    predictBooST(object::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64, Vector{Float64}},newx::Matrix{Float64})

# Arguments

- **`object`** Output object from `BooST` function
- **`newx`** Matrix of characteristics

Returns predicted values from a BooST model. If `newx = x` returns the fitted values.
"""
function predictBooST(object::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64, Vector{Float64}},newx::Matrix{Float64})::Vector{Float64}
    v::Float64 = object[5][1]
    y0::Float64 = object[4]
    rho::Vector{Float64} = object[7]
    model::Dict{Float64,Array{Float64,2}} = object[1]
    rhov::Vector{Float64} = rho*v
    N::Int64 = size(newx,1)
    M::Int64 = object[5][5]
    fittedaux::Matrix{Float64} = Matrix{Float64}(N,M)
    fitted::Vector{Float64} = Vector{Float64}(N)
    for i in 1:M
        fittedaux[:,i] = eval_tree(newx,model[i])
    end
    fitted = y0 + fittedaux*rhov
    return fitted
end

"""
    BooST(x::Matrix{Float64},y::Vector{Float64}; v::Float64 = 0.2, p::Float64 = 2/3, d_max::Int64 = 4, gamma::Vector{Float64} = collect(0.5:0.01:5), M::Int64 = 300, display::Bool = true, node_obs::Int64 = 1)

# Arguments

- **`x`** Matrix of characteristics
- **`y`** Response variable vector
- **`v`** Step size (shrinkage)
- **`p`** Proportion of variables tested in each new split (default 2/3)
- **`d_max`** number of splits
- **`gamma`** Transiction function intensity. Bigger numbers makes the transition less smoth. The default is a sequence of values (0.5:5) to be randomized in each new node. Multiple values may be supplied in a vector to increase the model randomness
- **`M`** Number of Boosting iterations to add
- **`display`** if `true`, prints iteration counter
- **`node_obs`** Equivalent to the minimum number of observations in a termina node for a discrete tree

Estimates the BooST and returns a tupple with the following objects below. The entire tupple is required to call the functions `BooSTMore`, `predictBooST` and `estimate_derivatives`.

# Outputs (Tupple order)

save_tree, phi, brmse, ybar, params, nvar, save_rho

- **`model`** Dict with information of each tree
- **`phi`** Fitted values
- **`brmse`** RMSE in each iteration
- **`ybar`** Mean of y
- **`params`** Parameters from call
- **`nvar`** Number of characteristics
- **`rho`** Vector with the parameter phi estimated in each iteration

"""
function BooST(x::Matrix{Float64},y::Vector{Float64}; v::Float64 = 0.2, p::Float64 = 2/3, d_max::Int64 = 4, gamma::Vector{Float64} = collect(0.5:0.01:5), M::Int64 = 300, display::Bool = true, node_obs::Int64 = 1)::Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64,Vector{Float64}}

    params::Tuple{Float64,Float64,Int64,Vector{Float64},Int64,Float64} = (v,p,d_max,gamma,M,node_obs)
    d_max::Int64 = d_max-1
    N::Int64 = length(y)
    ybar::Float64 = mean(y)
    phi::Vector{Float64} = ones(N)*ybar
    nvar::Int64 = size(x,2)

    brmse::Vector{Float64} = Vector{Float64}(M)
    save_rho::Vector{Float64} = Vector{Float64}(M)
    save_tree::Dict{Float64,Matrix{Float64}} = Dict{Float64,Matrix{Float64}}()
    sizehint!(save_tree,M)

    for i in 1:M
        u::Vector{Float64} = y .-phi
        (tr::Matrix{Float64},fitted::Vector{Float64}) = grow_tree(x,u,p,d_max,gamma,node_obs)::Tuple{Array{Float64,2},Array{Float64,1}}
        rho::Float64 = sum([fitted[j]*u[j] for j in 1:N])/sum([fitted[j]*fitted[j] for j in 1:N])
        [phi[j] = phi[j] + v*rho*fitted[j] for j in 1:N]
        save_tree[i] = tr
        save_rho[i] = rho
        sqerror::Vector{Float64} = [(y[j]-phi[j])^2 for j in 1:N]
        rmse::Float64 = sqrt(mean(sqerror))
        brmse[i] = rmse
        if display == true
            println("iteration ",i, " rmse = ", rmse)
        end

    end
    return save_tree, phi, brmse, ybar, params, nvar, save_rho
end

max_stars_repo_name: gabrielrvsc/BooSTjl

max_issues_repo_name: gabrielrvsc/BooSTjl

max_forks_repo_name: gabrielrvsc/BooSTjl


text:
##################
# map combinator #
##################

# used for type dispatch on the VectorTrace type (e.g. we will also have a UnfoldType)
struct MapType end

"""
    gen_fn = Map(kernel::GenerativeFunction)

Return a new generative function that applies the kernel independently for a vector of inputs.

The returned generative function has one argument with type `Vector{X}` for each argument of the input generative function with type `X`.
The length of each argument, which must be the same for each argument, determines the number of times the input generative function is called (N).
Each call to the input function is made under address namespace i for i=1..N.
The return value of the returned function has type `FunctionalCollections.PersistentVector{Y}` where `Y` is the type of the return value of the input function.
The map combinator is similar to the 'map' higher order function in functional programming, except that the map combinator returns a new generative function that must then be separately applied.
"""
struct Map{T,U} <: GenerativeFunction{PersistentVector{T},VectorTrace{MapType,T,U}}
    kernel::GenerativeFunction{T,U}
end

export Map

has_argument_grads(map_gf::Map) = has_argument_grads(map_gf.kernel)
accepts_output_grad(map_gf::Map) = accepts_output_grad(map_gf.kernel)

function (gen_fn::Map)(args...)
    (_, _, retval) = propose(gen_fn, args)
    retval
end

function get_args_for_key(args::Tuple, key::Int)
    map((arg) -> arg[key], args)
end

function get_prev_and_new_lengths(args::Tuple, prev_trace)
    new_length = length(args[1])
    prev_args = get_args(prev_trace)
    prev_length = length(prev_args[1])
    (new_length, prev_length)
end

include("assess.jl")
include("propose.jl")
include("simulate.jl")
include("generate.jl")
include("generic_update.jl")
include("update.jl")
include("regenerate.jl")
include("backprop.jl")

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


text:
# functions related to binomial distribution

import .RFunctions:
    binompdf,
    binomlogpdf,
    binomcdf,
    binomccdf,
    binomlogcdf,
    binomlogccdf,
    binominvcdf,
    binominvccdf,
    binominvlogcdf,
    binominvlogccdf


text:
using Test
using Pkg
using AutoViz
using Colors
using Random
using AutomotiveDrivingModels

@testset "Renderable" begin
    rw = gen_straight_roadway(3, 100.0)
    car = ArrowCar(0.0, 0.0, 0.0, id=1)
    car2 = ArrowCar(1.0, 1.0, 1.0, color=colorant"green", text="text")

    render([rw, car, "some text"])
    render([rw, car, car2, "some text"], camera=TargetFollowCamera(0, zoom=10.))
    render([rw, car, car2, TextOverlay(text=["overlay"], color=colorant"blue")])
    c = render([rw, car, car2], camera=SceneFollowCamera(zoom=10.))
end

@testset "write SVG, PDF, PNG" begin
    roadway = gen_stadium_roadway(4)
    c = @test_deprecated render(roadway)
    write("out.svg", c)
    @test isfile("out.svg")

    # write pdf
    camera = StaticCamera(position=(50,30), zoom=6.)
    c = render([roadway], camera=camera,
           surface=AutoViz.CairoPDFSurface(IOBuffer(), AutoViz.canvas_width(camera), AutoViz.canvas_height(camera)))
    write("out.pdf", c)

    # write png
    c = render([roadway], camera=camera,
           surface=AutoViz.CairoRGBSurface(AutoViz.canvas_width(camera), AutoViz.canvas_height(camera)))
    write("out.png", c)

    # try to write svg surface to pdf
    c = render([roadway], camera=camera)
    @test_throws ErrorException write("out.pdf", c)

    # try to write pdf surface to svg
    c = render([roadway], camera=camera,
           surface=AutoViz.CairoPDFSurface(IOBuffer(), AutoViz.canvas_width(camera), AutoViz.canvas_height(camera)))
    @test_throws ErrorException write("out.svg", c)

    # png should always work
    c = render([roadway], camera=camera)
    write("out2.png", c)
end

@testset "vehicle rendering" begin
    AutoViz.set_render_mode(:basic)
    @test AutoViz.rendermode == :basic

    roadway = gen_stadium_roadway(4)
    vehstate = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)

    def1 = VehicleDef()
    def2 = BicycleModel(def1)
    def3 = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)

    veh1 = Entity(vehstate, def1, 1)
    veh2 = Entity(vehstate, def2, 2)
    veh3 = Entity(vehstate, def3, 3)

    render([roadway, veh1, veh2, veh3])

    AutoViz.set_render_mode(:fancy)
    @test AutoViz.rendermode == :fancy

    render([roadway, veh1, veh2, veh3])
    scene = Frame([veh1])
    cam = TargetFollowCamera(1)
    update_camera!(cam, scene)
    render([Frame([veh1, veh2, veh3])], camera = cam)  # TODO: multiple dispatch not working on update_camera!
    render([Frame([veh1, veh2, veh3])], camera=StaticCamera(zoom=10.))
end

@testset "color theme" begin
    AutoViz.set_color_theme(OFFICETHEME)
    @test AutoViz.colortheme == OFFICETHEME

    roadway = gen_stadium_roadway(4)
    vehstate = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)

    def1 = VehicleDef()
    def2 = BicycleModel(def1)
    def3 = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)

    veh1 = Entity(vehstate, def1, 1)
    veh2 = Entity(vehstate, def2, 2)
    veh3 = Entity(vehstate, def3, 3)

    render([roadway, veh1, veh2, veh3])

    AutoViz.set_color_theme(LIGHTTHEME)
    @test AutoViz.colortheme == LIGHTTHEME

    render([roadway, veh1, veh2, veh3])

    AutoViz.set_color_theme(MONOKAY)
    @test AutoViz.colortheme == MONOKAY

    s = Frame([veh1])
    d = get_pastel_car_colors(s)
    @test length(d) == length(s)
end

@testset "doc examples" begin
    @testset "basics" begin
        include("../docs/lit/tutorials/basics.jl")
    end
    @testset "cameras" begin
        include("../docs/lit/tutorials/cameras.jl")
    end
    @testset "overlays" begin
        include("../docs/lit/tutorials/overlays.jl")
    end
end


text:
@info "test normal operator"
N = 512

Random.seed!(1234)
x = rand(N)
A = rand(N,N)
W = WeightingOp(rand(N))

y1 = adjoint(A)*W*A*x
y = normalOperator(A,W)*x

@test norm(y1 - y) / norm(y) ≈ 0 atol=0.01

y1 = adjoint(A)*A*x
y = normalOperator(A)*x

@test norm(y1 - y) / norm(y) ≈ 0 atol=0.01


text:
##
# Demo of using DualNumbers to differentiate with respect to
# changes in parameters
##

using DualNumbers, ApproxFun
# What is the derivative of the function (without differentiating the Chebyshev expansion)?
f=Fun(x->exp(dual(x,1)),[-1,1])
dualpart(f)
# check versus differentiate
norm(realpart(f)'-dualpart(f))

# What is the derivative of the first coefficient with respect to the left endpoint?
f=Fun(exp,Interval(dual(1.0,1),dual(2.0)))
dualpart(f.coefficients[1])
# check versus finite difference calculation:
h=0.00001;(Fun(exp,Interval(1.0+h,2.0)).coefficients[1]-Fun(exp,Interval(1.0,2.0)).coefficients[1])/h

# Or an ApproxFun calculation:
Fun(h->Fun(exp,Interval(1.0+h,2.0)).coefficients[1],[0.,.1])'(0.)

# What is the derivative of the first coefficient with respect to the exponential's constant?
f=Fun(x->exp(dual(x,x)),[-1,1])
dualpart(f.coefficients[1])
# check versus finite difference calculation:
h=0.00001;(Fun(x->exp((1+h)x),[-1,1]).coefficients[1]-Fun(exp,[-1,1]).coefficients[1])/h

max_stars_repo_path: examples/Dual Numbers.jl

max_issues_repo_path: examples/Dual Numbers.jl

max_forks_repo_path: examples/Dual Numbers.jl


text:
#Some useful functions

#Define the sigmoid function & its derivative
sigmoid(x)  = 1.0 / (1.0 + exp(-x));
dsigmoid(x)  = (1.0 / (1.0 + exp(-x))^2)*exp(-x);

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
function solve(r::Robot)
    mvs = []
    go_to_corner!(r, mvs) #Занять угол
    mark_corners!(r) #Покрасить углы
    return_back!(r, mvs) #Вернуться в начально положение
end

function MoveIfPossible(r::Robot, mvs::AbstractArray, side::HorizonSide)
    if !isborder(r, side)
        move!(r, side)
        push!(mvs, side)
    end
end

function go_to_corner!(r::Robot, mvs::AbstractArray)
    while isborder(r, Nord) == false || isborder(r, Ost) == false
            MoveIfPossible(r,mvs,Nord)
            MoveIfPossible(r,mvs,Ost)
    end
end

function moves!(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
        move!(r, side)
    end
end

function mark_corners!(r::Robot)
    for side in (Sud, West, Nord, Ost)
        moves!(r, side)
        putmarker!(r)
    end
end

function inverse(Side::HorizonSide)
    return HorizonSide((mod(Int(Side)+2, 4)))
end

function return_back!(r::Robot, mvs::AbstractArray)
    while length(mvs) > 0
        last_move = pop!(mvs)
        move!(r, inverse(last_move))
    end
end

main(r)


text:
using sem, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Arrow, DataFrames, LineSearches

miss20_dat = DataFrame(Arrow.Table("test/comparisons/dat_miss20_dat.arrow"))
miss20_par = DataFrame(Arrow.Table("test/comparisons/dat_miss20_par.arrow"))
miss20_par_mean = DataFrame(Arrow.Table("test/comparisons/dat_miss20_par_mean.arrow"))
miss20_mat = Matrix(miss20_dat)
miss20_start = DataFrame(Arrow.Table("test/comparisons/dat_miss20_start.arrow"))
miss20_start_mean = DataFrame(Arrow.Table("test/comparisons/dat_miss20_start_mean.arrow"))

three_path_dat = DataFrame(Arrow.Table("test/comparisons/three_path_dat.arrow"))

diff_fin =
    SemFiniteDiff(
        LBFGS(
            ;alphaguess = LineSearches.InitialHagerZhang(),
            linesearch = LineSearches.HagerZhang()
        ),
        Optim.Options(;f_tol = 1e-10, x_tol = 1.5e-8)
    )

## Model definition
@ModelingToolkit.variables x[1:31], mₓ, my[1:8], m[1:11]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

M = [mₓ, mₓ, mₓ, my[1:5]..., 3.0, my[7:8]..., 0.0, 0.0, 0.0]

M_free = [m[1:11]..., 0.0, 0.0, 0.0]

#S

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

# start_val = vcat(
#     vec(var(Matrix(three_path_dat), dims = 1))./2,
#     fill(0.05, 3),
#     fill(0.0, 6),
#     fill(1.0, 8),
#     fill(0, 3)
#     )

par_order = [collect(21:34); collect(15:20); 2;3; 5;6;7; collect(9:14);
    collect(35:45)]

par_order_mean = [collect(25:38); collect(15:20); 2;3; 5;6;7; collect(9:14); 21;
    collect(39:45)]

#= start_val_mean = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3),
    [mean(Matrix(three_path_dat), dims = 1)...][[1, 4:8..., 10:11...]]
    ) =#

start_val_mean = convert(Vector{Float64}, miss20_start_mean.est[par_order_mean])

#= start_val_mean_free = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3),
    [mean(Matrix(three_path_dat), dims = 1)...]
) =#

start_val_mean_free = convert(Vector{Float64}, miss20_start.est[par_order])

semobserved = SemObsMissing(miss20_mat)

imply = ImplySymbolic(A, S, F, [x..., m...], start_val_mean_free; M = M_free)
imply_mean = ImplySymbolic(A, S, F, [x..., mₓ, my[1:5]..., my[7:8]...], start_val_mean; M = M)

loss = Loss([SemFIML(semobserved, imply, [0.0], start_val_mean_free)])
loss_mean = Loss([SemFIML(semobserved, imply_mean, [0.0], start_val_mean)])

model_fin = Sem(semobserved, imply, loss, diff_fin)
model_fin_mean = Sem(semobserved, imply_mean, loss_mean, diff_fin)

solution_fin = sem_fit(model_fin)
solution_fin_mean = sem_fit(model_fin_mean)

pattern = model_fin.observed.patterns[55]
inv(model_fin.imply.imp_cov[pattern, pattern]) ≈ model_fin.loss.functions[1].inverses[55]
model_fin.loss.functions[1].logdets[55] ≈ logdet(model_fin.imply.imp_cov[pattern, pattern])

all(
    abs.(solution_fin.minimizer .- miss20_par.est[par_order]
        ) .< 0.05*abs.(miss20_par.est[par_order]))
all(
    abs.(solution_fin_mean.minimizer .- miss20_par_mean.est[par_order_mean]
        ) .< 0.05*abs.(miss20_par.est[par_order_mean]))

using FiniteDiff

grad_diff = FiniteDiff.finite_difference_gradient(model_fin, start_val_mean_free)

grad_fiml = sem.∇SemFIML(semobserved, imply, A, S, F, [x..., m...], start_val_mean_free; M = M_free)

diff_ana =
    SemAnalyticDiff(
        LBFGS(
            alphaguess = LineSearches.InitialHagerZhang(),
            linesearch = LineSearches.HagerZhang()
        ),
        Optim.Options(
            ;f_tol = 1e-10,
            x_tol = 1.5e-8),
            (grad_fiml,))

model_ana = Sem(semobserved, imply, loss, diff_ana)

grad = zeros(42)

model_ana(start_val_mean_free, grad)

all(
    abs.(grad .- grad_diff
        ) .< 0.001*abs.(grad_diff))
#solution_ana = sem_fit(model_ana)

solution_ana = sem_fit(model_ana)

all(
    abs.(solution_ana.minimizer .- miss20_par.est[par_order]
        ) .< 0.05*abs.(miss20_par.est[par_order]))

using ProfileView, Profile

ProfileView.@profview sem_fit(model_ana)

Profile.clear(); sem_fit(model_ana); ProfileView.view()


text:
module UnitfulAtomicHarmonic

import Unitful;
using Unitful: @unit, Dimension, Dimensions, dimension, Units, NoDims, NoUnits, uconvert, ustrip;
#import UnitfulAngles;

export ahunit, ahuconvert, ahustrip, initahu;

"""
    initahu(ω::Number, m::Number)
    initahu(ω::Unitful.Quantity, m::Unitful.Quantity)

Finishes the units initialisation by setting the harmonic frequency base unit
value to `ω` and the mass base unit to `m`. If neither `ω` of `m` have units,
`ω` is assumed to be in units of radians per second (`rad s⁻¹`) and `m` is
assumed to be in units of the atomic mass unit `u`.

Note that conversions to frequencies should use the `HzTn` unit defined in this
function if the frequency is to be expressed in terms of cycles per second
instead of radians per second.
"""
function initahu(ωh::Unitful.Quantity, mh::Unitful.Quantity)
  # Define base units.
  @eval @unit ħ_u "ħ"   ReducedPlanckConstant Unitful.ħ*Unitful.rad^-1                false;
  @eval @unit ω   "ω"   HarmonicFrequency     uconvert(Unitful.rad*Unitful.s^-1, $ωh) false;
  @eval @unit mₕ  "mₕ"  HarmonicMass          uconvert(Unitful.kg, $mh)               false;
  @eval @unit aₕ  "aₕ"  HarmonicLength        (1ħ_u / (1mₕ * 1ω))^(1/2)               false;

  # Define convenience units to aid conversions to non-angular quantities.
  @eval @unit τ   "τ"   Turn            2π*Unitful.rad  false;
  @eval @unit Tₕ  "Tₕ"  HarmonicPeriod  1τ/ω            false;

  # Hertz is usually expressed in terms of turns per second. The default Unitful
  # definition would make it in terms of radians per second instead.
  @eval @unit HzTn "Hz"  TurnHertz       1τ/Unitful.s true;

  # Simplify energy units.
  unit = :(ħ_u*ω);
  @eval ahunit(::typeof(dimension($unit))) = $unit

  # Only register symbols with Unitful when they are all defined.
  Unitful.register(UnitfulAtomicHarmonic);
end
# Assume any dimensionless numbers given are in angular frequency and atomic
# mass units.
initahu(ω::Number, m::Number) = initahu(ω*Unitful.rad*Unitful.s^-1, m*Unitful.u);

"""
  ahunit(x::Unitful.Quantity)
  ahunit(x::Unitful.Units)
  ahunit(x::Unitful.Dimensions)

Returns the appropriate atomic harmonic unit (a `Unitful.Units` object) for the
dimension of `x`.
"""
ahunit(x) = ahunit(dimension(x));

# Harmonic units for each re-defined dimension.
ahunit(x::Dimension{:Length})  = aₕ^x.power;
ahunit(x::Dimension{:Mass})    = mₕ^x.power;
ahunit(x::Dimension{:Time})    = Tₕ^x.power;

ahunit(::Dimension{D}) where D = throw(ArgumentError(string(
                                   "no atomic harmonic unit defined for the ",
                                   "dimension $D."
                                 )));

# Define ahunit() when multiple dimensions are passed.
@generated ahunit(::Dimensions{N}) where N = prod(ahunit, N);

# Define ahunit() for dimensionless quantities.
ahunit(::typeof(NoDims)) = NoUnits;

"""
    ahuconvert(x::Unitful.Quantity)

Convert a quantity to appropriate atomic harmonic units.
"""
ahuconvert(x) = uconvert(ahunit(x), x);

"""
    ahuconvert(u::Unitful.Units, x::Number)

Associate the passed units `u` with the passed number `x` which is in atomic
harmonic units.
"""
ahuconvert(u::Units, x::Number) = uconvert(u, x*ahunit(u));

"""
    ahustrip(x::Unitful.Quantity)

Return the value of `x` without its' units.
"""
ahustrip(x) = ustrip(ahuconvert(x));

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* File containg models of columns: distillation, stripping, absorbers
#* rectifier, ....
#*
#* The default nomenclature is:
#*		Type_Column_reboilertype_condensertyper
#*
#* where:
#*	Type = refluxed or reboiled or section
#*	Column = Stripping, Absorption, Rectifier, Distillation
#*	Reboiler type (if exists) = kettle or thermosyphon
#*	Condenser type (if exists) = with subccoling or without subcooling
#*
#*-----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: column.mso 511 2008-05-12 17:25:33Z paula $
#*---------------------------------------------------------------------
type Packed_Section_Column
	Packed_Section_Column()=begin
		new(
			Packed_Section_ColumnBasic(),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed stream",
				:PosX=>0,
				:PosY=>0.55
			)),
			[
				:(FeedStage.F= STAGES([1:_base_1.NumberOfStages]).Inlet.F),
				:(FeedStage.T = STAGES([1:_base_1.NumberOfStages]).Inlet.T),
				:(FeedStage.P = STAGES([1:_base_1.NumberOfStages]).Inlet.P),
				:(FeedStage.z = STAGES([1:_base_1.NumberOfStages]).Inlet.z),
				:(FeedStage.v = STAGES([1:_base_1.NumberOfStages]).Inlet.v),
				:(FeedStage.h = STAGES([1:_base_1.NumberOfStages]).Inlet.h),
				:(0*"mol/h"= STAGES([1:_base_1.NumberOfStages]).Inlet.F),
				:(300*"K" = STAGES([1:_base_1.NumberOfStages]).Inlet.T),
				:(1*"atm" = STAGES([1:_base_1.NumberOfStages]).Inlet.P),
				:(0.1 = STAGES([1:_base_1.NumberOfStages]).Inlet.z),
				:(0 = STAGES([1:_base_1.NumberOfStages]).Inlet.v),
				:(0*"J/mol" = STAGES([1:_base_1.NumberOfStages]).Inlet.h),
			],
			[
				"FeedTrayTop Inlet Flow","FeedTrayTop Inlet Temperature","FeedTrayTop Inlet Pressure","FeedTrayTop Inlet Composition","FeedTrayTop Inlet Vapour Fraction","FeedTrayTop Inlet Enthalpy","Inlet Tray - Flow Sealed","Inlet Tray -  Temperature","Inlet Tray -  Pressure","Inlet Tray -  Composition","Inlet Tray -  Vapour Fraction","Inlet Tray -  Enthalpy",
			],
			[:FeedStage,]
		)
	end
	_base_1::Packed_Section_ColumnBasic
	FeedStage::stream
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Packed_Section_Column
function set(in::Packed_Section_Column)
	_base_1.NumberOfFeeds = 1

end
function setEquationFlow(in::Packed_Section_Column)
	#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	# Equating Feed Tray Variables to Trays Variables
	#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	if i == FeedStageLocation(1)
		addEquation(1)
		addEquation(2)
		addEquation(3)
		addEquation(4)
		addEquation(5)
		addEquation(6)
	else
		addEquation(7)
		addEquation(8)
		addEquation(9)
		addEquation(10)
		addEquation(11)
		addEquation(12)
	end
end
function atributes(in::Packed_Section_Column,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/PackedSectionColumn"
	fields[:Brief]="Model of a packed column section."
	fields[:Info]="== Model of a packed column section containing ==
* NStages theoretical stages.

== Specify ==
* the feed stream of each tray (Inlet);
* the InletLiquid stream of the top tray;
* the InletVapour stream of the bottom tray;
* the total pressure drop (dP) of the section.

== Initial Conditions ==
* the stages temperature (OutletLiquid.T);
* the stages liquid holdup;
* (NoComps - 1) OutletLiquid (OR OutletVapour) compositions for each tray.
"
	drive!(fields,_)
	return fields
end
Packed_Section_Column(_::Dict{Symbol,Any})=begin
	newModel=Packed_Section_Column()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
@everywhere function create_model(scenario_id::PH.ScenarioID;
                                  opt=()->Ipopt.Optimizer(),
                                  opt_args=(print_level=0,)
                                  )

    model = JuMP.Model(opt)
    for (key,value) in pairs(opt_args)
        JuMP.set_optimizer_attribute(model, string(key), value)
    end

    scid = PH.value(scenario_id)

    c = [1.0, 10.0, 0.01]
    d = 7.0
    a = 16.0

    α = 1.0
    β = 1.0
    γ = 1.0
    δ = 1.0
    ϵ = 1.0

    s1 = 8.0
    s2 = 4.0
    s11 = 9.0
    s12 = 16.0
    s21 = 5.0
    s22 = 18.0

    stage1 = JuMP.@variable(model, x[1:3] >= 0.0)
    JuMP.@constraint(model, x[3] <= 1.0)
    obj = zero(JuMP.GenericQuadExpr{Float64,JuMP.VariableRef})
    JuMP.add_to_expression!(obj, sum(c.*x))

    # Second stage
    stage2 = Vector{JuMP.VariableRef}()
    if scid < 2
        vref = JuMP.@variable(model, y >= 0.0)
        JuMP.@constraint(model, α*sum(x) + β*y >= s1)
        JuMP.add_to_expression!(obj, d*y)
    else
        vref = JuMP.@variable(model, y >= 0.0)
        JuMP.@constraint(model, α*sum(x) + β*y >= s2)
        JuMP.add_to_expression!(obj, d*y)
    end
    push!(stage2, vref)

    # Third stage
    stage3 = Vector{JuMP.VariableRef}()
    if scid == 0
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s11)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))

    elseif scid == 1
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s12)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))

    elseif scid == 2
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s21)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))

    else
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s22)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))
    end
    append!(stage3, vref)

    JuMP.@objective(model, Min, obj)

    vdict = Dict{PH.StageID, Vector{JuMP.VariableRef}}(PH.stid(1) => stage1,
                                                       PH.stid(2) => stage2,
                                                       PH.stid(3) => stage3,
                                                       )

    return JuMPSubproblem(model, scenario_id, vdict)
end

function build_scen_tree()

    probs = [0.5*0.75, 0.5*0.25, 0.5*0.75, 0.5*0.25]

    tree = PH.ScenarioTree()

    for k in 1:2
        node2 = PH.add_node(tree, tree.root)
        for l in 1:2
            PH.add_leaf(tree, node2, probs[(k-1)*2 + l])
        end
    end
    return tree
end

function timeout_wait(t::Task, limit::Real=10, interval::Real=1)
    count = 0
    success = false
    while true
        sleep(interval)
        count += 1
        if istaskdone(t)
            success = true
            break
        elseif count >= limit
            break
        end
    end
    return success
end

function invert_map(my_map::Dict{A,B})::Dict{B,A} where {A,B}
    inv_map = Dict{B,A}()
    for (k,v) in pairs(my_map)
        inv_map[v] = k
    end
    @assert length(inv_map) == length(my_map)
    return inv_map
end

@everywhere function two_stage_model(scenario_id::PH.ScenarioID)

    model = JuMP.Model(()->Ipopt.Optimizer())
    JuMP.set_optimizer_attribute(model, "print_level", 0)
    JuMP.set_optimizer_attribute(model, "tol", 1e-12)
    JuMP.set_optimizer_attribute(model, "acceptable_tol", 1e-12)

    scen = PH.value(scenario_id)

    ref = JuMP.@variable(model, 0.0 <= x <= 1000.0)
    stage1 = [ref]
    push!(stage1, JuMP.@variable(model, 0.0 <= u <= 1.0))

    ref = JuMP.@variable(model, y >= 0.0)
    stage2 = [ref]

    val = scen == 0 ? 11.0 : 4.0

    JuMP.@constraint(model, x + y + u == val)

    c_y = scen == 0 ? 1.5 : 2.0
    JuMP.@objective(model, Min, 1.0*x + c_y * y)

    return PH.JuMPSubproblem(model, scenario_id, Dict(PH.stid(1) => stage1,
                                                      PH.stid(2) => stage2)
                             )
end

if isdefined(Main, :Xpress)
    function two_stage_int(scenario_id::PH.ScenarioID)
        model = JuMP.Model(()->Xpress.Optimizer())
        JuMP.set_optimizer_attribute(model, "OUTPUTLOG", 0)

        scen = PH.value(scenario_id)

        ref = JuMP.@variable(model, 4 <= k <= 10, Int)
        stage1 = [ref]

        ref = JuMP.@variable(model, y)
        stage2 = [ref]

        val = 2.0 * scen
        JuMP.@constraint(model, k + y == val)

        c_y = (scen % 2 == 0 ? 1.0 : 5.0)
        JuMP.@objective(model, Min, c_y * y^2 + 2.0*k)

        return PH.JuMPSubproblem(model, scenario_id, Dict(PH.stid(1) => stage1,
                                                          PH.stid(2) => stage2)
                                 )
    end
end

function build_var_map(n::Int,
                       m1::Int,
                       m2::Int)
    vmap = Dict{PH.ScenarioID, Dict{PH.VariableID, PH.VariableInfo}}()
    vval = Dict{PH.ScenarioID, Dict{PH.VariableID, Float64}}()
    for k in 1:n
        scid = PH.scid(k-1)

        vars = Dict{PH.VariableID, PH.VariableInfo}()
        vals = Dict{PH.VariableID, Float64}()

        # Stage 1 variables
        for j in 1:m1
            vid = PH.VariableID(scid, PH.stid(1), PH.index(j + k - 1))
            vars[vid] = PH.VariableInfo("a$j",false) # convert(Float64, j)
            vals[vid] = convert(Float64, j)
        end

        # Stage 2 variables
        for j in 1:m2
            vid = PH.VariableID(scid, PH.stid(2), PH.index(j-1))
            vars[vid] = PH.VariableInfo("b$j", false) # convert(Float64, k*j + n)
            vals[vid] = convert(Float64, k*j + n)
        end

        vmap[scid] = vars
        vval[scid] = vals
    end

    return (vmap, vval)
end

function fake_phdata(nv1::Int, nv2::Int; add_leaves=true)
    nscen = 2
    st = two_stage_tree(nscen)
    (var_map, var_val) = build_var_map(nscen, nv1, nv2)

    phd = PH.PHData(PH.ScalarPenaltyParameter(1.0),
                    st,
                    Dict(1=>copy(PH.scenarios(st))),
                    var_map,
                    TimerOutputs.TimerOutput()
                    )

    for (scid, sinfo) in pairs(phd.scenario_map)

        for vid in keys(sinfo.branch_vars)
            xhid = PH.convert_to_xhat_id(phd, vid)
            sinfo.branch_vars[vid] = var_val[scid][vid]
            phd.xhat[xhid].value = var_val[scid][vid]
        end

        if add_leaves
            # Create entries for leaf variables.  This is normally done at the end of the
            # solve call but since we aren't calling that here...
            for vid in keys(sinfo.leaf_vars)
                xhid = PH.convert_to_xhat_id(phd, vid)
                sinfo.leaf_vars[vid] = var_val[scid][vid]
                phd.xhat[xhid] = PH.HatVariable(var_val[scid][vid], vid, false)
            end
        end
    end

    return phd
end

function fake_worker_info(queue_size::Int=10)

    worker_results = RemoteChannel(()->Channel{PH.Message}(queue_size))
    worker_queue = RemoteChannel(()->Channel{PH.Message}(queue_size))

    return PH.WorkerInf(Dict(1=>worker_queue),
                        worker_results,
                        Dict(1=>Future())
                        )
end


text:
@testset "ChebyshevPoints" begin
    @test_throws Exception ChebyshevPoints{Int}(2)
    @test eltype(@inferred(ChebyshevPoints(2))) == Float64

    for T in Floats
        @test eltype(@inferred(ChebyshevPoints{T}(3))) == T
        @test typeof(@inferred(ChebyshevPoints{T}(3)[1])) == T

        @test ChebyshevPoints{T}(0) ≈ T[]
        @test ChebyshevPoints{T}(1) ≈ T[0]  atol = eps(real(T))
        @test ChebyshevPoints{T}(2) ≈ T[-1,1]
        @test ChebyshevPoints{T}(3) ≈ T[-1,0,1]
        @test ChebyshevPoints{T}(4) ≈ T[-1,-0.5,0.5,1]
        @test ChebyshevPoints{T}(5) ≈ T[-1,-1/sqrt(T(2)),0,1/sqrt(T(2)),1]
    end
end


text:
module BigRationals

import Base.GMP: Limb

export BigRational

"""
    BigRational <: Real

Arbitrary precision rational type.
"""
mutable struct BigRational <: Real
    num_alloc::Cint
    num_size::Cint
    num_d::Ptr{Limb}
    den_alloc::Cint
    den_size::Cint
    den_d::Ptr{Limb}

    function BigRational()
        b = MPQ.init!(new())
        finalizer(cglobal((:__gmpq_clear, :libgmp)), b)
        return b
    end
end

include("MPQ.jl")
include("functions.jl")

end


text:
# The model is the same Rijke tube configuration as in Tutorial 01. However,
# this time the FTF is customly defined. Let's start with the definition of the
# flame transfer function. The custom function must return the FTF and all its
# derivatives w.r.t. the complex freqeuncy. Therefore, the interface **must** be
# a function with two inputs: a complex number and an integer. For an n-τ-model
# it reads.
function FTF(ω::ComplexF64, k::Int=0)
    return n*(-1.0im*τ)^k*exp(-1.0im*ω*τ)
end
# Note, that I haven't defined `n` and `τ` here, but Julia is a quite generous
# programming language; it won't complain if we define them before making our
# first call to the function.

# Of course, it is sometimes hard to find a closed form expression for an FTF
# and all its derivatives. You need to specify at least these derivative orders
# that will be used by the methods applied to analyse the model. This is at
# least the first order derivative for beyn and the standard householder method.
# You might return derivative orders up to the order you need it using an if
# clause for each order. Anyway, the function may get complicated. The implicit
# method explained below is then a viable alternative.

# Now let's set-up the Rijke tube model. The data is the same as in tutorial 1
using WavesAndEigenvalues.Helmholtz
mesh=Mesh("Rijke_mm.msh",scale=0.001) #load mesh
dscrp=Dict() #initialize model descriptor
dscrp["Interior"]=(:interior, ()) #define resonant cavity
dscrp["Outlet"]=(:admittance, (:Y,1E15)) #specify outlet BC
γ=1.4 #ratio of specific heats
ρ=1.225 #density at the reference location upstream to the flame in kg/m^3
Tu=300.0    #K unburnt gas temperature
Tb=1200.0    #K burnt gas temperature
P0=101325.0 # ambient pressure in Pa
A=pi*0.025^2 # cross sectional area of the tube
Q02U0=P0*(Tb/Tu-1)*A*γ/(γ-1) #the ratio of mean heat release to mean velocity Q02U0
x_ref=[0.0; 0.0; -0.00101] #reference point
n_ref=[0.0; 0.0; 1.00] #directional unit vector of reference velocity
R=287.05 # J/(kg*K) specific gas constant (air)
speedofsound(x,y,z) = z<0. ? sqrt(γ*R*Tu) : sqrt(γ*R*Tb)
c=generate_field(mesh,speedofsound)
# btw this is the moment where we specify the values of n and τ
n=1 #interaction index
τ=0.001 #time delay

##
# but instead of passing n and τ and its values to the flame description we
# just pass the FTF
dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref,FTF)) #flame dynamics

#... discretize ...
L=discretize(mesh,dscrp,c)
# and done!
## Check the model
print(L)
# Independently of how you name it, your FTF will be displayed as `FTF`
# in the signature of the operator.

##
# Solving the problem shows that we get the same result as with the built-in
# n-τ-model in tutorial 1
sol,nn,flag=householder(L,340*2*pi,maxiter=20,tol=1E-11)
## changing the parameters
# Because the system parameters n and τ are defined in this outer scope,
# changing them will change the FTF and thus the model without rediscretization.
# For instance, you can completely deactivate the FTF by setting n to 0.
n=0
# Indeed the model now converges to purely acoustic mode
sol,nn,flag=householder(L,340*2*pi,maxiter=20,tol=1E-11)

# However, be aware that there is currently no mechanism keeping track of these
# parameters. It is the programmer's (that's you) responsibility to know the
# specification of the FTF when `sol` was computed.

# Also, note the parameters are defined in this scope. You cannot change it by
# iterating it in a for-loop or from within a function, due to julia's scoping
# rules.

## Something else than n-τ...
# Of course this is an academic example. In practice you will specify
# something very custom as FTF. Maybe a superposition of σ-τ-models or even
# a statespace model fitted to experimental data.

##TODO: , vectorfit
#implicit FTF modelling as in Silva et al ASME2020
dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref)) #no flame dynamics specified at all
H=discretize(mesh,dscrp,c)
sol_base,nn,flag=householder(H,340*2*pi,maxiter=20,tol=1E-11)
perturb_fast!(sol_base,H,:FTF,30)
## Newton-Raphson for finding flame response
ω0=sol_base.params[:ω]
ω(f,k=0)=sol(:FTF,f,15,15)

n=1.0
η0=0
for i=1:40
    Δη=η0-(FTF(ω(η0))-η0)/FTF(ω(η0,1))
    println(Δη)
end


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
    x1'(t) = a * x1(t) - b * x1(t) * x2(t),
    x2'(t) = -c * x2(t) + d * x1(t) * x2(t),
    y1(t) = x1(t)
)

res = res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=0, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
# Index constructors
NoIndex(t::Table) = NoIndex()

function UniqueIndex{names}(t::Table) where {names}
    # TODO check columns actually exist in table
    # TODO add uniqueness check? Maybe a keyword argument or `@checkbounds` or something
    return UniqueIndex{names}()
end

function SortIndex{names}(t::Table) where {names}
    order = sortperm(Project(names)(t))
    return SortIndex{names}(order)
end

function UniqueSortIndex{names}(t::Table) where {names}
    # TODO add uniqueness check? Maybe a keyword argument or `@checkbounds` or something
    order = sortperm(Project(names)(t))
    return UniqueSortIndex{names}(order)
end

function HashIndex{names}(t::Table) where {names}
    t_projected = Project(names)(t)
    d = Dict{_valuetype(eltype(t_projected)), Vector{Int}}()
    for i in keys(t_projected) # TODO make faster
        row = Tuple(t_projected[i])
        if haskey(d, row)
            push!(d[row], i)
        else
            d[row] = Int[i]
        end
    end
    return HashIndex{names}(d)
end

function UniqueHashIndex{names}(t::Table) where {names}
    t_projected = Project(names)(t)
    d = Dict{_valuetype(eltype(t_projected)), Int}()
    for i in keys(t_projected) # TODO make faster
        row = Tuple(t_projected[i])
        if haskey(d, row)
            error("Columns $names do not contain unique values")
        else
            d[row] = i
        end
    end
    return UniqueHashIndex{names}(d)
end

"""
    accelerate(table, AcceleratorType)

Returns a new table with an acceleration index of the specified type. The output table
may be faster at certain operations like `filter` with a relavent `Predicate`, however the
input table will remain unchanged.

See also `accelerate!`.
"""
function accelerate(t::Table{names}, ::Type{I}) where {I <: AbstractIndex, names}
    newindex = I(t)
    return Table{names}(getdata(t), (accelerators(t)..., newindex))
end

"""
    accelerate!(table, AcceleratorType)

Returns a new table with an acceleration index of the specified type. The output table may
be faster at certain operations like `filter` with a `Predicate`. The order of rows in the
input table may be changed in-place to optimize locality and further increase speed - for
example the rows may be ordered by one of the columns. However, the input table will not
have the knowledge of the acceleration index attached - take care to use the output table
after this operation.

See also `accelerate`.
"""
function accelerate!(t::Table{names}, ::Type{<:AbstractIndex})
    return accelerate(t, I)
end

function accelerate!(t::Table{<:Any,<:Any,Tuple{}}, ::Type{SortIndex{names}}) where {names}
    order = sortperm(Project(names)(t))
    foreach(col -> permute!(col, order), getdata(t))
    return Table{names}(getdata(t), (SortIndex{names}(order),))
end

function accelerate!(t::Table{<:Any,<:Any,Tuple{}}, ::Type{UniqueSortIndex{names}}) where {names}
    order = sortperm(Project(names)(t))
    foreach(col -> permute!(col, order), getdata(t))
    return Table{names}(getdata(t), (UniqueSortIndex{names}(order),))
end

# TODO tests for above


text:
using MathOptInterface, JuMP, Dualization, Test

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

const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}
const SQF{T} = MOI.ScalarQuadraticFunction{T}

MOIU.@model(
    TestModel,
    (MOI.ZeroOne, MOI.Integer),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        MOI.Semicontinuous,
        MOI.Semiinteger,
    ),
    (
        MOI.Reals,
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.ExponentialCone,
        MOI.DualExponentialCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.PositiveSemidefiniteConeSquare,
        MOI.RootDetConeTriangle,
        MOI.RootDetConeSquare,
        MOI.LogDetConeTriangle,
        MOI.LogDetConeSquare,
    ),
    (MOI.PowerCone, MOI.DualPowerCone, MOI.SOS1, MOI.SOS2),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Functions that are used inside tests
include("utils.jl")

# Problems database
include("Problems/Linear/linear_problems.jl")
include("Problems/Linear/conic_linear_problems.jl")
include("Problems/Quadratic/quadratic_problems.jl")
include("Problems/SOC/soc_problems.jl")
include("Problems/RSOC/rsoc_problems.jl")
include("Problems/SDP/sdp_triangle_problems.jl")
include("Problems/Exponential/exponential_cone_problems.jl")
include("Problems/Power/power_cone_problems.jl")

# Run tests to travis ci
include("Tests/test_structures.jl")
include("Tests/test_supported.jl")
include("Tests/test_objective_coefficients.jl")
include("Tests/test_dual_model_variables.jl")
include("Tests/test_dual_sets.jl")
include("Tests/test_dualize_conic_linear.jl")
include("Tests/test_dualize_linear.jl")
include("Tests/test_dualize_soc.jl")
include("Tests/test_dualize_rsoc.jl")
include("Tests/test_dualize_sdp.jl")
include("Tests/test_dualize_exponential.jl")
include("Tests/test_dualize_power.jl")
include("Tests/test_dualize_quadratic.jl")
include("Tests/test_dual_names.jl")
include("Tests/test_dot.jl")

include("Tests/test_partial_dual_linear.jl")
include("Tests/test_partial_dual_quadratic.jl")

#=
    Tests depending on solvers
=#

include("optimize_abstract_models.jl")

primal_linear_factory = []
primal_conic_factory = []
primal_power_cone_factory = []

dual_linear_factory = []
dual_conic_factory = []
dual_power_cone_factory = []

dual_linear_optimizer = []
dual_conic_optimizer = []
dual_power_cone_optimizer = []

primal_linear_optimizer = []
primal_conic_optimizer = []
primal_power_cone_optimizer = []

# Load & Test strong duality in linear/conic problems
# Comment the solver that are not available for development
include("Solvers/highs_test.jl")
include("Solvers/csdp_test.jl")
include("Solvers/scs_test.jl")

include("Tests/test_JuMP_dualize.jl")
include("Tests/test_MOI_wrapper.jl")
include("Tests/test_modify.jl")


text:
using CovarianceEstimation
using Statistics
using LinearAlgebra
using Test
using Random
using DelimitedFiles
using StatsBase

include("ref_lw_lshrink.jl")
include("ref_lw_nlshrink.jl")
include("legacy.jl")
Random.seed!(1234)

const CE = CovarianceEstimation

const X   = randn(3, 8)
const Z   = [2 -1 2; -1 2 -1; -1 -1 -1]
const X2s = [[1 0; 0 1; -1 0; 0 -1], [1 0; 0 1; 0 -1; -1 0]]
const Xa  = randn(5, 5)
const Xb  = randn(8, 3)
const Xc  = randn(15, 20)
const Xd  = randn(30, 50)
const Xe  = randn(50, 30)

const test_matrices = [X, Xa, Xb, Xc, Xd, Xe, Z]

function testTransposition(ce::CovarianceEstimator, X)
    @test cov(ce, X; dims=1) ≈ cov(ce, transpose(X); dims=2)
    @test cov(ce, X; dims=2) ≈ cov(ce, transpose(X); dims=1)
end

function testDims(ce::CovarianceEstimator, X)
    @test_throws ArgumentError cov(ce, X; dims=0)
    @test_throws ArgumentError cov(ce, X; dims=3)
end

function testUncorrelated(ce::CovarianceEstimator)
    for X2 ∈ X2s
        @test isdiag(cov(ce, X2))
    end
end

function testTranslation(ce::CovarianceEstimator, X)
    C1 = cov(ce, X)
    C2 = cov(ce, X .+ randn(1, size(X, 2)))
    @test C1 ≈ C2 atol = 1e-12 rtol = 1e-16
    C1t = cov(ce, X')
    C2t = cov(ce, X' .+ randn(1, size(X, 1)))
    @test C1t ≈ C2t atol = 1e-12 rtol = 1e-16
end

function _get_ref(name::AbstractString)
    path = endswith(pwd(), "CovarianceEstimation") ? "test" : ""
    return readdlm(joinpath(path, "test_matrices", "$name.csv"))
end

function _test_ref(ce::CovarianceEstimator, name::AbstractString, ref::AbstractString)
    X = _get_ref(name)
    ref_cov = _get_ref("$(name)_$ref")
    c = cov(ce, X)
    @test c ≈ ref_cov
    @test issymmetric(c)
end

function _test_refs(
    ce::CovarianceEstimator, names::AbstractVector{<:AbstractString}, ref::AbstractString,
)
    for name in names
        _test_ref(ce, name, ref)
    end
end

@testset "CovarianceEstimation" begin
    include("test_biweight.jl")
    include("test_linearshrinkage.jl")
    include("test_nonlinearshrinkage.jl")
end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export mass, set_mass!, charge, set_charge!, name, set_name!, type, set_type!
export vdw_radius, covalent_radius, atomic_number
export set_property!, property, properties_count, list_properties

__ptr(atom::Atom) = __ptr(atom.__handle)
__const_ptr(atom::Atom) = __const_ptr(atom.__handle)

"""
Create an atom with the given ``name`` and set the atom ``type`` to be the same
as ``name``.
"""
function Atom(name::String)
    ptr = @__check_ptr(lib.chfl_atom(pointer(name)))
    return Atom(CxxPointer(ptr, is_const=false))
end

"""
Get a copy of the ``atom`` at the given ``index`` from a ``frame``.
"""
function Atom(frame::Frame, index::Integer)
    ptr = @__check_ptr(lib.chfl_atom_from_frame(__ptr(frame), UInt64(index)))
    atom = Atom(CxxPointer(ptr, is_const=false))
    copy = deepcopy(atom)
    finalize(atom)
    return copy
end

"""
Get a copy of the ``atom`` at the given ``index`` from a ``topology``.
"""
function Atom(topology::Topology, index::Integer)
    ptr = @__check_ptr(lib.chfl_atom_from_topology(__ptr(topology), UInt64(index)))
    atom = Atom(CxxPointer(ptr, is_const=false))
    copy = deepcopy(atom)
    finalize(atom)
    return copy
end

"""
Get the mass of an ``atom`` in atomic mass units.
"""
function mass(atom::Atom)
    result = Ref{Float64}(0)
    __check(lib.chfl_atom_mass(__const_ptr(atom), result))
    return result[]
end

"""
Set the mass of an ``atom`` to ``mass``.

The mass must be in atomic mass units.
"""
function set_mass!(atom::Atom, mass)
    __check(lib.chfl_atom_set_mass(__ptr(atom), Float64(mass)))
    return nothing
end

"""
Get the charge of an ``atom`` in number of the electron charge *e*.
"""
function charge(atom::Atom)
    result = Ref{Float64}(0)
    __check(lib.chfl_atom_charge(__const_ptr(atom), result))
    return result[]
end

"""
Set the charge of an ``atom`` to ``charge``.

The charge must be in number of the electron charge *e*.
"""
function set_charge!(atom::Atom, charge)
    __check(lib.chfl_atom_set_charge(__ptr(atom), Float64(charge)))
    return nothing
end

"""
Get the name of an ``atom``.
"""
function name(atom::Atom)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_atom_name(__const_ptr(atom), buffer, size))
    )
end

"""
Set the name of an ``atom`` to ``name``.
"""
function set_name!(atom::Atom, name::String)
    __check(lib.chfl_atom_set_name(__ptr(atom), pointer(name)))
    return nothing
end

"""
Get the type of an ``atom``.
"""
function type(atom::Atom)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_atom_type(__const_ptr(atom), buffer, size))
    )
end

"""
Set the type of an ``atom`` to ``type``.
"""
function set_type!(atom::Atom, type::String)
    __check(lib.chfl_atom_set_type(__ptr(atom), pointer(type)))
    return nothing
end

"""
Get the full name of an ``atom`` from the atom type.

For example, the full name of an atom with type "He" is "Helium".
"""
function Base.fullname(atom::Atom)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_atom_full_name(
            __const_ptr(atom), buffer, size)
        )
    )
end

"""
Get the van der Waals radius of an ``atom`` from the atom type.

If the radius can not be found, this function returns 0.
"""
function vdw_radius(atom::Atom)
    radius = Ref{Float64}(0)
    __check(lib.chfl_atom_vdw_radius(__const_ptr(atom), radius))
    return radius[]
end

"""
Get the covalent radius of an ``atom`` from the atom type.

If the radius can not be found, returns 0.
"""
function covalent_radius(atom::Atom)
    radius = Ref{Float64}(0)
    __check(lib.chfl_atom_covalent_radius(__const_ptr(atom), radius))
    return radius[]
end

"""
Get the atomic number of an ``atom`` from the atom type.

If the atomic number can not be found, returns 0.
"""
function atomic_number(atom::Atom)
    number = Ref{UInt64}(0)
    __check(lib.chfl_atom_atomic_number(__const_ptr(atom), number))
    return number[]
end

"""
Set a named property for the given atom.
"""
function set_property!(atom::Atom, name::String, value)
    property = Property(value)
    __check(lib.chfl_atom_set_property(
        __ptr(atom), pointer(name), __const_ptr(property)
    ))
    return nothing
end

"""
Get a named property for the given atom.
"""
function property(atom::Atom, name::String)
    ptr = lib.chfl_atom_get_property(__const_ptr(atom), pointer(name))
    return extract(Property(CxxPointer(ptr, is_const=false)))
end

"""
Get the number of properties associated with an atom.
"""
function properties_count(atom::Atom)
    count = Ref{UInt64}(0)
    __check(lib.chfl_atom_properties_count(__const_ptr(atom), count))
    return count[]
end

"""
Get the names of all properties associated with an atom.
"""
function list_properties(atom::Atom)
    count = properties_count(atom)
    names = Array{Ptr{UInt8}}(undef, count)
    __check(lib.chfl_atom_list_properties(__const_ptr(atom), pointer(names), count))
    return map(unsafe_string, names)
end

"""
Make a deep copy of an ``atom``.
"""
function Base.deepcopy(atom::Atom)
    ptr = lib.chfl_atom_copy(__const_ptr(atom))
    return Atom(CxxPointer(ptr, is_const=false))
end


text:
# JuLIP.jl master file.

__precompile__()
module JuLIP

# warn("""This is a major breaking change for JuLIP, effectively removing
#         the dependency on ASE. If you have been using the JuLIP master branch
#         until now, then please checkout the branch `v0.2ase`. Access to ASE
#         will be restored asap, in a separate package. Documentation
#         on how to use the 'new JuLIP' will hopefully follow very soon.""")

using Reexport, NeighbourLists, StaticArrays, Parameters

# quickly switch between Matrices and Vectors of SVectors, etc
include("arrayconversions.jl")

# define types and abstractions of generic functions
include("abstractions.jl")

include("chemistry.jl")
@reexport using JuLIP.Chemistry

# the main atoms type
include("atoms.jl")

# how to build some simple domains
include("build.jl")
@reexport using JuLIP.Build

# a few auxiliary routines
include("utils.jl")

# interatomic potentials prototypes and some example implementations
include("Potentials.jl")
@reexport using JuLIP.Potentials

# submodule JuLIP.Constraints
include("Constraints.jl")
@reexport using JuLIP.Constraints

# basic preconditioning capabilities
include("preconditioners.jl")
@reexport using JuLIP.Preconditioners

# some solvers
include("Solve.jl")
@reexport using JuLIP.Solve

# experimental features
include("Experimental.jl")
@reexport using JuLIP.Experimental

# codes to facilitate testing
include("Testing.jl")

# # only try to import Visualise, it is not needed for the rest to work.
# try
#    # some visualisation options
#    if isdefined(Main, :JULIPVISUALISE)
#       if Main.JULIPVISUALISE == true
#          include("Visualise.jl")
#       end
#    end
# catch
#    JuLIP.julipwarn("""JuLIP.Visualise did not import correctly, probably because
#                `imolecule` is not correctly installed.""")
# end

end # module