text:
module mFEM
using Printf
using ProgressMeter

struct Mesh
  P::Array # be an information matrix consisting of the coordinates of all mesh nodes
  T::Array # be an information matrix consisting of the global node indices
  BNode::Array # be an information matrix consisting to the node indices of boundary nodes
                # e.g. [0,2,3]
  BEage::Array # e.g. [0 1;3 4;2 5], with counterclockwise
  BSurface::Array # e.g. [0 1 3; 2 3 5] with counterclockwise
  Integral::Function # (u::Function,n::Int)::Float64
                      # calc integral on the nth element of mesh
                      # u: (x::Tuple)::Float64
end
export Mesh

export Basis
struct Basis
  Nlb::Int # the number of local trial basis functions
  Pb::Array # the nodes corresponding to the finite element basis functions.
  Tb::Array # the nodes corresponding to the finite element basis functions.
  BNode::Array # be an information matrix consisting to the node indices of boundary nodes
                # e.g. [0,2,3]
  BEage::Array # e.g. [0 1;3 4;2 5], with counterclockwise
  BSurface::Array # e.g. [0 1 3; 2 3 5] with counterclockwise
  Dϕ::Function  # (k::Int, x::Tuple, n::Int, α::Int ) -> Real
                # return the kth deriavtive of (n,α) in mesh basis function at x
end

export stiffness_mat,load_vec,set_boundary_cond!
# - D(cDu) = f in Ω, f=g on ∂Ω
# c(x::Tuple)
function stiffness_mat(mesh::mFEM.Mesh,trial::mFEM.Basis,test::mFEM.Basis,c::Function)

  A = zeros(size(test.Pb,1),size(trial.Pb,1)) # TODO sparse matrix

  @showprogress "stiffness_mat " for n = 1:size(mesh.T,1)
    for α = 1:trial.Nlb, β = 1:test.Nlb

      r = mesh.Integral(n) do x::Tuple
        return c(x) * trial.Dϕ(1,x,n,α) .* test.Dϕ(1,x,n,β)
      end

      i = test.Tb[n,β]
      j = trial.Tb[n,α]

      A[i,j] += r
    end
  end

  return A
end

# f(x::Tuple)
function load_vec(mesh::mFEM.Mesh,test::mFEM.Basis,f::Function)

  b = zeros(size(test.Pb,1))
  @showprogress "load_vec " for n = 1:size(mesh.T,1)
    for  β = 1:test.Nlb

      r = mesh.Integral(n) do x::Tuple
        return f(x) * test.Dϕ(0,x,n,β)
      end

      i = test.Tb[n,β]

      b[i] += r
    end
  end
  return b
end

# g(x::Tuple, boundary_type::Int)
function set_boundary_cond!(trial::mFEM.Basis,boundarynodes::AbstractArray,
    A::AbstractArray,b::AbstractArray,g::Function)
  nA = size(A,2)
  nbn = size(boundarynodes,1)
  for nb = 1:nbn
    type,i = boundarynodes[nb,:]
    if type == -1 # Dirichlet
      A[i,:] = zeros(nA)
      A[i,i] = 1
      b[i] = g(tuple(trial.Pb[i,:]...),type)
    else
      error("unsupported type")
    end
  end
end

# submodules begin
include("./dim1.jl")
include("./dim2.jl")
# submodules end

# for test

# test()::Float64, error of method
function convergence_order(test::Function;Oh::Bool = true)
  N = [2^i for i=2:7]
  Error = zeros(Float64,size(N,1))
  Order = zeros(Float64,size(N,1))

  for i in 1:size(N,1)
    n = N[i]
    Error[i] = test(n)
  end

  for i in size(Error,1):-1:2
    Order[i] = Error[i-1] / Error[i]
  end

  println("   n\terror\t\torder")
  for i in 1:size(N,1)
    @printf "%4d %.10e\t%.2f" N[i] Error[i] Order[i]
    if Oh && Order[i] > 0
      @printf "\tO(h^%.6f)" log(2,Order[i])
    end
    @printf "\n"
  end
end

end # mFEM


text:
using NiGraphEmbedding
import LinearAlgebra, Statistics
using Test, Random
using NiLang, NiLang.AD
using Zygote
using Optim

@testset "gradient" begin
    Random.seed!(2)
    params = randn(5, 10)
    x = randn(5,10)
    @test embedding_loss(0.0, x)[1] ≈ embedding_loss(x)
    @test grad(Grad(embedding_loss)(Val(1), 0.0, x)[3]) ≈ Zygote.gradient(embedding_loss, x)[1]
end

@testset "train" begin
    Random.seed!(2)
    params = randn(5, 10)
    train(params, Adam(lr=0.01))
    # distances of connected bonds.
    d1s = [LinearAlgebra.norm(params[:,i]-params[:,j]) for (i,j) in NiGraphEmbedding.L1]
    # distances of disconnected bonds.
    d2s = [LinearAlgebra.norm(params[:,i]-params[:,j]) for (i,j) in NiGraphEmbedding.L2]
    @test abs(Statistics.var(d1s)) < 1e-5
    @test abs(Statistics.var(d2s)) < 1e-5

    params = randn(5, 10)
    opt = NewtonTrustRegion(; initial_delta = 1.0,
                    delta_hat = 20.0,
                    eta = 0.1,
                    rho_lower = 0.25,
                    rho_upper = 0.75)
    train(params, opt)

    d1s = [LinearAlgebra.norm(params[:,i]-params[:,j]) for (i,j) in NiGraphEmbedding.L1]
    # distances of disconnected bonds.
    d2s = [LinearAlgebra.norm(params[:,i]-params[:,j]) for (i,j) in NiGraphEmbedding.L2]
    @test abs(Statistics.var(d1s)) < 1e-5
    @test abs(Statistics.var(d2s)) < 1e-5
end

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_dashdraggableresponsive

"""
    ''_dashdraggableresponsive(;kwargs...)
    ''_dashdraggableresponsive(children::Any;kwargs...)
    ''_dashdraggableresponsive(children_maker::Function;kwargs...)

A DashDraggableResponsive component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `value` (String; optional): The value displayed in the input.
- `layout` (Array of Dicts; optional)
- `ncols` (Real; optional)
- `nrows` (Real; optional)
"""
function ''_dashdraggableresponsive(; kwargs...)
        available_props = Symbol[:children, :id, :value, :layout, :ncols, :nrows]
        wild_props = Symbol[]
        return Component("''_dashdraggableresponsive", "DashDraggableResponsive", "dash_draggable", available_props, wild_props; kwargs...)
end

''_dashdraggableresponsive(children::Any; kwargs...) = ''_dashdraggableresponsive(;kwargs..., children = children)
''_dashdraggableresponsive(children_maker::Function; kwargs...) = ''_dashdraggableresponsive(children_maker(); kwargs...)

max_stars_repo_name: tcgvn/dash-draggable

max_issues_repo_name: tcgvn/dash-draggable

max_forks_repo_name: tcgvn/dash-draggable


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

using Markdown
using InteractiveUtils

# ╔═╡ b562e4c8-d339-11ea-3903-a936ef2608e6
using JLD2

# ╔═╡ e58fb87e-d339-11ea-15e9-755175e9db02
using Plots

# ╔═╡ 5f25df2c-d33e-11ea-2349-b3cde1716f8f
using LaTeXStrings

# ╔═╡ 9ba366a6-d33f-11ea-294e-57f4c58b40c0
md"# Load saved data"

# ╔═╡ e93aba66-d339-11ea-0296-b360f581dba3
begin
	data = Dict{String,Any}()
	jldopen( "Nflation_example.jld2", "r") do f
		for k in keys(f)
			data[k] = read(f,k)
		end
	end
	data
end

# ╔═╡ adca77d4-d33f-11ea-2a01-31ef55bf0781
md"# Background evolution"

# ╔═╡ 477e130e-d33b-11ea-2c9b-8300969f6553
begin
	l = @layout [a b ; c d]
	ehp = plot(data["back_Ne"],data["eh"], ylabel=L"\epsilon_H", xlabel=L"N_e")
	evp = plot(data["back_Ne"],data["ev"], ylabel=L"\epsilon_V", xlabel=L"N_e")
	etp = plot(data["back_Ne"],data["et"], ylabel=L"\eta_H", xlabel=L"N_e")
	omp = plot(data["back_Ne"],data["omega"], ylabel=L"\omega^2/H^2", xlabel=L"N_e")
	plot(ehp,evp,etp,omp,layout=l, legend=false, linewidth=2)
	#plot(ehp,evp,omp,layout=l, legend=false, linewidth=2)
end

# ╔═╡ a6dc05f4-d340-11ea-22a7-39b8356c7cd8
md"# Planck compatibility?"

# ╔═╡ 57e277c6-d345-11ea-2cc1-4fb4bfc0be6f
md"This is a rough plot based on the Planck 2018 data release, real planck data should be used if you want to be more careful."

# ╔═╡ b6b376d8-d340-11ea-2d9f-2f68ab41a13b
begin
	planck_ns = 0.9668
	planck_ns_std = 0.0037
	planck_r = 0
	planck_r_std = 0.058/1.95996

	num_sigma = 5
	function ellipseShape(h,k,a,b)
		θ = range(0,2π,length=300)
		h .+ a*cos.(θ), k .+ b*sin.(θ)
	end
	plot(ellipseShape(planck_ns,planck_r,planck_ns_std,planck_r_std), seriestype = [:shape,], lw=0.5, c = :red, linecolor=:black, legend=false, fillalpha = 0.2,label="Planck 2018 (up to $(num_sigma)σ)")
	for i in 2:num_sigma
		plot!(ellipseShape(planck_ns,planck_r,i*planck_ns_std,i*planck_r_std), seriestype = [:shape,], lw=0.5, c = :red, linecolor=:black, legend=false, fillalpha = 0.2,label="")
	end
	scatter!([data["ns"]],[data["r"]],color=:green,leg=true,label="This model")
	plot!(ylim=[0,:auto])
end

# ╔═╡ a2f15f0e-d352-11ea-2f7b-e16f13cbdd76
md"# Sub and super-horizon evolution of 2-pt functions"

# ╔═╡ 79dabff8-d356-11ea-0363-d1d94cdbc5c1
md"The powerspectra are plotted starting when Bunch-Davies initial conditions are imposed, and horizon exit is marked as a vertical red line"

# ╔═╡ b6bb0206-d352-11ea-3110-63c5067874a8
begin
	Nend = data["pert_Ne"][end]
	plot(data["pert_Ne"] .- Nend, data["PzN"],label=L"P_\zeta",lw=2)
	plot!(data["pert_Ne"] .- Nend, data["PisoN"],label=L"P_\mathrm{iso}",lw=2)
	plot!(data["pert_Ne"] .- Nend, data["rN"].*data["PzN"],label=L"P_\mathrm{tens}",lw=2)
	plot!(yscale=:log10)
	vline!([-55.0],color=:red,lw=2,linestyle=:dash,label="")
	plot!(xlabel=L"N_e",ylabel=L"P_{\circ}(k_\star,N)")
end

# ╔═╡ Cell order:
# ╠═b562e4c8-d339-11ea-3903-a936ef2608e6
# ╠═e58fb87e-d339-11ea-15e9-755175e9db02
# ╠═5f25df2c-d33e-11ea-2349-b3cde1716f8f
# ╟─9ba366a6-d33f-11ea-294e-57f4c58b40c0
# ╠═e93aba66-d339-11ea-0296-b360f581dba3
# ╟─adca77d4-d33f-11ea-2a01-31ef55bf0781
# ╠═477e130e-d33b-11ea-2c9b-8300969f6553
# ╟─a6dc05f4-d340-11ea-22a7-39b8356c7cd8
# ╟─57e277c6-d345-11ea-2cc1-4fb4bfc0be6f
# ╠═b6b376d8-d340-11ea-2d9f-2f68ab41a13b
# ╟─a2f15f0e-d352-11ea-2f7b-e16f13cbdd76
# ╟─79dabff8-d356-11ea-0363-d1d94cdbc5c1
# ╠═b6bb0206-d352-11ea-3110-63c5067874a8


text:
function on_gcloud(;return_name::Bool=false,hostname::String="instance")
    # return if we're on google cloud, i.e. the word 'instance' is in the hostname
    if !("hostname" in readdir("/etc"))
        # if the "hostname" file doesn't exist we're not on gcloud
        return false
    end
    hname = open("/etc/hostname","r") do f
        readstring(f)
    end
    if return_name
        return hname
    end
    if contains(hname,hostname)
        return true
    else
        return false
    end
end
function send_mail(subject::String,body::String;to::String="brett.israelsen@gmail.com")
    # send an email via mutt, mutt has to be install and configured properly to do this.
    # this site might be helpful in doing that nickdesauliniers.github.io/blog/2016/06/18/mutt-gmail-ubuntu/
    run(pipeline(`echo $body`,`mutt -s $body $to`))
    println("sent mail")
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using DataFrames

import ...ComboLeg,
       ...CommissionReport,
       ...ConditionType,
       ...Contract,
       ...ContractDescription,
       ...ContractDetails,
       ...DeltaNeutralContract,
       ...Execution,
       ...FamilyCode,
       ...FaDataType,
       ...MarketDataType,
       ...Order,
       ...OrderState,
       ...SoftDollarTier,
       ...TickAttrib,
       ...TickAttribLast,
       ...TickAttribBidAsk,
       ...condition_map,
       ...ns

"""
    slurp(::Type{T}, it)

Utility functions to read from an iterator `it` and convert to types `T`.
"""
slurp(::Type{T}, it) where T<:Union{Bool,Int,Enum{Int32},Float64,String} = convert(T, pop(it))

slurp(::Type{T}, it) where T = T(take(it, fieldcount(T))...)

slurp(t, it) = convert.(t, take(it, length(t)))

slurp!(x::T, idx, it) where T = setfield!.(Ref(x), idx, slurp(fieldtype.(T, idx), it))

"""
    tagvalue2nt(x)

Convert a `String[]` into a `NamedTuple` like this:

    ["tag1", "value1", "tag2", "value2", ...] -> (tag1="value1", tag2="value2", ...)
"""
function tagvalue2nt(x)

  @assert iseven(length(x))

  s = collect(String, x)

  (; (Symbol(t) => v for (t, v) ∈ Iterators.partition(s, 2))...)
end

function unmask(T::Type{NamedTuple{M,NTuple{N,Bool}}}, mask) where {M,N}

  a = digits(Bool, mask, base=2, pad=N)

  length(a) == N || @error "unmask(): wrong attribs" T mask

  T(a)
end

function fill_df(eltypes, names, n, it)

  df = DataFrame([Vector{T}(undef, n) for T ∈ eltypes], names; copycols=false)

  for r ∈ 1:nrow(df), c ∈ 1:ncol(df)
    df[r, c] = pop(it)
  end

  df
end

"""
    process::Dict{Int,Function}

Collection of parsers indexed by message ID
"""
const process = Dict{Int,Function}(    # TODO Use a Tuple instead?

  # TICK_PRICE
   1 => function(it, w, ver)

          tickerId::Int,
          ticktype::Int,
          price::Float64,
          size::Float64,
          mask::Int = it

          w.tickPrice(tickerId, tickname(ticktype), price, size, unmask(TickAttrib, mask))
        end,

  # TICK_SIZE
   2 => function(it, w, ver)

          tickerId::Int,
          ticktype::Int,
          size::Float64 = it

          w.tickSize(tickerId, tickname(ticktype), size)
        end,

  # ORDER_STATUS
   3 => (it, w, ver) -> w.orderStatus(slurp((Int,String,Float64,Float64,Float64,Int,Int,Float64,Int,String,Float64), it)...),

  # ERR_MSG
   4 => (it, w, ver) -> w.error(slurp((Int,Int,String), it)...),

  # OPEN_ORDER
   5 => function(it, w, ver)

          o = Order()
          c = Contract()

          o.orderId = pop(it)

          slurp!(c, [1:8; 10:12], it)

          slurp!(o, 4:9, it)  # :action through :tif

          slurp!(o, (:ocaGroup,
                     :account,
                     :openClose,
                     :origin,
                     :orderRef,
                     :clientId,
                     :permId,
                     :outsideRth,
                     :hidden,
                     :discretionaryAmt,
                     :goodAfterTime), it)

          pop(it)    # Deprecated sharesAllocation

          slurp!(o, (:faGroup,
                     :faMethod,
                     :faPercentage,
                     :faProfile,
                     :modelCode,
                     :goodTillDate,
                     :rule80A,
                     :percentOffset,
                     :settlingFirm,
                     :shortSaleSlot,
                     :designatedLocation,
                     :exemptCode), it)

          slurp!(o, 43:48, it)    # :auctionStrategy through :stockRangeUpper

          slurp!(o, (:displaySize,
                     :blockOrder,
                     :sweepToFill,
                     :allOrNone,
                     :minQty,
                     :ocaType), it)

          pop(it)    # Deprecated eTradeOnly
          pop(it)    # Deprecated firmQuoteOnly
          pop(it)    # Deprecated nbboPriceCap

          slurp!(o, (:parentId,
                     :triggerMethod), it)

          slurp!(o, 51:54, it)    # :volatility through :deltaNeutralAuxPrice

          !isempty(o.deltaNeutralOrderType) && slurp!(o, 55:62, it)  # :deltaNeutralConId through :deltaNeutralDesignatedLocation

          slurp!(o, (:continuousUpdate,
                     :referencePriceType,
                     :trailStopPrice,
                     :trailingPercent,
                     :basisPoints,
                     :basisPointsType), it)

          c.comboLegsDescrip = pop(it)

          # ComboLegs
          n::Int = pop(it)

          for _ ∈ 1:n
            push!(c.comboLegs, slurp(ComboLeg, it))
          end

          # OrderComboLeg
          n = pop(it)
          append!(o.orderComboLegs, take(it, n))

          # SmartComboRouting
          n = pop(it)
          n > 0 && (o.smartComboRoutingParams = tagvalue2nt(take(it, 2n)))

          slurp!(o, (:scaleInitLevelSize,
                     :scaleSubsLevelSize,
                     :scalePriceIncrement), it)

          !isnothing(o.scalePriceIncrement) &&
          o.scalePriceIncrement > 0         && slurp!(o, 70:76, it) # scalePriceAdjustValue through scaleRandomPercent

          o.hedgeType = pop(it)

          !isempty(o.hedgeType) && (o.hedgeParam = pop(it))

          slurp!(o, (:optOutSmartRouting,
                     :clearingAccount,
                     :clearingIntent,
                     :notHeld), it)

          # DeltaNeutralContract
          slurp(Bool, it) && (c.deltaNeutralContract = slurp(DeltaNeutralContract, it))

          # AlgoStrategy
          o.algoStrategy = pop(it)

          if !isempty(o.algoStrategy)
            n = pop(it)
            n > 0 && (o.algoParams = tagvalue2nt(take(it, 2n)))
          end

          o.solicited,
          o.whatIf = it

          os = OrderState(take(it, 15)..., ns, ns)

          o.randomizeSize,
          o.randomizePrice = it

          o.orderType == "PEG BENCH" && slurp!(o, (:referenceContractId,
                                                   :isPeggedChangeAmountDecrease,
                                                   :peggedChangeAmount,
                                                   :referenceChangeAmount,
                                                   :referenceExchangeId), it)

          # Conditions
          n = pop(it)

          if n > 0
            for _ ∈ 1:n
              push!(o.conditions, slurp(condition_map[slurp(ConditionType, it)], it))
            end

            o.conditionsIgnoreRth,
            o.conditionsCancelOrder = it
          end

          slurp!(o, (:adjustedOrderType,
                     :triggerPrice,
                     :trailStopPrice,
                     :lmtPriceOffset,
                     :adjustedStopPrice,
                     :adjustedStopLimitPrice,
                     :adjustedTrailingAmount,
                     :adjustableTrailingUnit), it)

          o.softDollarTier = slurp(SoftDollarTier, it)

          slurp!(o, (:cashQty,
                     :dontUseAutoPriceForHedge,
                     :isOmsContainer,
                     :discretionaryUpToLimitPrice,
                     :usePriceMgmtAlgo), it)

          ver ≥ Client.DURATION && (o.duration = pop(it))

          ver ≥ Client.POST_TO_ATS && (o.postToAts = pop(it))

          ver ≥ Client.AUTO_CANCEL_PARENT && (o.autoCancelParent = pop(it))

          w.openOrder(o.orderId, c, o, os)
        end,

  # ACCT_VALUE
   6 => (it, w, ver) -> w.updateAccountValue(collect(String, take(it, 4))...),

  # PORTFOLIO_VALUE
   7 => function(it, w, ver)

          c = Contract()

          slurp!(c, [1:7; 9:12], it)

          w.updatePortfolio(c, collect(Float64, take(it, 6))..., slurp(String, it))
        end,

  # ACCT_UPDATE_TIME
   8 => (it, w, ver) -> w.updateAccountTime(slurp(String, it)),

  # NEXT_VALID_ID
   9 => (it, w, ver) -> w.nextValidId(slurp(Int, it)),

  # CONTRACT_DATA
  10 => function(it, w, ver)

          reqId::Int = pop(it)

          cd = ContractDetails()

          slurp!(cd.contract, [2:6; 8; 10; 11], it)

          cd.marketName,
          cd.contract.tradingClass,
          cd.contract.conId,
          cd.minTick,
          cd.mdSizeMultiplier,
          cd.contract.multiplier = it

          slurp!(cd, 4:8, it)
          cd.contract.primaryExchange = pop(it)
          slurp!(cd, 9:17, it)

          n::Int = pop(it)
          n > 0 && (cd.secIdList = tagvalue2nt(take(it, 2n)))

          slurp!(cd, (:aggGroup,
                      :underSymbol,
                      :underSecType,
                      :marketRuleIds,
                      :realExpirationDate,
                      :stockType), it)

          ver ≥ Client.FRACTIONAL_SIZE_SUPPORT && (cd.sizeMinTick = pop(it))

          w.contractDetails(reqId, cd)
        end,

  # EXECUTION_DATA
  11 => function(it, w, ver)

          reqId::Int,
          orderId = it

          c = Contract()
          slurp!(c, [1:8; 10:12], it)

          e = Execution(orderId, take(it, 17)...)

          w.execDetails(reqId, c, e)
        end,

  # MARKET_DEPTH
  12 => (it, w, ver) -> w.updateMktDepth(slurp((Int,Int,Int,Int,Float64,Float64), it)...),

  # MARKET_DEPTH_L2
  13 => (it, w, ver) -> w.updateMktDepthL2(slurp((Int,Int,String,Int,Int,Float64,Float64,Bool), it)...),

  # NEWS_BULLETINS
  14 => (it, w, ver) -> w.updateNewsBulletin(slurp((Int,Int,String,String), it)...),

  # MANAGED_ACCTS
  15 => (it, w, ver) -> w.managedAccounts(slurp(String, it)),

  # RECEIVE_FA
  16 => (it, w, ver) -> w.receiveFA(slurp((FaDataType,String), it)...),

  # HISTORICAL_DATA
  17 => function(it, w, ver)

          reqId::Int = pop(it)

          pop(it) # Ignore startDate
          pop(it) # Ignore endDate

          n::Int = pop(it)

          df = fill_df([String,Float64,Float64,Float64,Float64,Float64,Float64,Int],
                       [:time, :open, :high, :low, :close, :volume, :wap, :count],
                       n, it)

          w.historicalData(reqId, df)
        end,

  # BOND_CONTRACT_DATA
  18 => function(it, w, ver)

          reqId::Int = pop(it)

          cd = ContractDetails()

          slurp!(cd.contract, 2:3, it)

          slurp!(cd, (:cusip,
                      :coupon,
                      :maturity,
                      :issueDate,
                      :ratings,
                      :bondType,
                      :couponType,
                      :convertible,
                      :callable,
                      :putable,
                      :descAppend), it)

          cd.contract.exchange,
          cd.contract.currency,
          cd.marketName,
          cd.contract.tradingClass,
          cd.contract.conId = it

          slurp!(cd, (:minTick,
                      :mdSizeMultiplier,
                      :orderTypes,
                      :validExchanges,
                      :nextOptionDate,
                      :nextOptionType,
                      :nextOptionPartial,
                      :notes,
                      :longName,
                      :evRule,
                      :evMultiplier), it)

          n::Int = pop(it)
          n > 0 && (cd.secIdList = tagvalue2nt(take(it, 2n)))

          cd.aggGroup,
          cd.marketRuleIds = it

          w.bondContractDetails(reqId, cd)
        end,

  # SCANNER_PARAMETERS
  19 => (it, w, ver) -> w.scannerParameters(slurp(String, it)),

  # SCANNER_DATA
  20 => function(it, w, ver)

          tickerId::Int,
          n::Int = it

          rank = Vector{Int}(undef, n)
          cd = [ContractDetails() for _ ∈ 1:n]
          distance =   Vector{String}(undef, n)
          benchmark =  Vector{String}(undef, n)
          projection = Vector{String}(undef, n)
          legsStr =    Vector{String}(undef, n)

          for i ∈ 1:n

            rank[i] = pop(it)

            slurp!(cd[i].contract, [1:6; 8; 10; 11], it)

            cd[i].marketName,
            cd[i].contract.tradingClass,
            distance[i],
            benchmark[i],
            projection[i],
            legsStr[i] = it
          end

          w.scannerData(tickerId, rank, cd, distance, benchmark, projection, legsStr)
        end,

  # TICK_OPTION_COMPUTATION
  21 => function(it, w, ver)

          tickerId::Int,
          ticktype::Int = it

          tickAttrib = if ver ≥ Client.PRICE_BASED_VOLATILITY
                         slurp(Int, it)
                       end

          v = collect(Union{Float64,Nothing}, take(it, 8))

          # (impliedVol, optPrice, pvDividend, undPrice) == -1 means NA
          replace!(view(v, [1, 3, 4, 8]), -1 => nothing)

          # (delta, gamma, vega, theta) == -2 means NA
          replace!(view(v, [2, 5, 6, 7]), -2 => nothing)

          w.tickOptionComputation(tickerId, tickname(ticktype), tickAttrib, v...)
        end,

  # TICK_GENERIC
  45 => function(it, w, ver)

          tickerId::Int,
          ticktype::Int,
          value::Float64 = it

          w.tickGeneric(tickerId, tickname(ticktype), value)
        end,

  # TICK_STRING
  46 => function(it, w, ver)

          tickerId::Int,
          ticktype::Int,
          value::String = it

          w.tickString(tickerId, tickname(ticktype), value)
        end,

  # TICK_EFP
  47 => function(it, w, ver)

          tickerId::Int,
          ticktype::Int = it

          w.tickEFP(tickerId, tickname(ticktype), slurp((Float64,String,Float64,Int,String,Float64,Float64), it)...)
        end,

  # CURRENT_TIME
  49 => function(it, w, ver)

          # TODO: Convert to [Zoned]DateTime
          w.currentTime(slurp(Int, it))
        end,

  # REAL_TIME_BARS
  50 => (it, w, ver) -> w.realtimeBar(slurp((Int,Int,Float64,Float64,Float64,Float64,Float64,Float64,Int), it)...),

  # FUNDAMENTAL_DATA
  51 => (it, w, ver) -> w.fundamentalData(slurp((Int,String), it)...),

  # CONTRACT_DATA_END
  52 => (it, w, ver) -> w.contractDetailsEnd(slurp(Int, it)),

  # OPEN_ORDER_END
  53 => (it, w, ver) -> w.openOrderEnd(),

  # ACCT_DOWNLOAD_END
  54 => (it, w, ver) -> w.accountDownloadEnd(slurp(String, it)),

  # EXECUTION_DATA_END
  55 => (it, w, ver) -> w.execDetailsEnd(slurp(Int, it)),

  # DELTA_NEUTRAL_VALIDATION
  56 => function(it, w, ver)

          reqId::Int = pop(it)

          w.deltaNeutralValidation(reqId, slurp(DeltaNeutralContract, it))
        end,

  # TICK_SNAPSHOT_END
  57 => (it, w, ver) -> w.tickSnapshotEnd(slurp(Int, it)),

  # MARKET_DATA_TYPE
  58 => (it, w, ver) -> w.marketDataType(slurp((Int,MarketDataType), it)...),

  # COMMISSION_REPORT
  59 => (it, w, ver) -> w.commissionReport(slurp(CommissionReport, it)),

  # POSITION_DATA
  61 => function(it, w, ver)

          account::String = pop(it)

          c = Contract()
          slurp!(c, [1:8; 10:12], it)

          w.position(account, c, collect(Float64, take(it, 2))...)
        end,

  # POSITION_END
  62 => (it, w, ver) -> w.positionEnd(),

  # ACCOUNT_SUMMARY
  63 => function(it, w, ver)

          reqId::Int = pop(it)

          w.accountSummary(reqId, collect(String, take(it, 4))...)
        end,

  # ACCOUNT_SUMMARY_END
  64 => (it, w, ver) -> w.accountSummaryEnd(slurp(Int, it)),

  # VERIFY_MESSAGE_API
  65 => (it, w, ver) -> w.verifyMessageAPI(slurp(String, it)),

  # VERIFY_COMPLETED
  66 => (it, w, ver) -> w.verifyCompleted(slurp((Bool,String), it)...),

  # DISPLAY_GROUP_LIST
  67 => (it, w, ver) -> w.displayGroupList(slurp((Int,String), it)...),

  # DISPLAY_GROUP_UPDATED
  68 => (it, w, ver) -> w.displayGroupUpdated(slurp((Int,String), it)...),

  # VERIFY_AND_AUTH_MESSAGE_API
  69 => (it, w, ver) -> w.verifyAndAuthMessageAPI(collect(String, take(it, 2))...),

  # VERIFY_AND_AUTH_COMPLETED
  70 => (it, w, ver) -> w.verifyAndAuthCompleted(slurp((Bool,String), it)...),

  # POSITION_MULTI
  71 => function(it, w, ver)

          reqId::Int,
          account::String = it

          c = Contract()
          slurp!(c, [1:8; 10:12], it)

          position::Float64,
          avgCost::Float64,
          modelCode::String = it

          w.positionMulti(reqId, account, modelCode, c, position, avgCost)
        end,

  # POSITION_MULTI_END
  72 => (it, w, ver) -> w.positionMultiEnd(slurp(Int, it)),

  # ACCOUNT_UPDATE_MULTI
  73 => (it, w, ver) -> w.accountUpdateMulti(slurp(Int, it), collect(String, take(it, 5))...),

  # ACCOUNT_UPDATE_MULTI_END
  74 => (it, w, ver) -> w.accountUpdateMultiEnd(slurp(Int, it)),

  # SECURITY_DEFINITION_OPTION_PARAMETER
  75 => function(it, w, ver)

          args = slurp((Int,String,Int,String,String), it)

          ne::Int = pop(it)
          expirations = collect(String, take(it, ne))

          ns::Int = pop(it)
          strikes = collect(Float64, take(it, ns))

          w.securityDefinitionOptionalParameter(args..., expirations, strikes)
        end,

  # SECURITY_DEFINITION_OPTION_PARAMETER_END
  76 => (it, w, ver) -> w.securityDefinitionOptionalParameterEnd(slurp(Int, it)),

  # SOFT_DOLLAR_TIERS
  77 => function(it, w, ver)

          reqId::Int,
          n::Int = it

          w.softDollarTiers(reqId, [slurp(SoftDollarTier, it) for _ ∈ 1:n])
        end,

  # FAMILY_CODES
  78 => function(it, w, ver)

          n::Int = pop(it)

          w.familyCodes([slurp(FamilyCode, it) for _ ∈ 1:n])
        end,

  # SYMBOL_SAMPLES
  79 => function(it, w, ver)

          reqId::Int,
          n::Int = it

          cd = map(1:n) do _

                  c = Contract()

                  slurp!(c, (:conId,
                             :symbol,
                             :secType,
                             :primaryExchange,
                             :currency), it)

                  nd::Int = pop(it)

                  ContractDescription(c, collect(String, take(it, nd)))
                end

          w.symbolSamples(reqId, cd)
        end,

  # MKT_DEPTH_EXCHANGES
  80 => function(it, w, ver)

          n::Int = pop(it)

          df = fill_df([String,String,String,String,Union{Int,Nothing}],
                       [:exchange, :secType, :listingExch, :serviceDataType, :aggGroup],
                       n, it)

          w.mktDepthExchanges(df)
        end,

  # TICK_REQ_PARAMS
  81 => (it, w, ver) -> w.tickReqParams(slurp((Int,Float64,String,Int), it)...),

  # SMART_COMPONENTS
  82 => function(it, w, ver)

          reqId::Int,
          n::Int = it

          df = fill_df([Int,String,String], [:bit, :exchange, :exchangeLetter], n, it)

          w.smartComponents(reqId, df)
        end,

  # NEWS_ARTICLE
  83 => (it, w, ver) -> w.newsArticle(slurp((Int,Int,String), it)...),

  # TICK_NEWS
  84 => function (it, w, ver)

          args = slurp((Int,Int,String,String,String,String), it)

          # TODO: Convert timeStamp to [Zoned]DateTime
          w.tickNews(args...)
        end,

  # NEWS_PROVIDERS
  85 => function(it, w, ver)

          n::Int = pop(it)

          df = fill_df([String,String], [:providerCode, :providerName], n, it)

          w.newsProviders(df)
        end,

  # HISTORICAL_NEWS
  86 => (it, w, ver) -> w.historicalNews(slurp((Int,String,String,String,String), it)...),

  # HISTORICAL_NEWS_END
  87 => (it, w, ver) -> w.historicalNewsEnd(slurp((Int,Bool), it)...),

  # HEAD_TIMESTAMP
  88 => (it, w, ver) -> w.headTimestamp(slurp((Int,String), it)...),

  # HISTOGRAM_DATA
  89 => function(it, w, ver)

          reqId::Int,
          n::Int = it

          df = fill_df([Float64,Float64], [:price, :size], n, it)

          w.histogramData(reqId, df)
        end,

  # HISTORICAL_DATA_UPDATE
  90 => function(it, w, ver)

          reqId::Int = pop(it)

          bar = NamedTuple{(:count,
                            :time,
                            :open,
                            :close,
                            :high,
                            :low,
                            :wap,
                            :volume)}(slurp((Int,String,Float64,Float64,Float64,Float64,Float64,Float64), it))

          w.historicalDataUpdate(reqId, bar)
        end,

  # REROUTE_MKT_DATA_REQ
  91 => (it, w, ver) -> w.rerouteMktDataReq(slurp((Int,Int,String), it)...),

  # REROUTE_MKT_DEPTH_REQ
  92 => (it, w, ver) -> w.rerouteMktDepthReq(slurp((Int,Int,String), it)...),

  # MARKET_RULE
  93 => function(it, w, ver)

          marketRuleId::Int,
          n::Int = it

          df = fill_df([Float64,Float64], [:lowEdge, :increment], n, it)

          w.marketRule(marketRuleId, df)
        end,

  # PNL
  94 => function(it, w, ver)

          reqId::Int,
          dailyPnL::Float64,
          unrealizedPnL::Float64,
          realizedPnL::Float64 = it

          w.pnl(reqId, dailyPnL, unrealizedPnL, realizedPnL)
        end,

  # PNL_SINGLE
  95 => function(it, w, ver)

          reqId::Int,
          pos::Int,
          dailyPnL::Float64,
          unrealizedPnL::Union{Float64,Nothing},
          realizedPnL::Union{Float64,Nothing},
          value::Float64 = it

          w.pnlSingle(reqId, pos, dailyPnL, unrealizedPnL, realizedPnL, value)

        end,

  # HISTORICAL_TICKS
  96 => function(it, w, ver)

          reqId::Int,
          n::Int = it

          df = fill_df([Int,Int,Float64,Float64], [:time, :ignore, :price, :size], n, it)

          select!(df, Not(:ignore))

          # TODO: Convert df.time to [Zoned]DateTime

          done::Bool = pop(it)

          w.historicalTicks(reqId, df, done)
        end,

  # HISTORICAL_TICKS_BID_ASK
  97 => function(it, w, ver)

          reqId::Int,
          n::Int = it

          df = fill_df([Int,Int,Float64,Float64,Float64,Float64], [:time, :mask, :priceBid, :priceAsk, :sizeBid, :sizeAsk], n, it)

          # TODO: Convert df.time to [Zoned]DateTime
          # TODO: Unmask df.mask

          done::Bool = pop(it)

          w.historicalTicksBidAsk(reqId, df, done)
        end,

  # HISTORICAL_TICKS_LAST
  98 => function(it, w, ver)

          reqId::Int,
          n::Int = it

          df = fill_df([Int,Int,Float64,Float64,String,String], [:time, :mask, :price, :size, :exchange, :specialConditions], n, it)

          # TODO: Convert df.time to [Zoned]DateTime
          # TODO: Unmask df.mask

          done::Bool = pop(it)

          w.historicalTicksLast(reqId, df, done)
        end,

  # TICK_BY_TICK
  99 => function(it, w, ver)

          reqId::Int,
          ticktype::Int,
          time::Int = it

          local mask::Int  # To avoid "multiple type declarations" error

          if ticktype ∈ (1, 2)

            price::Float64,
            size::Float64,
            mask,
            exchange::String,
            specialConditions::String = it

            w.tickByTickAllLast(reqId, ticktype, time, price, size, unmask(TickAttribLast, mask), exchange, specialConditions)

          elseif ticktype == 3

            bidPrice::Float64,
            askPrice::Float64,
            bidSize::Float64,
            askSize::Float64,
            mask = it

            w.tickByTickBidAsk(reqId, time, bidPrice, askPrice, bidSize, askSize, unmask(TickAttribBidAsk, mask))

          elseif ticktype == 4

            w.tickByTickMidPoint(reqId, time, slurp(Float64, it))

          else
            @warn "TICK_BY_TICK: unknown ticktype" T=ticktype
          end
        end,

  # ORDER_BOUND
 100 => (it, w, ver) -> w.orderBound(collect(Int, take(it, 3))...),

  # COMPLETED_ORDER
 101 => function(it, w, ver)

          o = Order()
          c = Contract()

          slurp!(c, [1:8; 10:12], it)

          slurp!(o, 4:9, it)  # :action through :tif

          slurp!(o, (:ocaGroup,
                     :account,
                     :openClose,
                     :origin,
                     :orderRef,
                     :permId,
                     :outsideRth,
                     :hidden,
                     :discretionaryAmt,
                     :goodAfterTime,
                     :faGroup,
                     :faMethod,
                     :faPercentage,
                     :faProfile,
                     :modelCode,
                     :goodTillDate,
                     :rule80A,
                     :percentOffset,
                     :settlingFirm,
                     :shortSaleSlot,
                     :designatedLocation,
                     :exemptCode), it)

          slurp!(o, 44:48, it)    # :startingPrice through :stockRangeUpper

          slurp!(o, (:displaySize,
                     :sweepToFill,
                     :allOrNone,
                     :minQty,
                     :ocaType,
                     :triggerMethod), it)

          slurp!(o, 51:54, it)    # :volatility through :deltaNeutralAuxPrice

          !isempty(o.deltaNeutralOrderType) && slurp!(o, [55; 60:62], it)  # :deltaNeutralConId through :deltaNeutralDesignatedLocation

          slurp!(o, (:continuousUpdate,
                     :referencePriceType,
                     :trailStopPrice,
                     :trailingPercent), it)

          c.comboLegsDescrip = pop(it)

          # ComboLegs
          n::Int = pop(it)

          for _ ∈ 1:n
            push!(c.comboLegs, slurp(ComboLeg, it))
          end

          # OrderComboLeg
          n = pop(it)
          append!(o.orderComboLegs, take(it, n))

          # SmartComboRouting
          n = pop(it)
          n > 0 && (o.smartComboRoutingParams = tagvalue2nt(take(it, 2n)))

          slurp!(o, (:scaleInitLevelSize,
                     :scaleSubsLevelSize,
                     :scalePriceIncrement), it)

          !isnothing(o.scalePriceIncrement) &&
          o.scalePriceIncrement > 0         && slurp!(o, 70:76, it) # scalePriceAdjustValue through scaleRandomPercent

          o.hedgeType = pop(it)

          !isempty(o.hedgeType) && (o.hedgeParam = pop(it))

          slurp!(o, (:clearingAccount,
                     :clearingIntent,
                     :notHeld), it)

          # DeltaNeutralContract
          slurp(Bool, it) && (c.deltaNeutralContract = slurp(DeltaNeutralContract, it))

          # AlgoStrategy
          o.algoStrategy = pop(it)

          if !isempty(o.algoStrategy)
            n = pop(it)
            n > 0 && (o.algoParams = tagvalue2nt(take(it, 2n)))
          end

          o.solicited,
          ostatus::String,     # OrderState.status
          o.randomizeSize,
          o.randomizePrice = it

          o.orderType == "PEG BENCH" && slurp!(o, (:referenceContractId,
                                                   :isPeggedChangeAmountDecrease,
                                                   :peggedChangeAmount,
                                                   :referenceChangeAmount,
                                                   :referenceExchangeId), it)

          # Conditions
          n = pop(it)

          if n > 0
            for _ ∈ 1:n
              push!(o.conditions, slurp(condition_map[slurp(ConditionType, it)], it))
            end

            o.conditionsIgnoreRth,
            o.conditionsCancelOrder = it
          end

          slurp!(o, (:trailStopPrice,
                     :lmtPriceOffset,
                     :cashQty,
                     :dontUseAutoPriceForHedge,
                     :isOmsContainer), it)

          slurp!(o, 119:126, it)    # :autoCancelDate through :parentPermId

          os = OrderState(ostatus, fill(ns, 9)..., fill(nothing, 3)..., ns, ns, take(it, 2)...)

          w.completedOrder(c, o, os)
         end,

  # COMPLETED_ORDERS_END
 102 => (it, w, ver) -> w.completedOrdersEnd(),

  # REPLACE_FA_END
 103 => (it, w, ver) -> w.replaceFAEnd(slurp((Int,String), it)...),

 # WSH_META_DATA
 104 => (it, w, ver) -> w.wshMetaData(slurp((Int,String), it)...),

 # WSH_EVENT_DATA
 105 => (it, w, ver) -> w.wshEventData(slurp((Int,String), it)...)

)


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff,
    BenchmarkTools

## Observed Data
three_path_dat = Feather.read("test/comparisons/three_path_dat.feather")
three_path_par = Feather.read("test/comparisons/three_path_par.feather")
three_path_start = Feather.read("test/comparisons/three_path_start.feather")

semobserved = SemObsCommon(data = Matrix(three_path_dat))

f_tol = 2.86e-13
diff_fin = SemFiniteDiff(BFGS(), Optim.Options(;f_tol = f_tol))

## Model definition
@ModelingToolkit.variables x[1:31]

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]

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)
    )

loss = Loss([SemML(semobserved, [0.0], similar(start_val))])

imply = ImplySymbolic(A, S, F, x, start_val)

model_fin = Sem(semobserved, imply, loss, diff_fin)

#solution_fin = sem_fit(model_fin)

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

#= all(
    abs.(solution_fin.minimizer .- three_path_par.est[par_order]
        ) .< 0.05*abs.(three_path_par.est[par_order])) =#

#@benchmark sem_fit(model_fin)

using NLopt, FiniteDiff

function (model::Sem{A, B, C, D} where {A, B, C, D <: SemFiniteDiff})(par::Vector, grad::Vector)
    if length(grad) > 0
        FiniteDiff.finite_difference_gradient!(grad, model, par)
    end
    return model(par)
end

opt = Opt(:LD_LBFGS, 31)

#cache = FiniteDiff.GradientCache(start)

opt.min_objective = (x,y) -> model_fin(x,y)

opt.ftol_rel = f_tol

(optf,optx,ret) = NLopt.optimize(opt, start_val)

all(
    abs.(optx .- three_path_par.est[par_order]
        ) .< 0.05*abs.(three_path_par.est[par_order]))


text:
#
# Edit `Port for Communication with the Julia Process` to 9999
#
using Atom
using Juno; Juno.connect("host.docker.internal", 9999)


text:
using NOAu111
using PyCall
using NQCBase
using NQCModels
using NonadiabaticMolecularDynamics
using UnitfulAtomic, Unitful
using CairoMakie
using LinearAlgebra

build = pyimport("ase.build")

slab = build.fcc111("Au", size=(2,4,4), vacuum=10.0)
no = build.molecule("NO")
build.add_adsorbate(slab, no, 5, "hcp")
build.add_vacuum(slab, 30)

c, atoms, R = extract_parameters_and_positions(slab)

model = NOAu(atoms.types, c, R)
slab_plane = model.H11.image.slab_plane

bond_length_ref = austrip(1.15077u"Å")
R[3,end-1] = 40
R[3,end] = 40 + bond_length_ref
const reference = ustrip.(auconvert.(u"kJ", eigvals(potential(model, R))[1]) .* Unitful.Na)

function plot_figure_a!(ax, bond, lims)
    bond_length = austrip(bond)
    zN = austrip.((collect(lims[1]:0.01:lims[2]) .+ slab_plane) .*u"Å" )
    traj = Matrix[]
    for z in zN
        positions = copy(R)
        positions[3,end-1] = z
        positions[3,end] = z + bond_length
        push!(traj, positions)
    end
    potentials = potential.(model, traj)
    x = ustrip.(auconvert.(u"Å", zN)) .- slab_plane
    transform_and_plot_potentials!(ax, potentials, lims, x)
    return ax
end

function plot_figure_e!(ax, lims)
    z = austrip((1.6 + slab_plane)u"Å")
    rs = austrip.(collect(lims[1]:0.01:lims[2]).*u"Å")
    traj = Matrix[]
    positions = copy(R)
    positions[3, end-1] = z
    for r in rs
        p = copy(positions)
        p[3,end] = z + r
        push!(traj, p)
    end
    potentials = potential.(model, traj)
    x = ustrip.(auconvert.(u"Å", rs))
    transform_and_plot_potentials!(ax, potentials, lims, x)
    return ax
end

function transform_and_plot_potentials!(ax, potentials, lims, x)
    potentials = [ustrip.(auconvert.(u"kJ", V) .* Unitful.Na) for V in potentials]
    ground_state = [eigvals(V)[1] for V in potentials] .- reference
    v11 = [V[1,1] for V in potentials] .- reference
    v12 = [V[1,2] for V in potentials]
    v22 = [V[2,2] for V in potentials] .- reference
    lines!(ax, x, v11, label="V11", color=:green, linestyle=:dash)
    lines!(ax, x, v22, label="V22", color=:red, linestyle=:dashdotdot)
    lines!(ax, x, -v12, label="V12", color=:purple, linestyle=:dashdot)
    lines!(ax, x, ground_state, label="E0", color=:blue)
    xlims!(ax, lims[1], lims[2])
end

theme = Theme(linewidth=3)
set_theme!(theme)
f = Figure(resolution=(700, 1200), backgroundcolor = RGBf0(0.98, 0.98, 0.98))

ax1 = f[1,1] = Axis(f, yticks=LinearTicks(5))
label_a = f[1,1,TopRight()] = Label(f, "A", textsize=24)
label_a.padding = (-50, 0, -50, 0)

ax2 = f[2,1] = Axis(f)
label_b = f[2,1,TopRight()] = Label(f, "C", textsize=24)
label_b.padding = (-50, 0, -50, 0)

ax3 = f[3,1] = Axis(f)
label_c = f[3,1,TopRight()] = Label(f, "E", textsize=24)
label_c.padding = (-50, 0, -50, 0)

plot_figure_a!(ax1, 1.191u"Å", (1.0, 2.5))
ylims!(ax1, -50, 400)

plot_figure_a!(ax2, 1.6u"Å", (0.5, 3.0))
ylims!(ax2, -20, 800)

plot_figure_e!(ax3, (0.8, 2.0))
ylims!(ax3, -100, 2000)

ax1.xlabel = "N atom height / Å"
ax2.ylabel = "Energy /kJ mol⁻¹"
ax2.xlabel = "N atom height / Å"
ax3.xlabel = "N-O bond length / Å"

CairoMakie.save("plot/fig2.png", f)
f


text:
module EvoTrees

export init_evotree, grow_evotree!, grow_tree, fit_evotree, predict,
    EvoTreeRegressor, EvoTreeCount, EvoTreeClassifier, EvoTreeGaussian,
    EvoTreeRModels, importance, Random

using Base.Threads: @threads
using Statistics
using StatsBase: sample, sample!, quantile
using Random
using Distributions
using StaticArrays
using CategoricalArrays
using CUDA
using BSON
using NetworkLayout
using RecipesBase
import MLJModelInterface
import MLJModelInterface: fit, update, predict, schema
import Base: convert

include("models.jl")
include("structs.jl")
include("loss.jl")
include("eval.jl")
include("predict.jl")
include("find_split.jl")
include("fit.jl")
include("importance.jl")
include("plot.jl")
include("MLJ.jl")

include("gpu/structs_gpu.jl")
include("gpu/loss_gpu.jl")
include("gpu/eval_gpu.jl")
include("gpu/predict_gpu.jl")
include("gpu/find_split_gpu.jl")
include("gpu/fit_gpu.jl")

function convert(::Type{GBTree}, m::GBTreeGPU)
    EvoTrees.GBTree([EvoTrees.Tree(Array(tree.feat),
            Array(tree.cond_bin),
            Array(tree.cond_float),
            Array(tree.gain),
            Array(tree.pred),
            Array(tree.split)) for tree in m.trees],
        m.params,
        m.metric,
        m.K,
        m.levels)
end

function save(model::GBTree, path)
    BSON.bson(path, Dict(:model => model))
end

function save(model::GBTreeGPU, path)
    m = convert(GBTree, model)
    save(m, path)
end

function load(path)
    m = BSON.load(path, @__MODULE__)
    return m[:model]
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ aeae4550-720d-11eb-3ef5-1f23537ec599
using GMT

# ╔═╡ 0745ac40-720d-11eb-1acd-2d4920b5c475
md"# Contourf examples

GMT does not actually have a *contourf* module like Matlab for example, but we can obtain the same result using *grdview*, *grdcontour* and *pscontour*. However, to make things the Julia wrapper wrapped up a module called *contourf* that makes it really easy to use. To show how it works let's start by creating an example grid and a CPT.
"

# ╔═╡ b2cdb8a0-720d-11eb-200b-cd2c056c8375
G = GMT.peaks();

# ╔═╡ d747f28e-720d-11eb-0e68-a7242445823b
C = makecpt(T=(-7,9,2));

# ╔═╡ df241fc0-720d-11eb-13d2-c311aa8f5518
md"Now if we pass those two to the *contourf* module we get an annotated plot where the annotations come from the color CPT.
"

# ╔═╡ 0faaea20-720e-11eb-0ff8-1105205be52c
contourf(G, C, show=true)

# ╔═╡ 29a06c20-720e-11eb-312a-ff244a4d7a71
md"If we want to just draw some contours and not annotate them, we pass an array with the contours to be drawn.
"

# ╔═╡ 5ca0c200-720e-11eb-3a76-5db53780ef22
contourf(G, C, contour=[-2, 0, 2, 5], show=true)

# ╔═╡ 6a4f90c0-720e-11eb-29c7-1529f5116d73
md"### What if one has an *x,y,z* file instead of a grid?
That is also simple, let's simulate it with synthetic data.
"

# ╔═╡ a9d59eae-720e-11eb-050a-d1737498f9a0
d = [0 2 5; 1 4 5; 2 0.5 5; 3 3 9; 4 4.5 5; 4.2 1.2 5; 6 3 1; 8 1 5; 9 4.5 5];

# ╔═╡ d12ce0e0-720e-11eb-05c4-3b65ee0d2da2
contourf(d, limits=(-0.5,9.5,0,5), pen=0.25, labels=(line=(:min,:max),), show=1)

# ╔═╡ b22ad6fe-720f-11eb-2de4-37064019d3f6
md"In the above since we did not specify a CPT the program picked the GMT's default one. But if we want use another one it's only a question of creating and passed it in.
"

# ╔═╡ f583bede-720f-11eb-0ff0-617f73c13447
cpt = makecpt(range=(0,10,1), cmap=:batlow);

# ╔═╡ feda3000-720f-11eb-2b2f-d546d65f4abe
contourf(d, contours=cpt, limits=(-0.5,9.5,0,5), pen=0.25, labels=(line=(:min,:max),), show=true)

# ╔═╡ Cell order:
# ╟─0745ac40-720d-11eb-1acd-2d4920b5c475
# ╠═aeae4550-720d-11eb-3ef5-1f23537ec599
# ╠═b2cdb8a0-720d-11eb-200b-cd2c056c8375
# ╠═d747f28e-720d-11eb-0e68-a7242445823b
# ╟─df241fc0-720d-11eb-13d2-c311aa8f5518
# ╠═0faaea20-720e-11eb-0ff8-1105205be52c
# ╠═29a06c20-720e-11eb-312a-ff244a4d7a71
# ╠═5ca0c200-720e-11eb-3a76-5db53780ef22
# ╟─6a4f90c0-720e-11eb-29c7-1529f5116d73
# ╠═a9d59eae-720e-11eb-050a-d1737498f9a0
# ╠═d12ce0e0-720e-11eb-05c4-3b65ee0d2da2
# ╟─b22ad6fe-720f-11eb-2de4-37064019d3f6
# ╠═f583bede-720f-11eb-0ff0-617f73c13447
# ╠═feda3000-720f-11eb-2b2f-d546d65f4abe


text:
Base.show(io::IO, ::MIME"text/latex", x::Symbolic) = print(io, latexify(x))

@latexrecipe function f(op::Symbolic{Number})
    expr = Expr(:latexifymerge, op.name)
    return expr
end

@latexrecipe function f(op::Add{Number})
    expr = Expr(:latexifymerge)
    expr.args = _to_expression(op)
    return expr
end

@latexrecipe function f(op::Mul{Number})
    expr = Expr(:latexifymerge)
    expr.args = _to_expression(op)
    return expr
end

@latexrecipe function f(op::Div{Number})
    expr = Expr(:latexifymerge)
    expr.args = _to_expression(op)
    return expr
end

@latexrecipe function f(op::Pow{Number})
    expr = Expr(:latexifymerge)
    expr.args = _to_expression(op)
    return expr
end

@latexrecipe function f(op::Symbolic{QOperator})
    expr = Expr(:latexifymerge)
    expr.args = _to_expression(op)
    return expr
end

@latexrecipe function f(my_term::Term{QOperator})
    my_expr = Expr(:latexifymerge)
    my_expr.args = _to_expression(my_term)
    return my_expr
end

_to_expression(x::Number) = x
_to_expression(op::Symbolic{Number}) = String(op.name)

function _to_expression(op::Symbolic{QOperator})
    metadt = getmetadata(op, QOperatorMeta)
    op_type = metadt.type
    if op_type == BosonicCreate()
        return ["\\hat{", String(op.name)[1:end-2], "}^\\dagger"]
    elseif op_type == PauliSM()
        return ["\\hat{", op.name, "}_{-}"]
    elseif op_type == PauliSP()
        return ["\\hat{", String(op.name)[1:end-2], "}_{+}"]
    else
        return ["\\hat{", op.name, "}"]
    end
end

function _to_expression(my_term::Term{QOperator})
    args = arguments(my_term)
    my_args = []
    my_f = operation(my_term)
    if my_f == (*)
        for op in args
            if isterm_op(+)(op) || isterm_op(-)(op)
                append!(my_args, "\\left(")
                append!(my_args, _to_expression(op))
                append!(my_args, "\\right)")
            else
                if !isequal(op, 1)
                    if isequal(op, -1)
                        append!(my_args, "-")
                    else
                        append!(my_args, _to_expression(op))
                    end
                end
            end
        end
    elseif my_f == (/) && length(args) == 2
        if istree(args[1])
            append!(my_args, ["\\frac{1}{", args[2], "}"])
            append!(my_args, "\\left(")
            append!(my_args, _to_expression(args[1]))
            append!(my_args, "\\right)")
        else
            append!(my_args, "\\frac{")
            append!(my_args, _to_expression(args[1]))
            append!(my_args, "}{")
            append!(my_args, _to_expression(args[2]))
            append!(my_args, "}")
        end
    elseif my_f == (^) && length(args) == 2
        if istree(args[1])
            append!(my_args, "\\left(")
            append!(my_args, _to_expression(args[1]))
            append!(my_args, "\\right)^{")
            append!(my_args, _to_expression(args[2]))
            append!(my_args, "}")
        else
            append!(my_args, "{")
            append!(my_args, _to_expression(args[1]))
            append!(my_args, "}^{")
            append!(my_args, _to_expression(args[2]))
            append!(my_args, "}")
        end
    elseif my_f == (+)
        append!(my_args, _to_expression(args[1]))
        for op in args[2:end]
            if istree(op)
                if isterm_op(*)(op)
                    if isequal(arguments(op)[1],-1)
                        append!(my_args, "-")
                        append!(my_args, _to_expression(*(arguments(op)[2:end]...)))
                    else
                        append!(my_args, "+")
                        append!(my_args, _to_expression(op))
                    end
                else
                    append!(my_args, "+")
                    append!(my_args, _to_expression(op))
                end
            else
                if isequal(op, -1)
                    append!(my_args, "-1")
                else
                    append!(my_args, "+")
                    append!(my_args, _to_expression(op))
                end
            end
        end
    elseif my_f == (sin)
        append!(my_args, "\\sin \\left(")
        append!(my_args, _to_expression(args[1]))
        append!(my_args, "\\right)")
    elseif my_f == (cos)
        append!(my_args, "\\cos \\left(")
        append!(my_args, _to_expression(args[1]))
        append!(my_args, "\\right)")
    elseif my_f == (exp)
        append!(my_args, "e^{")
        append!(my_args, _to_expression(args[1]))
        append!(my_args, "}")
    end

    return my_args
end

function _to_expression(my_term::Add{Number})
    args = SymbolicUtils.arguments(my_term)
    my_args = []
    append!(my_args, _to_expression(args[1]))
    for op in args[2:end]
        if isterm_op(*)(op)
            if arguments(op)[1] == -1
                append!(my_args, "-")
                append!(my_args, _to_expression(*(arguments(op)[2:end]...)))
            else
                append!(my_args, "+")
                append!(my_args, _to_expression(op))
            end
        else
            append!(my_args, "+")
            append!(my_args, _to_expression(op))
        end
    end
    return my_args
end

function _to_expression(my_term::Mul{Number})
    args = SymbolicUtils.arguments(my_term)
    my_args = []
    for op in args
        if isterm_op(+)(op) || isterm_op(-)(op)
            append!(my_args, "\\left(")
            append!(my_args, _to_expression(op))
            append!(my_args, "\\right)")
        else
            append!(my_args, _to_expression(op))
        end
    end
    return my_args
end

function _to_expression(my_term::Div{Number})
    args = SymbolicUtils.arguments(my_term)
    my_args = []
    if length(args) == 2
        append!(my_args, ["\\frac{", args[1], "}{", args[2], "}"])
    end
    return my_args
end

function _to_expression(my_term::Pow{Number})
    args = SymbolicUtils.arguments(my_term)
    my_args = []
    if istree(args[1])
        append!(my_args, "\\left(")
        append!(my_args, _to_expression(args[1]))
        append!(my_args, "\\right)^{")
        append!(my_args, _to_expression(args[2]))
        append!(my_args, "}")
    else
        append!(my_args, "{")
        append!(my_args, _to_expression(args[1]))
        append!(my_args, "}^{")
        append!(my_args, _to_expression(args[2]))
        append!(my_args, "}")
    end
    return my_args
end

# function _to_expression(op::Symbolic{QOperator})
#     if !istree(op)
#         metadt = getmetadata(op, QOperatorMeta)
#         op_type = metadt.type
#         if op_type == BosonicCreate()
#             return ["\\hat{", String(op.name)[1:end-2], "}^\\dagger"]
#         elseif op_type == PauliSM()
#             return ["\\hat{", op.name, "}_{-}"]
#         elseif op_type == PauliSP()
#             return ["\\hat{", String(op.name)[1:end-2], "}_{+}"]
#         else
#             return ["\\hat{", op.name, "}"]
#         end
#     else
#         args = arguments(op)
#         my_args = []
#         if ismul(op)
#             for op_tmp in args
#                 if isadd(op_tmp)
#                     append!(my_args, "\\left(")
#                     append!(my_args, _to_expression(op_tmp))
#                     append!(my_args, "\\right)")
#                 else
#                     append!(my_args, _to_expression(op_tmp))
#                 end
#             end
#         elseif isdiv(op)
#             if istree(args[1])
#                 append!(my_args, ["\\frac{1}{", args[2], "}"])
#                 append!(my_args, "\\left(")
#                 append!(my_args, _to_expression(args[1]))
#                 append!(my_args, "\\right)")
#             else
#                 append!(my_args, "\\frac{")
#                 append!(my_args, _to_expression(args[1]))
#                 append!(my_args, "}{")
#                 append!(my_args, _to_expression(args[2]))
#                 append!(my_args, "}")
#             end
#         elseif ispow(op)
#             if istree(args[1])
#                 append!(my_args, "\\left(")
#                 append!(my_args, _to_expression(args[1]))
#                 append!(my_args, "\\right)^{")
#                 append!(my_args, _to_expression(args[2]))
#                 append!(my_args, "}")
#             else
#                 append!(my_args, "{")
#                 append!(my_args, _to_expression(args[1]))
#                 append!(my_args, "}^{")
#                 append!(my_args, _to_expression(args[2]))
#                 append!(my_args, "}")
#             end
#         elseif isadd(op)
#             append!(my_args, _to_expression(args[1]))
#             for op_tmp in args[2:end]
#                 if istree(op_tmp)
#                     if ismul(op_tmp)
#                         if arguments(op_tmp)[1] == -1
#                             append!(my_args, "-")
#                             append!(my_args, _to_expression(*(arguments(op_tmp)[2:end]...)))
#                         else
#                             append!(my_args, "+")
#                             append!(my_args, _to_expression(op_tmp))
#                         end
#                     else
#                         append!(my_args, "+")
#                         append!(my_args, _to_expression(op_tmp))
#                     end
#                 else
#                     if op_tmp == -1
#                         append!(my_args, "-1")
#                     else
#                         append!(my_args, "+")
#                         append!(my_args, _to_expression(op_tmp))
#                     end
#                 end
#             end
#         end
#         return my_args
#     end
# end


text:
using Test
using GraphSearchZero

macro testblock(ex)
    quote
        try
            $(esc(eval(ex)))
            true
        catch err
            isa(err, ErrorException) ? false : rethrow(err)
        end
    end
end

include("test_search.jl")


text:
module GraphCircuits
using LightGraphs

export find_circuits

include("johnson.jl")

include("hadwick_james.jl")

"""

"""
find_circuits(g::DiGraph) = findcircuits(g, Val{:johnson})

end # module


text:
"""
Trait for `ODEOperator` that tells us whether the operator depends on the solution
(including its time derivatives), it is an affine operator that depends on time
or it is a constant operator (affine and time-indepedendent)
"""
abstract type OperatorType end
struct Nonlinear <: OperatorType end
struct Affine  <: OperatorType end
struct Constant  <: OperatorType end
struct ConstantMatrix  <: OperatorType end

"""
It represents the operator in an implicit N-th order ODE, i.e., A(t,u,∂tu,∂t^2u,...,∂t^Nu)
where the implicit PDE reads A(t,u,∂tu,∂t^2u,...,∂t^Nu) = 0, when ∂t^iu is the
i-th time derivative of u, with i=0,..,N. The trait `{C}` determines whether the
operator is fully nonlinear, affine or constant in time.
"""
abstract type ODEOperator{C<:OperatorType} <: GridapType end

"""
It represents an _affine_ operator in an implicit ODE, i.e., an ODE operator of
the form A(t,u,∂tu,...,∂t^Nu) = A_N(t)∂t^Nu + ...A_1(t)∂tu + A_0(t)u + f(t)
"""
const AffineODEOperator = ODEOperator{Affine}

"""
It represents a constant operator in an implicit ODE, i.e., an ODE operator of
the form A(t,u,∂tu,...,∂t^Nu) = A_N∂t^Nu + ...A_1∂tu + A_0u + f
"""
const ConstantODEOperator = ODEOperator{Constant}

"""
It represents an affine operator in an implicit ODE with constant matrix, but
time-dependent right-hand side, i.e., an ODE operator of
the form A(t,u,∂tu,...,∂t^Nu) = A_N∂t^Nu + ...A_1∂tu + A_0u + f(t)
"""
const ConstantMatrixODEOperator = ODEOperator{ConstantMatrix}

"""
Returns the `OperatorType`, i.e., nonlinear, affine, or constant in time
"""
OperatorType(::ODEOperator{C}) where {C} = C

"""
Returns the order of the ODE operator
"""
function get_order(::ODEOperator)
  @abstractmethod
end

"""
It provides A(t,u,∂tu,...,∂t^Nu) for a given (t,u,∂tu,...,∂t^Nu)
"""
function residual!(
  r::AbstractVector,
  op::ODEOperator,
  t::Real,
  u::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  ode_cache)
  @abstractmethod
end

"""
"""
function allocate_residual(
  op::ODEOperator,
  u::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  ode_cache)
  @abstractmethod
end

"""
It adds contribution to the Jacobian with respect to the i-th time derivative,
with i=0,...,N. That is, adding γ_i*[∂A/∂(∂t^iu)](t,u,∂tu,...,∂t^Nu) for a
given (t,u,∂tu,...,∂t^Nu) to a given matrix J, where γ_i is a scaling coefficient
provided by the `ODESolver`, e.g., 1/Δt for Backward Euler; It represents
∂(δt^i(u))/∂(u), in which δt^i(⋅) is the approximation of ∂t^i(⋅) in the solver.
Note that for i=0, γ_i=1.0.
"""
function jacobian!(
  J::AbstractMatrix,
  op::ODEOperator,
  t::Real,
  u::Tuple{Vararg{AbstractVector}},
  i::Integer,
  γᵢ::Real,
  ode_cache)
  @abstractmethod
  # Add values to J
end

"""
Add the contribution of all jacobians ,i.e., ∑ᵢ γ_i*[∂A/∂(∂t^iu)](t,u,∂tu,...,∂t^Nu)
"""
function jacobians!(
  J::AbstractMatrix,
  op::ODEOperator,
  t::Real,
  u::Tuple{Vararg{AbstractVector}},
  γ::Tuple{Vararg{Real}},
  ode_cache)
  @abstractmethod
  # Add values to J
end

"""
"""
function allocate_jacobian(op::ODEOperator,u::AbstractVector,ode_cache)
  @abstractmethod
end

"""
Allocates the cache data required by the `ODESolution` for a given `ODEOperator`
"""
allocate_cache(op::ODEOperator) = @abstractmethod

#@fverdugo to be used as `cache = update_cache!(cache,op,t)`
update_cache!(cache,op::ODEOperator,t::Real) = @abstractmethod

"""
Tests the interface of `ODEOperator` specializations
"""
function test_ode_operator(op::ODEOperator,t::Real,u::AbstractVector,u_t::AbstractVector)
  cache = allocate_cache(op)
  cache = update_cache!(cache,op,0.0)
  r = allocate_residual(op,u,cache)
  residual!(r,op,t,(u,u_t),cache)
  J = allocate_jacobian(op,u,cache)
  jacobian!(J,op,t,(u,u_t),1,1.0,cache)
  jacobian!(J,op,t,(u,u_t),2,1.0,cache)
  jacobians!(J,op,t,(u,u_t),(1.0,1.0),cache)
  true
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

function prepareValueReference(md::fmi3ModelDescription, vr::fmi3ValueReferenceFormat)
    tvr = typeof(vr)
    if isa(vr, AbstractArray{fmi3ValueReference,1})
        return vr
    elseif tvr == fmi3ValueReference
        return [vr]
    elseif tvr == String
        return [fmi3StringToValueReference(md, vr)]
    elseif isa(vr, AbstractArray{String,1})
        return fmi3StringToValueReference(md, vr)
    elseif tvr == Int64
        return [fmi3ValueReference(vr)]
    elseif isa(vr, AbstractArray{Int64,1})
        return fmi3ValueReference.(vr)
    elseif tvr == Nothing
        return Array{fmi3ValueReference,1}()
    end

    @assert false "prepareValueReference(...): Unknown value reference structure `$tvr`."
end
function prepareValueReference(fmu::FMU3, vr::fmi3ValueReferenceFormat)
    prepareValueReference(fmu.modelDescription, vr)
end
function prepareValueReference(comp::FMU3Instance, vr::fmi3ValueReferenceFormat)
    prepareValueReference(comp.fmu.modelDescription, vr)
end

"""
Returns an array of ValueReferences coresponding to the variable names.
"""
function fmi3StringToValueReference(md::fmi3ModelDescription, names::AbstractArray{String})
    vr = Array{fmi3ValueReference}(undef,0)
    for name in names
        reference = fmi3StringToValueReference(md, name)
        if reference === nothing
            @warn "Value reference for variable '$name' not found, skipping."
        else
            push!(vr, reference)
        end
    end
    vr
end

"""
Returns the model variable(s) fitting the value reference.
"""
function fmi3ModelVariablesForValueReference(md::fmi3ModelDescription, vr::fmi3ValueReference)
    ar = []
    for modelVariable in md.modelVariables
        if modelVariable.valueReference == vr
            push!(ar, modelVariable)
        end
    end
    ar
end

"""
Returns the ValueReference coresponding to the variable name.
"""
function fmi3StringToValueReference(md::fmi3ModelDescription, name::String)
    reference = nothing
    if haskey(md.stringValueReferences, name)
        reference = md.stringValueReferences[name]
    else
        @warn "No variable named '$name' found."
    end
    reference
end

function fmi3StringToValueReference(fmu::FMU3, name::Union{String, AbstractArray{String}})
    fmi3StringToValueReference(fmu.modelDescription, name)
end

"""
Returns an array of variable names matching a fmi3ValueReference.
"""
function fmi3ValueReferenceToString(md::fmi3ModelDescription, reference::fmi3ValueReference)
    [k for (k,v) in md.stringValueReferences if v == reference]
end
function fmi3ValueReferenceToString(md::fmi3ModelDescription, reference::Int64)
    fmi3ValueReferenceToString(md, fmi3ValueReference(reference))
end

function fmi3ValueReferenceToString(fmu::FMU3, reference::Union{fmi3ValueReference, Int64})
    fmi3ValueReferenceToString(fmu.modelDescription, reference)
end


text:
type MatchCapts
    data::Dict{Symbol,Any}
end
newmatchcapts() = MatchCapts(Dict{Symbol,Any}())

# NB, If we use this code, don't use tasks. Theyr
# are inefficient.

# This type corresponds to the expression:
# PatternTest[Pattern[name,Blank[head]],test]

# abstract AbstractPattern
# type BlankT <: AbstractPattern
#     head::Symbol
# end
# type PatternT1 <: AbstractPattern
#     name::Symbol
#     head::Symbol
# end
# function toBlankPat(mx::Mxpr{:Blank})
#     if length(mx) > 0
#         BlankT(mx[1])
#     else
#         BlankPat(:Any)
#     end
# end
# function toBlankPat(mx::Mxpr{:Pattern})
#     if length(mx) > 0
#         BlankT(mx[1])
#     else
#         BlankPat(:Any)
#     end
# end

# Match item to a single term in args of an AC expression.
# item can potentially match any of the terms. Find each match
# sequentially and return captures from that match and the remaining
# arguments so they can be matched by the next item.
single_ac_matches(item,args) = @task _single_ac_matches(item,args)
function _single_ac_matches(item,args)
    for i in 1:length(args)
        (goodflag,capts) = find_match(item,args[i])
        if goodflag
            restargs = copy(args)
            splice!(restargs,i)
            produce (capts, restargs)
        end
    end
end

# Patterns, such as  x_ can match multiple flat, orderess terms.
# I think, only x_, x__, x___, _, __, and ___ can match.
# Other qualifiers, heads and tests, mean only one term can match.
# This iterable (task) returns all partitions of args into two
# (a,b). First all subsets a of length 1 are chosen, then of
# length 2, etc. This is the correct order, since x_ should match
# the minimum number of terms possible.
mult_ac_matches(item,args) = @task _mult_ac_matches(item,args)
function _mult_ac_matches(item,args)
    for n in 1:length(args)
        for x in combinations(args,n) # all subsets of args of length n
            produce(x,setdiff(args,x)) # x and complement of x in args
        end
    end
end

function find_match(item,args)
    return (true, ())
end


text:
# Component Wise Weno and Mapped Weno schemes
# Flux splitting using local (LLF) or global Lax-Flux splitting (GLF).
#
# Spectral Mapped Weno
#
# Based on:
# C.-W. Shu, “High order ENO and WENO schemes for computational fluid
# dynamics,” in High-Order Methods for Computational Physics, Lecture Notes
#in Computational Science and Engineering vol. 9, 439–582. New York:
# Springer Verlag, 1999

# A. Henrick, T. Aslam, J. Powers, Mapped weighted essentially non-oscillatory
# schemes: Achiving optimal order near critical points

# R. Bürger, R. Donat, P. Mulet, C. Vega, On the implementation of WENO schemes
# for a class of polydisperse sedimentation models. October 21, 2010.

mutable struct FVCompWENOScheme{oType, aType, cType} <: AbstractFVAlgorithm
  order :: oType
  splitting :: Symbol #Splitting strategy local Lax-Friedrichs (LLF) or global
  α :: aType
  rec_scheme :: cType
  αf :: Function
end

function FVCompWENOScheme(;order=5, splitting = :GLF, αf = nothing)
  if αf == nothing
    αf = maxfluxρ
  end
  FVCompWENOScheme(order, splitting,0.0, WENO_Reconstruction(order), αf)
end

mutable struct FVCompMWENOScheme{oType, aType, cType} <: AbstractFVAlgorithm
  order :: oType
  splitting :: Symbol
  α :: aType
  rec_scheme :: cType
  αf :: Function
end

function FVCompMWENOScheme(;order=5, splitting = :GLF, αf = nothing)
  if αf == nothing
    αf = maxfluxρ
  end
  FVCompMWENOScheme(order, splitting, 0.0, MWENO_Reconstruction(order), αf)
end

struct FVSpecMWENOScheme{cType} <: AbstractFVAlgorithm
  order :: Int
  rec_scheme :: cType
end

FVCompWENOSchemes = Union{FVCompWENOScheme,FVCompMWENOScheme}

function update_dt(alg::FVCompWENOScheme{T1,T2,T3},u,Flux,
    CFL,mesh) where {T1,T2,T3}
  alg.α = alg.αf(u,Flux,mesh)
  @assert (abs(alg.α) > eps(T2))
  dx = cell_volume(mesh, 1)
  CFL*dx/alg.α
 end

 function update_dt(alg::FVCompMWENOScheme{T1,T2,T3},u,Flux,
  CFL,mesh) where {T1,T2,T3}
  alg.α = alg.αf(u,Flux,mesh)
  @assert (abs(alg.α) > eps(T2))
  dx = cell_volume(mesh, 1)
  CFL*dx/alg.α
end

"""
FVSpecMWENOScheme(;order=5)
Initialize Spectral Mapped Weno algorithm with
WENO reconstruction of `order` default WENO5
"""
function FVSpecMWENOScheme(;order=5)
  FVSpecMWENOScheme(order, MWENO_Reconstruction(order))
end
##############################################################
#Component Wise WENO Scheme and Mapped Component Wise WENO Scheme
##############################################################
function update_flux_value(j,k,fminus,fplus,mesh,alg::FVCompWENOSchemes,i)
      fm = get_cellvals(fminus,mesh,(i,j-k+1:j+k+1)...)
      fp = get_cellvals(fplus,mesh,(i,j-k:j+k)...)
      return sum(reconstruct(fm, fp, alg.rec_scheme))
end
function update_flux_value(j,k,fminus,fplus,mesh,alg::FVCompWENOSchemes)
  fm = get_cellvals(fminus,mesh,(j-k+1:j+k+1,)...)
  fp = get_cellvals(fplus,mesh,(j-k:j+k,)...)
  return sum(reconstruct(fm, fp, alg.rec_scheme))
end
"""
compute_fluxes!(fluxes, Flux, u, mesh, dt, M, alg::FVCompWENOScheme, ::Type{Val{true}})
Numerical flux of Component Wise WENO algorithm Scheme in 1D
"""
function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::FVCompWENOSchemes, nonscalar::Bool, ::Type{Val{true}})
    order = alg.order; splitting = alg.splitting
    k = Int((order + 1)/2)-1
    # Flux splitting
    if splitting == :GLF
      fminus, fplus = glf_splitting(u, alg.α, Flux, nonscalar, Val{true})
    elseif splitting == :LLF
      fminus, fplus = llf_splitting(u, mesh, Flux, nonscalar)
    else
      throw("Splitting strategy not supported...")
    end
    #Compute numerical fluxes
    Threads.@threads for j in node_indices(mesh)
      if nonscalar
        for i in 1:size(fluxes,1)
          fluxes[i,j] = update_flux_value(j-1,k,fminus,fplus,mesh,alg,i)
        end
      else
          fluxes[j] = update_flux_value(j-1,k,fminus,fplus,mesh,alg)
      end
    end
end

function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::FVCompWENOSchemes, nonscalar::Bool, ::Type{Val{false}})
  order = alg.order; splitting = alg.splitting
  k = Int((order + 1)/2)-1
  # Flux splitting
  if splitting == :GLF
    fminus, fplus = glf_splitting(u, alg.α, Flux, nonscalar, Val{false})
  elseif splitting == :LLF
    fminus, fplus = llf_splitting(u, mesh, Flux, nonscalar)
  else
    throw("Splitting strategy not supported...")
  end
  #Compute numerical fluxes
  for j in node_indices(mesh)
    if nonscalar
      for i in 1:size(fluxes,1)
        fluxes[i,j] = update_flux_value(j-1,k,fminus,fplus,mesh,alg,i)
      end
    else
        fluxes[j] = update_flux_value(j-1,k,fminus,fplus,mesh,alg)
    end
  end
end

###############################################################
#Characteristic Wise WENO algorithm (Spectral)
#################################################################
"""
compute_fluxes!(fluxes, Flux, u, mesh, dt, M, alg::FVSpecMWENOScheme, ::Type{Val{true}})
Numerical flux of Mapped Spectral WENO algorithm Scheme in 1D
"""
function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::FVSpecMWENOScheme, nonscalar::Bool, ::Type{Val{false}})
    N = getncells(mesh.mesh)
    order = alg.order; rec_scheme = alg.rec_scheme
    k = Int((order + 1)/2)-1
    if nonscalar
      M = size(u,1)
      save_case = fill(0.0,M,N+1)
      αj = fill(0.0,M,N+1)
      RMats = Vector{typeof(eigvecs(Flux.Dflux(u[:,1])))}(undef,0)
      LMats = Vector{eltype(RMats)}(undef,0)
      gk = fill!(similar(u), zero(eltype(u)))
      for j in node_indices(mesh)
        @inbounds ul = cellval_at_left(j,u,mesh)
        @inbounds ur = cellval_at_right(j,u,mesh)
        MatJf = Flux.Dflux(0.5*(ul+ur))
        Rj = eigvecs(MatJf);  Lj = inv(Rj)
        push!(RMats,Rj); push!(LMats,Lj)
        λl = eigvals(Flux.Dflux(ul)); λr = eigvals(Flux.Dflux(ur))
        αj[:,j] = maximum(abs,[λl λr],dims=2)
        if j < N+1
          gk[:,j] = Flux(u[:,j])
        end
        for i in 1:M
          if λl[i]*λr[i] <= 0
            save_case[i,j] = 1
          else
            if λl[i] > 0 && λr[i] > 0
              save_case[i,j] = 2
            else
              save_case[i,j] = 3
            end
          end
        end
      end
    #WEno Reconstrucion
    gklloc = fill(0.0,M,k*2+1);gkrloc = fill(0.0,M,k*2+1)
    gmloc = fill(0.0,M,k*2+1);gploc = fill(0.0,M,k*2+1)
    for j = 0:N
      for (ll,l) in enumerate((j-k):(j+k))
        gkl = get_cellvals(gk,mesh,(:,l+1)...)
        gkr = get_cellvals(gk,mesh,(:,l)...)
        ul = get_cellvals(u,mesh,(:,l+1)...)
        ur = get_cellvals(u,mesh,(:,l)...)
        gklloc[:,ll] = LMats[j+1]*gkl
        gkrloc[:,ll] =  LMats[j+1]*gkr
        gmloc[:,ll] = 0.5*LMats[j+1]*(gkl-αj[:,j+1].*ul)
        gploc[:,ll] = 0.5*LMats[j+1]*(gkr+αj[:,j+1].*ur)
      end
      for i = 1:M
        if save_case[i,j+1] == 1
          fluxes[i,j+1] = sum(reconstruct(gmloc[i,:],gploc[i,:], rec_scheme))
        elseif save_case[i,j+1] == 2
          fluxes[i,j+1] = reconstruct(gklloc[i,:],gkrloc[i,:],rec_scheme)[2]
        else
          fluxes[i,j+1] = reconstruct(gklloc[i,:],gkrloc[i,:],rec_scheme)[1]
        end
      end
      fluxes[:,j+1] = RMats[j+1]*fluxes[:,j+1]
    end
  else
    save_case = fill(0.0,N+1)
    αj = fill(0.0,N+1)
    gk = fill!(similar(u), zero(eltype(u)))
    for j in node_indices(mesh)
      @inbounds ul = cellval_at_left(j,u,mesh)
      @inbounds ur = cellval_at_right(j,u,mesh)
      λl = Flux.Dflux(ul); λr = Flux.Dflux(ur)
      αj[j] = maximum(abs,[λl,λr])
      if j < N+1
        gk[j] = Flux(u[j])
      end
      if λl*λr <= 0
        save_case[j] = 1
      else
        if λl > 0 && λr > 0
          save_case[j] = 2
        else
          save_case[j] = 3
        end
      end
    end
    #WEno Reconstrucion
    gklloc = fill(0.0,k*2+1);gkrloc = fill(0.0,k*2+1)
    gmloc = fill(0.0,k*2+1);gploc = fill(0.0,k*2+1)
    for j in 0:N
      for (ll,l) in enumerate((j-k):(j+k))
        gkl = get_cellvals(gk,mesh,(l+1,)...)
        gkr = get_cellvals(gk,mesh,(l,)...)
        ul = get_cellvals(u,mesh,(l+1,)...)
        ur = get_cellvals(u,mesh,(l,)...)
        gklloc[ll] = gkl
        gkrloc[ll] = gkr
        gmloc[ll] = 0.5*(gkl-αj[j+1]*ul)
        gploc[ll] = 0.5*(gkr+αj[j+1]*ur)
      end
      if save_case[j+1] == 1
        fluxes[j+1] = sum(reconstruct(gmloc[:],gploc[:], rec_scheme))
      elseif save_case[j+1] == 2
        fluxes[j+1] = reconstruct(gklloc[:],gkrloc[:],rec_scheme)[2]
      else
        fluxes[j+1] = reconstruct(gklloc[:],gkrloc[:],rec_scheme)[1]
      end

    end
  end
end


text:
include("../../../sim-study/configs/test-sim-6-8-5/tsne/tsne_imports.jl")

function get_output_dir(path_to_output)
  joinpath(splitdir(path_to_output)[1], "output.bson")
end

### Main ###

# Results directory
results_dir = "$(ENV["SCRATCH_DIR"])/cytof/results/repfam/patients-data/run11"

# Where to put results
csvdir = "img"

# Loop through different settings
for phi in (0, 1, 100, 10000)  # NOTE: mind this! (0, 1, 25, 50, 100, 10000)
  # Get simulation name
  simname = "phi$(phi)"
  println("Processing: ", simname)

  # Get path to simuldate data.
  path_to_output = joinpath(results_dir, simname, "output.bson")

  # Get output
  output = BSON.load(path_to_output)

  # Get relevant parameters
  samples =  [s[1] for s in output[:samples][1]]

  # Get best lam (for each sample)
  best_lam, N = let
    _best_lam = get_best_lam(samples)
    vcat(_best_lam...), length.(_best_lam)
  end

  # number of samples
  I = length(N)

  # Sample index
  sample_ind = vcat([fill(i, N[i]) for i in 1:I]...)

  # Write output to file.
  open("$(csvdir)/fam-$(simname)-clusterings.csv", "w") do io
    writedlm(io, [best_lam sample_ind])
  end
end


text:
using Documenter, JuliaForBioinformatics

include("literate.jl")

makedocs(;
    modules=[JuliaForBioinformatics],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "Introduction" => [
            "01_Introduction.md"
            "02_DataStructures.md"
            "03_Files_Homework.md"
        ],
        "Advance" => [
            "04_HowJuliaWorks.md"
            "05_DefiningTypes.md"
            "06_Strings.md"
            "07_Stats.md"
            "08_Packages_Homework.md"
        ],
        "Bioinformatics" => [
            "09_MSA.md"
            "10_Information.md"
            "11_PDB.md"
        ],
        "API" => "api.md"
    ],
    # repo="https://github.com/diegozea/JuliaForBioinformatics.jl/blob/{commit}{path}#L{line}",
    sitename="JuliaForBioinformatics",
    authors="Diego Javier Zea"
)

deploydocs(;
    repo="github.com/diegozea/JuliaForBioinformatics",
)

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
module CStatic

export @cstatic

# ref: https://github.com/JuliaLang/julia/issues/15056#issuecomment-183937358
macro cstatic(exprs...)
    prologue = Expr(:block)
    epilogue = Expr(:block)
    ret = length(exprs) > 2 ? Expr(:tuple) : first(exprs).args[1]
    for expr in exprs[1:end-1]
        Meta.isexpr(expr, :(=)) || throw(ArgumentError("syntax mismatch! should be `@cstatic x1=1 x2=2 ... xn=n expr`."))
        local_sym, init_expr = expr.args
        global_sym = gensym("static_$local_sym")
        Base.eval(__module__, Expr(:(=), global_sym, init_expr))
        push!(prologue.args, Expr(:global, global_sym))
        push!(prologue.args, Expr(:local, Expr(:(=), local_sym, global_sym)))
        push!(epilogue.args, Expr(:(=), global_sym, local_sym))
        length(exprs) > 2 && push!(ret.args, local_sym)
    end
    push!(epilogue.args, ret)
    push!(prologue.args, exprs[end])
    append!(prologue.args, epilogue.args)
    return esc(Expr(:let, Expr(:block), prologue))
end

end # module


text:
"""
    fresnels(x)

Compute the Fresnel S integral of type 0 (integral of `sin(πt²/2)`).
"""
fresnels(x) = (1 + im) / 4 * (erf((1 + im)√π / 2 * x) - im * erf((1 - im)√π / 2 * x))

"""
    fresnelc(x)

Compute the Fresnel C integral of type 0 (integral of `cos(πt²/2)`).
"""
fresnelc(x) = (1 - im) / 4 * (erf((1 + im)√π / 2 * x) + im * erf((1 - im)√π / 2 * x))


text:
"""
    install(pkg_name)

Lorem ipsum dolor sit amet.
"""
function install(args::Dict)
  pkg_name = args["<pkg_name>"]

  install(pkg_name)
end

function install(pkg_name)
  Pkg.add(pkg_name)
end


text:
"""
function to find ground state MPS for MPO H given input guess A in left
canonical form (normalisation not necessary). Full matrix calc.
"""

function ground_state_search_full(A,H,maxsweeps)

    na = length(A)
    At = deepcopy(A)
    FL = Array{Array{Complex{Float64},3},1}(na)
    FR = Array{Array{Complex{Float64},3},1}(na)

    FL[1] = ones(1,1,1)
    for jj = 1:na-1
        FL[jj+1] = update_le(A[jj],H[jj],A[jj],FL[jj])
    end
    FR[na] = ones(1,1,1)

    for sweep = 1:maxsweeps

        #sweep left
        for jj = na:-1:2
            HM = reshape(O_prod_LR(FL[jj],H[jj],FR[jj]),
                    size(FL[jj],1)*size(H[jj],2)*size(FR[jj],1),
                    size(FL[jj],1)*size(H[jj],2)*size(FR[jj],1))
            ae,av = eigs(HM; nev = 1, which = :SM, v0 = At[jj][:])
            println(ae)
            Q,R = qr(reshape(av,size(FL[jj],1),size(H[jj],2)*size(FR[jj],1))')
            At[jj] = reshape(Q',size(Q,2),size(At[jj],2),size(At[jj],3))
            FR[jj-1] = update_re(At[jj],H[jj],At[jj],FR[jj])
        end
            HM = reshape(O_prod_LR(FL[1],H[1],FR[1]),
                    size(FL[1],1)*size(H[1],2)*size(FR[1],1),
                    size(FL[1],1)*size(H[1],2)*size(FR[1],1))
            ae,av = eigs(HM; nev = 1, which = :SM, v0 = At[1][:])
            At[1] = reshape(av,size(FL[1],1),size(H[1],2),size(FR[1],1))

        #sweep right
        for jj = 1:na-1;
            HM = reshape(O_prod_LR(FL[jj],H[jj],FR[jj]),
                    size(FL[jj],1)*size(H[jj],2)*size(FR[jj],1),
                    size(FL[jj],1)*size(H[jj],2)*size(FR[jj],1))
            ae,av = eigs(HM; nev = 1, which = :SM, v0 = At[jj][:])
            println(ae)
            Q,R = qr(reshape(av,size(FL[jj],1)*size(H[jj],2),size(FR[jj],1)));
            At[jj] = reshape(Q,size(At[jj],1),size(At[jj],2),size(Q,2))
            FL[jj+1] = update_le(At[jj],H[jj],At[jj],FL[jj]);
        end
            HM = reshape(O_prod_LR(FL[na],H[na],FR[na]),
                    size(FL[na],1)*size(H[na],2)*size(FR[na],1),
                    size(FL[na],1)*size(H[na],2)*size(FR[na],1))
            ae,av = eigs(HM; nev = 1, which = :SM, v0 = At[na][:])
            At[na] = reshape(av,size(FL[na],1),size(H[na],2),size(FR[na],1))

    end

    no = sum(abs2(At[na][:]))
    return At, no

end

"""
function to find ground state MPS for MPO H given input guess A in left
canonical form (normalisation not necessary). Matrix as linear map.
"""

function ground_state_search!(At,H,maxsweeps,converge_rat,eigs_tol)

    na = length(At)
    FL = Array{Array{Complex{Float64},3},1}(na)
    FR = Array{Array{Complex{Float64},3},1}(na)
    ae = 0.0
    aestart = 0.0

    FL[1] = ones(1,1,1)
    for jj = 1:na-1
        FL[jj+1] = update_le(At[jj],H[jj],At[jj],FL[jj])
    end
    FR[na] = ones(1,1,1)

    for sweep = 1:maxsweeps

        println("sweeping left")
        #sweep left
        for jj = na:-1:2
            HM = LinearMap(x->A_prod_LOR_eigs(FL[jj],x,H[jj],FR[jj]),
            	size(At[jj][:],1),Complex128; ishermitian=true)
            ae,av = eigs(HM; nev = 1, which = :SR, v0 = At[jj][:],
                maxiter = 1000, tol = eigs_tol)
            println("site ", jj, " evalue is ", ae[1])
            if jj == na
                aestart = ae[1]
            end
            #Q,R = qr(reshape(av,size(FL[jj],1),size(H[jj],2)*size(FR[jj],1))')
            #av = qrfact!(reshape(av,size(FL[jj],1),size(H[jj],2)*size(FR[jj],1))')
            #av = full(av[:Q])
            av,rv = qr(reshape(av,size(FL[jj],1),size(H[jj],2)*size(FR[jj],1))')
            At[jj] = reshape(av',size(av,2),size(At[jj],2),size(At[jj],3))
            FR[jj-1] = update_re(At[jj],H[jj],At[jj],FR[jj])
            mo = At[jj-1]
            @tensor mo[-1,-2,-3] := mo[-1,-2,1]*rv[-3,1]
            At[jj-1] = mo
        end

        #HM = LinearMap(x->A_prod_LOR_eigs(FL[1],x,H[1],FR[1]),
        #    	size(At[1][:],1),Complex128; ishermitian=true)
        #ae,av = eigs(HM; nev = 1, which = :SR, v0 = At[1][:],
        #        maxiter = 1000, tol = eigs_tol)
        #println("site ", 1, " evalue is ", ae[1])
        #At[1] = reshape(av,size(FL[1],1),size(H[1],2),size(FR[1],1))

        println("sweeping right")
        #sweep right
        for jj = 1:na-1;
            HM = LinearMap(x->A_prod_LOR_eigs(FL[jj],x,H[jj],FR[jj]),
            	size(At[jj][:],1),Complex128; ishermitian=true)
            ae,av = eigs(HM; nev = 1, which = :SR, v0 = At[jj][:],
                maxiter = 1000, tol = eigs_tol)
            println("site ", jj, " evalue is ", ae[1])
            #Q,R = qr(reshape(av,size(FL[jj],1)*size(H[jj],2),size(FR[jj],1)));
            #av = qrfact!(reshape(av,size(FL[jj],1)*size(H[jj],2),size(FR[jj],1)))
            #av = full(av[:Q])
            av, rv = qr(reshape(av,size(FL[jj],1)*size(H[jj],2),size(FR[jj],1)))
            At[jj] = reshape(av,size(At[jj],1),size(At[jj],2),size(av,2))
            FL[jj+1] = update_le(At[jj],H[jj],At[jj],FL[jj]);
            mo = At[jj+1]
            @tensor mo[-1,-2,-3] := rv[-1,1]*mo[1,-2,-3]
            At[jj+1] = mo
        end

        HM = LinearMap(x->A_prod_LOR_eigs(FL[na],x,H[na],FR[na]),
            	size(At[na][:],1),Complex128; ishermitian=true)
        ae,av = eigs(HM; nev = 1, which = :SR, v0 = At[na][:],
                maxiter = 1000, tol = eigs_tol)
        println("site ", na, " evalue is ", ae[1])
        At[na] = reshape(av,size(FL[na],1),size(H[na],2),size(FR[na],1))

        if abs(ae[1]) > converge_rat*abs(aestart) || abs(ae[1])<1e-10
            return
        end

    end

    no = sum(abs2(At[na][:]))
    return ae

end


text:
using .MyExample
using Test

@testset "MyExample.jl" begin
    #2x+3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13

end


text:
__precompile__(false)

module RoutingNetworks

using LightGraphs, SFML, NearestNeighbors, Colors
import JLD, LibExpat

#main
export Node, Road, Network

#creation
export osm2network, subsetNetwork, removeNodes, singleNodes, inPolygon, roadTypeSubset
export stronglyConnected, intersections, queryOsmBox, queryOsmPolygon, updateProjection!
export loadTemplate, saveTemplate, isTemplate
export rectNetwork, squareNetwork, urbanNetwork

#routing
export RoutingPaths, roadDistances, getPathTimes, getPath, getPathEdges, shortestPaths!
export pathTime, pathEdgesTime, getPathWithTimes, longestPathTime
export parallelShortestPaths!, traveltime, pathRoads
export roadTypeRouting, randomTimeRouting, fixedSpeedTimes, maxSpeedTimes, uniformTimes, randomTimes
export allPairTimes

#visualization
export NetworkVisualizer, NetworkViz, RoutingViz
export visualize, visualInit, visualEvent, visualUpdate, visualRedraw, copyVisualData, redraw!
export VizColors, RoutingColors, RoadTypeColors, FadedColors, SpeedColors, RelativeSpeedColors

# projection onto network
export NetworkProjector, NetworkPosition, NodeProjector, NodePosition, RoadProjector, RoadPosition
export x, y, lat, lon, node, distToProj
#tools
export boundingBox, scaleBox, distanceGeo, distanceCoord, pointInsidePolygon, toENU, nRoads, nNodes

include("network.jl")

include("networkPosition.jl")

include("tools/geometry.jl")
include("tools/geodesy.jl")
include("tools/osmparser.jl")

include("creation/osm2network.jl")
include("creation/query.jl")
include("creation/subsets.jl")
include("creation/templates.jl")
include("creation/rect.jl")
include("creation/urban.jl")

include("routing/routingpaths.jl")
include("routing/shortestpaths.jl")
include("routing/tools.jl")
include("routing/virtualtimes.jl")

include("visualization/vizcolors.jl")
include("visualization/visualize.jl")
include("visualization/networkviz.jl")
include("visualization/routingviz.jl")

end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
include((@__DIR__)*"/../preamble.jl")

include("reflecting_model/model_def.jl")
include("default_params.jl")

pth = mkpath("hitting_times_model/reflecting_model/transient_distribution/point_mass/datat2/errors")

ic_string = "point_mass"

let
    idx = [1
    1001
    2001
    3001
    4001
    5001
    6001
    7001
    8001
    9001
   10001]
    ks_error = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    l1_error = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    cell_probs_error = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    ks_error_lwr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    l1_error_lwr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    cell_probs_error_lwr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    ks_error_upr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    l1_error_upr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    cell_probs_error_upr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])

    sim_cdf = CSV.read("hitting_times_model/reflecting_model/transient_distribution/point_mass/datat2/"*
                    "sims_cdf_evaluated.csv",DataFrame)
    for order in orders
        ks_error_row = zeros(4)
        l1_error_row = zeros(4)
        cell_probs_error_row = zeros(4)
        ks_error_row_lwr = zeros(4)
        l1_error_row_lwr = zeros(4)
        cell_probs_error_row_lwr = zeros(4)
        ks_error_row_upr = zeros(4)
        l1_error_row_upr = zeros(4)
        cell_probs_error_row_upr = zeros(4)
        for (c,model_string) in enumerate(string.([approx_types...]).*string.(1:4))
            approx_data = CSV.read("hitting_times_model/reflecting_model/transient_distribution/point_mass/datat2/approx_cdf/order_"*string(order)*"model_"*model_string*".csv",DataFrame)

            ks_error_row[c] = max(
                maximum(abs.(approx_data.phase_1-sim_cdf.phase_1)),
                maximum(abs.(approx_data.phase_2-sim_cdf.phase_2)),
            )
            l1_error_row[c] = sum(abs.(approx_data.phase_1[1:end-1]+approx_data.phase_1[2:end]-
                    sim_cdf.phase_1[1:end-1]-sim_cdf.phase_1[2:end]).*diff(approx_data.x))/2 +
                sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                    sim_cdf.phase_2[1:end-1]-sim_cdf.phase_2[2:end]).*diff(approx_data.x))/2

            cell_probs_error_row[c] = sum(abs.(diff(approx_data.phase_1[idx])-diff(sim_cdf.phase_1[idx]))) +
                sum(abs.(diff(approx_data.phase_2[idx])-diff(sim_cdf.phase_2[idx])))

            ks_error_samples = zeros(n_boot)
            l1_error_samples = zeros(n_boot)
            cell_probs_error_samples = zeros(n_boot)
            for n in 1:n_boot
                cdf_boot = CSV.read("hitting_times_model/reflecting_model/transient_distribution/point_mass/datat2/bootstrap/sample_"*string(n)*".csv",DataFrame)

                ks_error_samples[n] = max(
                    maximum(abs.(approx_data.phase_1-cdf_boot.phase_1)),
                    maximum(abs.(approx_data.phase_2-cdf_boot.phase_2)),
                )
                l1_error_samples[n] = sum(abs.(approx_data.phase_1[1:end-1]+approx_data.phase_1[2:end]-
                        cdf_boot.phase_1[1:end-1]-cdf_boot.phase_1[2:end]).*diff(approx_data.x))/2 +
                    sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                        cdf_boot.phase_2[1:end-1]-cdf_boot.phase_2[2:end]).*diff(approx_data.x))/2

                cell_probs_error_samples[n] = sum(abs.(diff(approx_data.phase_1[idx])-diff(cdf_boot.phase_1[idx]))) +
                    sum(abs.(diff(approx_data.phase_2[idx])-diff(cdf_boot.phase_2[idx])))
            end
            ks_q5 = quantile(ks_error_samples,qtiles[1])
            ks_q95 = quantile(ks_error_samples,qtiles[2])
            l1_q5 = quantile(l1_error_samples,qtiles[1])
            l1_q95 = quantile(l1_error_samples,qtiles[2])
            cp_q5 = quantile(cell_probs_error_samples,qtiles[1])
            cp_q95 = quantile(cell_probs_error_samples,qtiles[2])
            ks_error_row_lwr[c] = ks_q5
            ks_error_row_upr[c] = ks_q95
            l1_error_row_lwr[c] = l1_q5
            l1_error_row_upr[c] = l1_q95
            cell_probs_error_row_lwr[c] = cp_q5
            cell_probs_error_row_upr[c] = cp_q95
        end
        push!(ks_error_lwr,ks_error_row_lwr)
        push!(ks_error_upr,ks_error_row_upr)
        push!(l1_error_lwr,l1_error_row_lwr)
        push!(l1_error_upr,l1_error_row_upr)
        push!(cell_probs_error_lwr,cell_probs_error_row_lwr)
        push!(cell_probs_error_upr,cell_probs_error_row_upr)
        push!(ks_error,ks_error_row)
        push!(l1_error,l1_error_row)
        push!(cell_probs_error,cell_probs_error_row)
    end
    CSV.write(pth*"/ks_lwr_"*ic_string*".csv",ks_error_lwr)
    CSV.write(pth*"/ks_upr_"*ic_string*".csv",ks_error_upr)
    CSV.write(pth*"/l1_lwr_"*ic_string*".csv",l1_error_lwr)
    CSV.write(pth*"/l1_upr_"*ic_string*".csv",l1_error_upr)
    CSV.write(pth*"/cell_probs_lwr_"*ic_string*".csv",cell_probs_error_lwr)
    CSV.write(pth*"/cell_probs_upr_"*ic_string*".csv",cell_probs_error_upr)
    CSV.write(pth*"/ks_"*ic_string*".csv",ks_error_lwr)
    CSV.write(pth*"/l1_"*ic_string*".csv",l1_error_lwr)
    CSV.write(pth*"/cell_probs_"*ic_string*".csv",cell_probs_error)
end


text:
using Base.Test
using AGRN

@testset "Chromosome test" begin
    config = AGRN.Config()
    nin = rand(1:10)
    nout = rand(1:10)
    nreg = rand(1:10)
    grn = GRN(nin, nout, nreg, config)
    set_input!(grn, zeros(nin))
    step!(grn, 10)
    @test all(grn.cons[1:nin] .== 0.0)
    new_grn = GRN(grn)
    @test all(new_grn.cons .!= grn.cons)
    reset!(grn)
    @test all(new_grn.cons .== grn.cons)
    for i in 1:100
        inputs = rand(nin)
        step!(grn, inputs)
        step!(new_grn, inputs)
        @test all(new_grn.cons .== grn.cons)
    end
    @test distance(grn, new_grn, config) == 0.0
end

@testset "Import/export test" begin
    config = AGRN.Config()
    nin = rand(1:10)
    nout = rand(1:10)
    nreg = rand(1:10)
    grn = GRN(nin, nout, nreg, config)
    step!(grn, 10)
    grn_str = JSON.json(grn)
    clone = GRN(grn_str)
    @test all(grn.ids .== clone.ids)
    @test all(grn.enh .== clone.enh)
    @test all(grn.inh .== clone.inh)
    @test grn.beta == clone.beta
    @test grn.delta == clone.delta
    @test any(grn.cons .!= clone.cons)
end