text:
###
# This file contains a detailed example, demonstrating the typical workflow
# Graft aims to support.
# The dataset used here was constructed by splicing together two separate datasets:
#
# 1. `SOCR Data MLB HeightsWeights`: Heights, ages and weights of Baseball players (Vertex Data). References:
#   * Jarron M. Saint Onge, Patrick M. Krueger, Richard G. Rogers. (2008) Historical trends in height,
#   weight, and body mass: Data from U.S. Major League Baseball players, 1869-1983, Economics & Human
#   Biology, Volume 6, Issue 3, Symposium on the Economics of Obesity, December 2008, Pages 482-488,
#   ISSN 1570-677X, DOI: 10.1016/j.ehb.2008.06.008.
#   * Jarron M. Saint Onge, Richard G. Rogers, Patrick M. Krueger. (2008) Major League Baseball Players'
#   Life Expectancies, Southwestern Social Science Association, Volume 89, Issue 3, pages 817–830,
#   DOI: 10.1111/j.1540-6237.2008.00562.x.
# 2. Advogato Trust Network : Edge weights between 0 and 1. References:
#   * Advogato network dataset -- KONECT, July 2016. [http](http://konect.uni-koblenz.de/networks/advogato)
#   * Paolo Massa, Martino Salvetti, and Danilo Tomasoni. Bowling alone and trust decline in social network
#   sites. In Proc. Int. Conf. Dependable, Autonomic and Secure Computing, pages 658--663, 2009.
#
# The dataset has 6541 vertices, 51127 edges.
# Vertex properties: Age, Height(cm), Weight(kg)
# Edge properties  : Trust(float)
###

using Graft
using StatsBase
import LightGraphs

## Load and summarize the graph.

# Load the graph
download(
 "https://raw.githubusercontent.com/pranavtbhat/Graft.jl/gh-pages/Datasets/baseball.txt",
 joinpath(Pkg.dir("Graft"), "examples/baseball.txt")
)
g = loadgraph(joinpath(Pkg.dir("Graft"), "examples/baseball.txt"))

# Get the graph's size
size(g)

# Get an edge iterator for g
edges(g)

# List vertex labels
encode(g)

# Split the graph into vertex and edge descriptors
V,E = g;

# Display the vertex table
V

# Display the edge table
E

## Run some metadata queries

# Find the average BMI of baseball players
@query(g |> eachvertex(v.Weight / (v.Height / 100) ^ 2)) |> mean

# Find the median height of baseball players in their 20s
@query(g |> filter(v.Age < 30,v.Age >= 20) |> eachvertex(v.Height * 0.0328084)) |> median

# Find the mean age difference in strong relationships
@query(g |> filter(e.Trust > 0.8) |> eachedge(s.Age - t.Age)) |> abs |> mean

# Find fred's 3 hop neighborhood (friends and friends-of-friends and so on)
fred_nhood = hopgraph(g, "fred", 3)

# See how well younger players in fred's neighborhood trust each other
@query(fred_nhood |> filter(v.Age > 30) |> eachedge(e.Trust)) |> mean

## Export data to LightGraphs.jl

# Find the 2 hop neighborhood of 2 separate vertices (multi seed traversal)
sg = hopgraph(g, ["nikolay", "jbert"], 3)

# Generate an edge distance property on the inverse of normalized-trust
dists = @query(sg |> eachedge(1 / e.Trust ));
seteprop!(sg, :, dists, :Dist);

# Trim edges of very high distance
sg = @query(sg |> filter(e.Dist < 10))

# Export the graph's adjacency matrix
M = export_adjacency(sg)
lg = LightGraphs.DiGraph(M)

# Export the edge distance property
D = export_edge_property(sg, :Dist)

# Compute betweenness centrailty
centrality = LightGraphs.betweenness_centrality(lg)

# Set the centrality as a vertex property
setvprop!(sg, :, centrality, :Centrality)

# Apply all pair shortest paths on the graph
apsp = LightGraphs.floyd_warshall_shortest_paths(lg, D).dists;

# Add the new shortest paths as a property to the graph
eit = edges(sg);
seteprop!(sg, :, [apsp[e.second,e.first] for e in eit], :Shortest_Dists);

# Show new vertex descriptor
VertexDescriptor(sg)

# Show the new edge descriptor
EdgeDescriptor(sg)


text:
using MasterEquationSolve
using MasterEquationSolve: qt, qobj_to_jl, mesolve, mesolve_qutip, myqload
using Printf
using Test

# make the destroy operators in a tensor product space of dimensions dims
function destroy_ops(dims)
    ops = qt.qeye.(dims)
    return ntuple(length(dims)) do i
        ops[i] = qt.destroy(dims[i])
        A = qt.tensor(ops)
        ops[i] = qt.qeye(dims[i])
        A
    end
end

@testset "MasterEquationSolve.jl" begin
    ωe = 1.5
    ωc = 1.45
    g = 0.5
    γc = 0.4
    γdephasing = 0.5
    ts = LinRange(0,20,400)

    Nmol = 4
    Nc = Nmol+1

    (acav, σs...) = destroy_ops([Nc;fill(2,Nmol)])

    H0 = ωc * acav.dag() * acav + sum(ωe*σ.dag()*σ + g*(σ.dag()*acav + σ*acav.dag()) for σ in σs)

    # This is the dephasing operator (2 σp σm - 1)
    makedeph(σ) = sqrt(γdephasing)*(2*σ.dag()*σ-1)
    c_ops = [sqrt(γc)*acav, makedeph.(σs)...];
    e_ops = [σs[1].dag()*σs[1], acav.dag()*acav]
    ψ0 = qt.tensor([qt.basis(Nc); qt.basis.(fill(2,Nmol),1)])

    infile = "mesolve_runtests_in"
    outdir = "mesolve_runtests_out"
    qt.qsave((H0,c_ops,ψ0,ts,e_ops),infile)

    sol_qt = qt.mesolve(H0,ψ0,ts,c_ops,options=qt.Options(atol=1e-12,rtol=1e-10))
    sol_qt_e = qt.mesolve(H0,ψ0,ts,c_ops,e_ops=e_ops,options=qt.Options(atol=1e-12,rtol=1e-10))

    H = qobj_to_jl(H0)
    J = qobj_to_jl.(c_ops)
    ρ0 = qobj_to_jl(ψ0,true);
    e_ops_j = qobj_to_jl.(e_ops)

    for backend = (:CPU,:CUDA)
        sol, sv = mesolve(H,J,ρ0,ts;backend=backend,saveddata=:full)
        @test all([s ≈ q.full() for (s,q) in zip(sv.saveval,sol_qt.states)])

        sol, sv = mesolve(H,J,ρ0,ts;backend=backend,saveddata=e_ops_j)
        @test reduce(hcat,sol_qt_e.expect) ≈ transpose(reduce(hcat,sv.saveval))

        rm(outdir,force=true,recursive=true)
        sol, sv = mesolve_qutip(infile,outdir; backend=backend,saveddata=:full)
        @test all([s ≈ q.full() for (s,q) in zip(sv.saveval,sol_qt.states)])
        sv2 = [myqload(@sprintf("%s/%08d",outdir,i)) for i in 1:length(ts)]
        @test ts == [s[1] for s in sv2]
        @test all([s[2] ≈ q.full() for (s,q) in zip(sv2,sol_qt.states)])

        rm(outdir,force=true,recursive=true)
        # not passing saveddata corresponds to taking it from the saved qutip file
        sol, sv = mesolve_qutip(infile,outdir; backend=backend)
        @test reduce(hcat,sol_qt_e.expect) ≈ transpose(reduce(hcat,sv.saveval))
        sv2 = [myqload(@sprintf("%s/%08d",outdir,i)) for i in 1:length(ts)]
        @test ts == [s[1] for s in sv2]
        @test reduce(hcat,sol_qt_e.expect) ≈ transpose(reduce(hcat,[s[2] for s in sv2]))
    end

    rm(outdir,force=true,recursive=true)
    rm(infile*".qu",force=true)
end

ext: jl


text:
module PlutoDevMacros

using MacroTools
using Requires

# This hack is the only way I found to access @skip_as_script as PlutoDevMacros.@skip_as_script
# @eval const $(Symbol("@skip_as_script")) = $(Symbol("@plutohookskip"))

# export @only_in_nb, @only_out_nb, include_mapexpr, @skip_as_script
# export notebook_to_source

include("../notebooks/basics.jl") # @only_in_nb, @only_out_nb, is_notebook_local, plutodump, @current_pluto_cell_id, @current_pluto_notebook_file
include("../notebooks/mapexpr.jl") # hasexpr, default_exprlist, include_mapexpr
include("../notebooks/plutoinclude_macro.jl") # hasexpr, default_exprlist, include_mapexpr
include("../notebooks/pluto_traits.jl") # This defines and exports the @plutotraits macro

# function __init__()
# 	@require WhereTraits="c9d4e05b-6318-49cb-9b56-e0e2b0ceadd8" include("../notebooks/pluto_traits.jl") # This defines and exports the @plutotraits macro
# end

end # module


text:
const pi_token = Token{:reserved}("pi")

divpi(n::Int) = (pi_token, Token{:reserved}("/"), qasm_int(n))

# Intrinsics
# // 3-parameter 2-pulse single qubit gate
# gate u3(theta,phi,lambda) q { U(theta,phi,lambda) q; }
# // 2-parameter 1-pulse single qubit gate
# gate u2(phi,lambda) q { U(pi/2,phi,lambda) q; }
# // 1-parameter 0-pulse single qubit gate
# gate u1(lambda) q { U(0,0,lambda) q; }
# // Pauli gate: bit-flip
# gate x a { u3(pi,0,pi) a; }
# // Pauli gate: bit and phase flip
# gate y a { u3(pi,pi/2,pi/2) a; }
# // Pauli gate: phase flip
# gate z a { u1(pi) a; }
# // Clifford gate: Hadamard
# gate h a { u2(0,pi) a; }
# // Clifford gate: sqrt(Z) phase gate
# gate s a { u1(pi/2) a; }
# // Clifford gate: conjugate of sqrt(Z)
# gate sdg a { u1(-pi/2) a; }
# // C3 gate: sqrt(S) phase gate
# gate t a { u1(pi/4) a; }
# // C3 gate: conjugate of sqrt(S)
# gate tdg a { u1(-pi/4) a; }

qasm_gate_name(::Type{Intrinsics.XGate}) = "x"
qasm_gate_name(::Type{Intrinsics.YGate}) = "y"
qasm_gate_name(::Type{Intrinsics.ZGate}) = "z"
qasm_gate_name(::Type{Intrinsics.HGate}) = "h"
qasm_gate_name(::Type{Intrinsics.SGate}) = "s"
qasm_gate_name(::Type{Intrinsics.TGate}) = "t"
qasm_gate_name(::Type{<:Intrinsics.Rx}) = "rx"
qasm_gate_name(::Type{<:Intrinsics.Ry}) = "ry"
qasm_gate_name(::Type{<:Intrinsics.Rz}) = "rz"

function intrinsic_qasm(::Type{Intrinsics.XGate}, qargs)
    UGate(pi_token, qasm_int(0), pi_token, qargs[1])
end

function intrinsic_qasm(::Type{Intrinsics.YGate}, qargs)
    UGate(pi_token, divpi(2), divpi(2), qargs[1])
end

function intrinsic_qasm(::Type{Intrinsics.ZGate}, qargs)
    UGate(qasm_int(0), qasm_int(0), pi_token, qargs[1])
end

function intrinsic_qasm(::Type{Intrinsics.HGate}, qargs)
    UGate(divpi(2), qasm_int(0), pi_token, qargs[1])
end

function intrinsic_qasm(::Type{Intrinsics.SGate}, qargs)
    UGate(qasm_int(0), qasm_int(0), divpi(2), qargs[1])
end

function intrinsic_qasm(::Type{Intrinsics.TGate}, qargs)
    UGate(qasm_int(0),qasm_int(0),divpi(4), qargs[1])
end

# // Rotation around X-axis
# gate rx(theta) a { u3(theta,-pi/2,pi/2) a; }
# // rotation around Y-axis
# gate ry(theta) a { u3(theta,0,0) a; }
# // rotation around Z axis
# gate rz(phi) a { u1(phi) a; }
# NOTE: theta can be a variable or exp
function intrinsic_qasm(::Type{<:Intrinsics.Rx}, theta, qargs)
    UGate(theta, Neg(divpi(2)), divpi(2), qargs[1])
end

function intrinsic_qasm(::Type{<:Intrinsics.Ry}, theta, qargs)
    UGate(qasm_int(0), theta, qasm_int(0), qargs[1])
end

function intrinsic_qasm(::Type{<:Intrinsics.Rz}, theta, qargs)
    UGate(qasm_int(0), qasm_int(0), theta, qargs[1])
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SimpleAFMachine <: Machine
        R::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Xd_pp::Float64
        Xq_pp::Float64
        Td0_p::Float64
        Tq0_p::Float64
        Td0_pp::Float64
        Tq0_pp::Float64
        MVABase::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of 4-states simplified Anderson-Fouad (SimpleAFMachine) model.
 The derivative of stator fluxes (ψd and ψq) is neglected and ωψd = ψd and
 ωψq = ψq is assumed (i.e. ω=1.0). This is standard when transmission network
 dynamics is neglected.
 If transmission dynamics is considered use the full order Anderson Fouad model.

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: (0, nothing)
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_pp::Float64`: Sub-Transient reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq_pp::Float64`: Sub-Transient reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: (0, nothing)
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: (0, nothing)
- `Td0_pp::Float64`: Time constant of sub-transient d-axis voltage, validation range: (0, nothing)
- `Tq0_pp::Float64`: Time constant of sub-transient q-axis voltage, validation range: (0, nothing)
- `MVABase::Float64`: Nominal Capacity in MVA, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SimpleAFMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Sub-Transient reactance after EMF in d-axis per unit"
    Xd_pp::Float64
    "Sub-Transient reactance after EMF in q-axis per unit"
    Xq_pp::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    "Time constant of sub-transient d-axis voltage"
    Td0_pp::Float64
    "Time constant of sub-transient q-axis voltage"
    Tq0_pp::Float64
    "Nominal Capacity in MVA"
    MVABase::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SimpleAFMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, MVABase, ext=Dict{String, Any}(), )
    SimpleAFMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, MVABase, ext, [:eq_p, :ed_p, :eq_pp, :ed_pp], 4, InfrastructureSystemsInternal(), )
end

function SimpleAFMachine(; R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, MVABase, ext=Dict{String, Any}(), )
    SimpleAFMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, MVABase, ext, )
end

# Constructor for demo purposes; non-functional.
function SimpleAFMachine(::Nothing)
    SimpleAFMachine(;
        R=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Xd_pp=0,
        Xq_pp=0,
        Td0_p=0,
        Tq0_p=0,
        Td0_pp=0,
        Tq0_pp=0,
        MVABase=0,
        ext=Dict{String, Any}(),
    )
end

"""Get SimpleAFMachine R."""
get_R(value::SimpleAFMachine) = value.R
"""Get SimpleAFMachine Xd."""
get_Xd(value::SimpleAFMachine) = value.Xd
"""Get SimpleAFMachine Xq."""
get_Xq(value::SimpleAFMachine) = value.Xq
"""Get SimpleAFMachine Xd_p."""
get_Xd_p(value::SimpleAFMachine) = value.Xd_p
"""Get SimpleAFMachine Xq_p."""
get_Xq_p(value::SimpleAFMachine) = value.Xq_p
"""Get SimpleAFMachine Xd_pp."""
get_Xd_pp(value::SimpleAFMachine) = value.Xd_pp
"""Get SimpleAFMachine Xq_pp."""
get_Xq_pp(value::SimpleAFMachine) = value.Xq_pp
"""Get SimpleAFMachine Td0_p."""
get_Td0_p(value::SimpleAFMachine) = value.Td0_p
"""Get SimpleAFMachine Tq0_p."""
get_Tq0_p(value::SimpleAFMachine) = value.Tq0_p
"""Get SimpleAFMachine Td0_pp."""
get_Td0_pp(value::SimpleAFMachine) = value.Td0_pp
"""Get SimpleAFMachine Tq0_pp."""
get_Tq0_pp(value::SimpleAFMachine) = value.Tq0_pp
"""Get SimpleAFMachine MVABase."""
get_MVABase(value::SimpleAFMachine) = value.MVABase
"""Get SimpleAFMachine ext."""
get_ext(value::SimpleAFMachine) = value.ext
"""Get SimpleAFMachine states."""
get_states(value::SimpleAFMachine) = value.states
"""Get SimpleAFMachine n_states."""
get_n_states(value::SimpleAFMachine) = value.n_states
"""Get SimpleAFMachine internal."""
get_internal(value::SimpleAFMachine) = value.internal


text:
#This is the main file. See Readme.md for how to run the program.

#It is assumed that the working directory is the directory where main.jl is located and contains file Project.toml.

#Activate project environment and install all dependent packages.
#(same versions as the ones used in development)

#Comment this out during development (debugger is a lot faster this way)
import Pkg

#Can comment this out during debugging to save time
Pkg.activate(".")
Pkg.instantiate()

#Load necessary packages
using Optim, Parameters, QuantEcon, BenchmarkTools,JLD, Interpolations,Plots
using .Threads #so we don't have to write Threads.@threads every time

#Import the local module brexit_investment which contains all new functions, data types, etc. used in the project
import brexit_investment
using brexit_investment

#Run start-up script (see the file for details of what it does)
#Most importantly it loads parameters from file and collects these in array par.
#Some outputs are N_S (number of parametrisations/equilibria), N_th (num of threads), par (array of parameters).
include("./src/startup.jl")

#=
The main body of the program follows. The algorithm proceeds in the following stages:

(0) Prep work. Generate objects used throughout such as objects for storing output.

(1) Compute stationary equilibria for the initial state (pre-referendum) and each of the long-run post-Brexit states.

(2) Compute transition paths.

(3) Compute various statistics, plot results, save results
=#

#*************(0) Prep work**********************

#=
#Set up storage for stationary equilibria and transition paths

#SE - array of stationary equilibria (one for each element of par_diff).

#TP -transition paths
 This is a matrix of stationary equilibria struct. Calling it a stationary equilibrium is not accurate but it is essentially a path of distribution, value functions, etc. over time, where in each period we store the same objects that would be stored in a stationary equilibrium.

#It is assumed that:
1) the number of grid points etc. are the same for all parametrisations, values from par[1] are used.
2) the length of all transition paths is the same and equal to the value in par[1].T_max)
=#

if ((@isdefined loadData) && (@isdefined loadFolder) && loadData)
    #loading SE, TP from files
    println("Loading initial guess from folder results/$loadFolder.")
    SE = load("results/$loadFolder/SE.jld","SE")
    TP = load("results/$loadFolder/TP.jld","TP")
    dataloaded = true
else
    #Initialisation with default values

    #Initial guess for value function 0: this leads to a corner solution in the first iteration (choose capital equal to the first point on the grid). This is usually fine but may not be, depending on the optimisation subroutine that we choose.
    #Ultimately, the initial guess does not matter (checked numerically)
    SE = fill(stat_equil(N_kh = par[1].N_kh,N_k = par[1].N_k,N_z = par[1].N_z),N_S);

    #Alternative - initial guess for value function is log of capital. This avoids the corner solution in the first iteration.
    #SE = fill(stat_equil(N_kh = par[1].N_kh,N_k = par[1].N_k,N_z = par[1].N_z,V=log.(fillV(par[1].N_k,par[1].N_z,par[1].k_gr))),N_S);

    TP = fill(stat_equil(N_kh = par[1].N_kh,N_k = par[1].N_k,N_z = par[1].N_z),par[1].T_max,N_S);
    dataloaded = false
end

#JIT compilation acceleration (first call the functions with tiny arguments so the compiler precompiles functions). Subsequent callSEs are fast.
partiny = pars(N_z = 2,N_k = 2,VFI_maxiter=1,SE_maxiter=1)
SEtiny = stat_equil(N_z=2,N_k=2,N_kh=2)
SE_compute!(partiny,SEtiny)

#************(1) Stationary equilibrium ***************
#Compute stationary equilibria for each set of parameters.
println("Computing Stationary Equilibria...\n")
for i=1:N_S
println("________________________________
Equilibrium $i out of $N_S")
    #SE[i] is initial guess in the first iteration. In the following iteration, the initial guess is the stationary equilibrium computed in the first call (SE[1]), unless we loaded the guess from file (which should be better)
    if(!dataloaded)
        SE[i] = SE_compute!(par[i],SE[i])
    else
        SE[i] = SE_compute!(par[i],SE[1])
    end

    #In development, uncomment this to only compute one stationary equilibrium to save time
    #print("Only computing the first stationary equilibrium./n")
    #break
end

#************(2) Transition paths************************

#As initial guess for value function, we use the value functions from stationary equilibria.

#************(3) Statistics, Plots*****************

#Generate plots

#(The historgram plot)
#range for plots (grid points indices)
a1 = 1
a2 = par[1].N_k
zind = Int(round(par[1].N_z/2)) #index of shock for which we plot the distribution
zval = par[1].shock_mc.state_values[zind]
p = plot(par[1].k_gr[a1:a2],SE[1].μ[a1:a2,1],title = "μ", xlabel="k", label = "z = $zval",lw = 2)
display(p)
savefig("$foldername/mu.pdf")

#Plot of the policy function and value function
a1 = 1
a2 = par[1].N_k
zind = Int(round(par[1].N_z/2)) #index of shock for which we plot the function
zval = par[1].shock_mc.state_values[zind]
pV = plot(par[1].k_gr[a1:a2],SE[1].V[a1:a2,1],title = "V", xlabel="k", label = "z = $zval",lw = 2)
display(pV)
savefig("$foldername/V.pdf")

#Also plot ξc (investment cost threshold):
a1 = 1
a2 = par[1].N_k
zind = Int(round(par[1].N_z/2)) #index of shock for which we plot the function
zval = par[1].shock_mc.state_values[zind]
pξ = plot(par[1].k_gr[a1:a2],SE[1].ξc[a1:a2,1],title = "ξ^c (investment cost threshold)", xlabel="k", label = "z = $zval",lw = 2)
display(pξ)
savefig("$foldername/xi.pdf")

#Saving results
#(function saveAll is defined in brexTools.jl, see there for details)
saveAll(foldername,SE,TP)

#=

=#


text:
@testset "nc_dims" begin
    f = "data/temp_HI.nc"
    # f = "temp_HI.nc"
    dims = nc_dims(f)
    @test names(dims) == ["lon", "lat", "time"]
end
# ds = nc_open(f)

@testset "ncatt_put" begin
    f = "data/temp_HI.nc"
    # f = "temp_HI.nc"
    # ds = nc_open(f)
    # close(ds)

    atts = Dict("a" => 1, "b" => 2)
    ncatt_put(f, atts)
    # nc_atts(f) == atts
    nc_info(f)

    ks = keys(atts) |> collect
    ncatt_del(f, ks)
    nc_info(f)
    # @test names(dims) == ["lon", "lat", "time"]
end


text:
@testset verbose = true "Point processes" begin
    @testset "Interface" begin
        pp = NaiveMultivariatePoissonProcess([1., 2., 3.])
        h = rand(pp, 0., 100.)
        pp_init = NaiveMultivariatePoissonProcess(ones(3))
        pp_est1 = fit(pp_init, h)
        error1 = intensity(pp_est1) - intensity(pp)
        @test maximum(error1) < 0.1
        pp_est2 = fit(pp_init, h, adtype=GalacticOptim.AutoForwardDiff(), alg=Optim.LBFGS())
        error2 = intensity(pp_est2) - intensity(pp)
        @test maximum(error2) < 0.1
    end
end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials [MDP]
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen [cuongng@mit.edu, exapde@gmail.com]
#****************************************************************************

function boundingbox2D(pbc::Array{Int64,1}, a::Array{Float64,1}, b::Array{Float64,1}, r::Array{Float64,1})
    # a = [a1, a2] and b = [b1,b2] are two vectors defining a parallelogram
    # r is the distance from the original parallelogram to the bounding parallelogram
    # w are four vertices defining the bounding parallelogram

    if length(r)==1
        r = [r; r];
    end

    a = a[:]; b = b[:];

    a1 = a[1]; a2 = a[2];
    b1 = b[1]; b2 = b[2];
    norma = norm(a);
    normb = norm(b);

    # vertices of the parallelogram defined by a and b
    v1 = [0.0; 0.0];
    v2 = a;
    v3 = a+b;
    v4 = b;

    # a2*x - a1*y = 0 -> y = a2/a1 * x [1st principal axis]
    # a2*x - a1*y = da a parallel line to 1st principal axis
    # The distance between the two lines is da/norma. As da/norma = r, we get
    da = r[2]*norma;
    # a2*x - a1*y = da is a line which is below and parallel to the 1st axis at a distance r

    # apply the same formula to the second axis b
    db = r[1]*normb;
    # b2*x - b1*y = -db is a line which is left and parallel to the 2nd axis at a distance r

    # intersection of a2*x - a1*y = da and b2*x - b1*y = -db
    w1 = [a2 -a1; b2 -b1]\[da; -db]; # the 1st vertex of the bounding parallelogram

    # find e = [e1,e2] such that e1*a + e2*b = w1
    e = [a b]\w1;

    # distance between w1 and e[1]*a
    sb = norm(w1-e[1]*a);
    # distance between w1 and e[2]*b
    sa = norm(w1-e[2]*b);

    # # length of the bounding parallelogram along the 1st axis
    # l1 = norm[a] + 2*sa;
    # # length of the bounding parallelogram along the 2nd axis
    # l2 = norm[b] + 2*sb;

    # length of the bounding parallelogram along the 1st axis
    l1 = norm(a) + 2*sa*pbc[1];
    # length of the bounding parallelogram along the 2nd axis
    l2 = norm(b) + 2*sb*pbc[2];

    # the 1st vertex of  the bounding parallelepiped
    w1 = pbc[1]*e[1]*a + pbc[2]*e[2]*b;
    # the 2nd vertex of  the bounding parallelogram
    w2 = w1 + l1*a/norm(a);
    # the 3rd vertex of  the bounding parallelogram
    w3 = v3 - w1;
    # the 4th vertex of  the bounding parallelogram
    w4 = w1 + l2*b/norm(b);

    v = [v1 v2 v3 v4];
    w = [w1 w2 w3 w4];

    return v, w

end

function boundingbox3D(pbc::Array{Int64,1}, a::Array{Float64,1}, b::Array{Float64,1}, c::Array{Float64,1}, r::Array{Float64,1})
    # a = (a1, a2, a3), b = (b1,b2,b3), c=(c1,c2,c3) are three vectors defining a parallelepiped
    # r is the distance from the original parallelepiped to the bounding parallelepiped
    # w are 8 vertices defining the bounding parallelepiped

    if length(r)==1
        r = [r; r; r];
    end

    a = a[:]; b = b[:]; c = c[:];

    a1 = a[1]; a2 = a[2]; a3 = a[3];
    b1 = b[1]; b2 = b[2]; b3 = b[3];
    c1 = c[1]; c2 = c[2]; c3 = c[3];
    norma = norm(a);
    normb = norm(b);
    normc = norm(c);

    # vertices of the parallelogram defined by a, b, and c
    v1 = [0; 0; 0];
    v2 = a;
    v3 = a+b;
    v4 = b;
    v5 = c;
    v6 = a+c;
    v7 = a+b+c;
    v8 = b+c;

    # the 1st plane defined by a and b
    p1 = [a2*b3-a3*b2; a3*b1-a1*b3; a1*b2-a2*b1];
    normp1 = sqrt(p1[1]*p1[1] + p1[2]*p1[2] + p1[3]*p1[3]);

    # Since the distance between (0,0,0) and the 1st plane p1[1]*x + p1[2]*y + p1[3]*z = d1
    # is equal to r[3], we have
    d1 = -sign(p1[3])*r[3]*normp1;

    # the 2nd plane defined by b and c
    p2 = [b2*c3-b3*c2; b3*c1-b1*c3; b1*c2-b2*c1];
    normp2 = sqrt(p2[1]*p2[1] + p2[2]*p2[2] + p2[3]*p2[3]);

    # Since the distance between (0,0,0) and the 2nd plane p2[1]*x + p2[2]*y + p2[3]*z = d2
    # is equal to r[1], we have
    d2 = -sign(p2[1])*r[1]*normp2;

    # the 3rd plane defined by c and a
    p3 = [c2*a3-c3*a2; c3*a1-c1*a3; c1*a2-c2*a1];
    normp3 = sqrt(p3[1]*p3[1] + p3[2]*p3[2] + p3[3]*p3[3]);

    # Since the distance between (0,0,0) and the 2nd plane p3(1)*x + p3(2)*y + p3(3)*z = d3
    # is equal to r[2], we have
    d3 = -sign(p3[2])*r[2]*normp3;

    # intersection of the above three planes
    A = reshape([p1; p2; p3],3,3)';
    w1 = A\[d1; d2; d3];

    # find e = (e1,e2,e3) such that e1*a + e2*b + e3*c = w1
    e = [a b c]\w1;

    # distance between w1 and the point e2*b + e3*c
    sa = norm(w1 - e[2]*b - e[3]*c);

    # distance between w1 and the point e1*a + e3*c
    sb = norm(w1 - e[1]*a - e[3]*c);

    # distance between w1 and the point e1*a + e2*b
    sc = norm(w1 - e[1]*a - e[2]*b);

    # length of the bounding parallelepiped along the 1st axis
    la = norma + 2*sa*pbc[1];
    # length of the bounding parallelepiped along the 2nd axis
    lb = normb + 2*sb*pbc[2];
    # length of the bounding parallelepiped along the 3rd axis
    lc = normc + 2*sc*pbc[3];

    # the 1st vertex of  the bounding parallelepiped
    w1 = pbc[1]*e[1]*a + pbc[2]*e[2]*b + pbc[3]*e[3]*c;
    # the 2nd vertex of  the bounding parallelepiped
    w2 = w1 + la*a/norma;
    # the 4th vertex of  the bounding parallelepiped
    w4 = w1 + lb*b/normb;
    # the 3rd vertex of  the bounding parallelepiped
    w3 = w2 + w4 - w1;
    # the 5th vertex of  the bounding parallelepiped
    w5 = w1 + lc*c/normc;
    # the 6th vertex of  the bounding parallelepiped
    w6 = w5 + la*a/norma;
    # the 8th vertex of  the bounding parallelepiped
    w8 = w5 + lb*b/normb;
    # the 7th vertex of  the bounding parallelepiped
    w7 = w6 + w8 - w5;

    v = [v1 v2 v3 v4 v5 v6 v7 v8];
    w = [w1 w2 w3 w4 w5 w6 w7 w8];

    return v, w

end

function boundingbox(pbc::Array{Int64,1}, rmax::Array{Float64,1}, a::Array{Float64,1}, b::Array{Float64,1}, c=nothing)
    # pbc periodic boundary conditions
    # a = (a1, a2, a3), b = (b1,b2,b3), c=(c1,c2,c3) are three vectors defining a parallelepiped
    # rmax is the distance from the original parallelepiped to the bounding parallelepiped
    # v are 8 vertices defining the original parallelepiped
    # w are 8 vertices defining the bounding parallelepiped

    if c===nothing
        v, w = boundingbox2D(pbc, a, b, rmax);
    else
        v, w = boundingbox3D(pbc, a, b, c, rmax);
    end

    return v, w
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module Arrow

immutable Bool
    value::UInt64
end

# Time types
immutable Second end
immutable Millisecond end
immutable Microsecond end
immutable Nanosecond end

# Timestamp type with time unit `P` and timezone `Z`
immutable Timestamp{P,Z}
    value::Int64
end

const UNIXEPOCH_TS = Dates.value(DateTime(1970)) #Rata Die milliseconds for 1970-01-01T00:00:00

scale(::Type{Second}, x) = 1000 * x.value
scale(::Type{Millisecond}, x) = x.value
scale(::Type{Microsecond}, x) = div(x.value,1000)
scale(::Type{Nanosecond}, x) =  div(x.value,1000000)

function unix2datetime{P}(::Type{P}, x)
    rata = UNIXEPOCH_TS + scale(P, x)
    return DateTime(Dates.UTM(rata))
end
function datetime2unix(x::DateTime)
    return Dates.value(x) - UNIXEPOCH_TS
end

Base.convert{P,Z}(::Type{DateTime}, x::Timestamp{P,Z}) = unix2datetime(P, x)
Base.show(io::IO, x::Timestamp) = show(io, convert(DateTime,x))

immutable Date
    value::Int32
end

const UNIXEPOCH_DT = Dates.value(Dates.Date(1970))
function unix2date(x)
    rata = UNIXEPOCH_DT + x.value
    return Dates.Date(Dates.UTD(rata))
end
date2unix(x::Dates.Date) = Int32(Dates.value(x) - UNIXEPOCH_DT)

Base.convert(::Type{Dates.Date}, x::Arrow.Date) = unix2date(x.value)
Base.show(io::IO, x::Arrow.Date) = show(io, convert(Dates.Date,x))

# Exact Time type with time unit `P`
immutable Time{P}
    value::Int64
end

end # module


text:
using GymSpaces: TupleSpace, Discrete

#=
Algorithmic environments have the following traits in common:

- A 1-d "input tape" or 2-d "input grid" of characters
- A target string which is a deterministic function of the input characters

Agents control a read head that moves over the input tape. Observations consist
of the single character currently under the read head. The read head may fall
off the end of the tape in any direction. When this happens, agents will observe
a special blank character (with index=0) until they get back in bounds.

Actions consist of 3 sub-actions:
    - Direction to move the read head (left or right, plus up and down for 2-d envs)
    - Whether to write to the output tape
    - Which character to write (ignored if the above sub-action is 0)

An episode ends when:
    - The agent writes the full target string to the output tape.
    - The agent writes an incorrect character.
    - The agent runs out the time limit. (Which is fairly conservative.)

Reward schedule:
    write a correct character: +1
    write a wrong character: -.5
    run out the clock: -1
    otherwise: 0

In the beginning, input strings will be fairly short. After an environment has
been consistently solved over some window of episodes, the environment will
increase the average length of generated strings. Typical env specs require
leveling up many times to reach their reward threshold.
=#

abstract type AlgorithmicEnv <: AbstractEnv end
abstract type TapeAlgorithmicEnv <: AlgorithmicEnv end
abstract type GridAlgorithmicEnv <: AlgorithmicEnv end

include("vis/algorithmicenv.jl")

#=
Called between episodes. Update our running record of episode rewards and,
if appropriate, 'level up' minimum input length.
=#
function check_level_up!(env::AlgorithmicEnv)
    push!(env.reward_shortfalls, env.episode_total_reward - length(env.target))
    if length(env.reward_shortfalls) == env.reward_shortfalls.capacity &&
       minimum(env.reward_shortfalls) >= env.MIN_REWARD_SHORTFALL_FOR_PROMOTION &&
       env.min_length < 30
       env.min_length += 1
       env.reward_shortfalls = CircularBuffer{Float32}(10)
   end
end

function reset!(env::AlgorithmicEnv)
    check_level_up!(env)
    env.last_action = nothing
    env.last_reward = 0
    env.read_head_position = env.READ_HEAD_START
    env.write_head_position = 1
    env.episode_total_reward = 0f0
    env.time = 0
    len = rand(0:2) + env.min_length
    generate_input_data!(env, len)
    target_from_input_data!(env)
    _get_obs(env)
end

function step!(env::AlgorithmicEnv, action)
    @assert action ∈ env.action_space "$(action) is not a valid action for this environment."
    env.last_action = action
    inp_act = action[1:1]
    out_act = action[2:2]
    pred = action[3:3]
    done = false
    reward = 0f0
    env.time += 1  # Clock ticks one time
    @assert 1 <= env.write_head_position  # Ensure that we have someplace to write to
    if out_act[1] == 2
        correct = false  # initialising 'correct'
        try
            correct = pred[1] == env.target[env.write_head_position]
        catch BoundsError
            @warn "It looks like you're calling step() even though this " *
            "environement has already returned done=true. You should " *
            "always call reset!() once you recieve done=true. Any " *
            "further steps are undefined behaviour."
            correct = false
        end
        reward += correct ? 1f0 : -5f-1  # Calculate reward.
        !correct && (done = true)  # Bail immediately if wrong character is written.

        env.write_head_position += 1
        env.write_head_position > length(env.target) && (done = true)
    end
    move!(env, inp_act)
    # If an agent takes more than this many timesteps, end the episode
    # immediately and return negative reward
    timelimit = time_limit(env)
    if env.time > timelimit
        reward = 0reward - 1f0
        done = true
    end
    obs = _get_obs(env)
    env.last_reward = reward
    env.episode_total_reward += reward
    return (obs, reward, done, Dict())
end

generate_input_data!(env::TapeAlgorithmicEnv, size_) =
    env.input_data = [rand(1:env.base) for _ in 1:size_]
generate_input_data!(env::GridAlgorithmicEnv, size_) =
    env.input_data = [[rand(1:env.base) for _  in 1:env.rows] for __ in 1:size_]

function move!(env::TapeAlgorithmicEnv, movement)
    named = env.MOVEMENTS[movement[1]]
    env.read_head_position += named == :right ? 1 : -1
end

function move!(env::GridAlgorithmicEnv, movement)
    named = env.MOVEMENTS[movement[1]]
    x, y = env.read_head_posiion
    if named == :left
        x -= 1
    elseif named == :right
        x += 1
    elseif named == :up
        y -= 1
    elseif named == :down
        y += 1
    else
        throw(ArgumentError("Unrecognized direction: $(named)"))
    end
    env.read_head_posiion = [x, y]
end

# return character under read_head/pos. If the read_head/pos is out-of-bounds, then return null(0).
_get_obs(env::AlgorithmicEnv) = _get_obs(env, env.read_head_position)
function _get_obs(env::TapeAlgorithmicEnv, pos::Integer)
    try
        return [env.input_data[pos]]
    catch BoundsError
        return [env.base + 1]
    end
end

function _get_obs(env::GridAlgorithmicEnv, pos::Array)
    x, y = pos
    try
        return [env.input_data[x][y]]
    catch BoundsError
        return [env.base + 1]
    end
end

time_limit(env::TapeAlgorithmicEnv) = length(env.input_data) + length(env.target) + 4


text:
# benchmarks inner loop acquisition functions
using LinearAlgebra
using GaussianDistributions
using SARA
using SARA: random_sampling, uncertainty_sampling, inner_sampling,
        integrated_uncertainty_sampling, integrated_uncertainty, iSARA
using ForwardDiff
using Base.Threads: @threads
using CovarianceFunctions
const Kernel = CovarianceFunctions

# computes the mean derivative of the Gaussian process conditioned on (x, y)
function get_gradient(G::Gaussian, x::AbstractVector, y::AbstractVector, σ::Real)
    C = G | (x, y, σ^2)
    D = GaussianDistributions.derivative(C)
    ∂y = mean(D, x)
end

############### generate sampling paths with different policies
function sampling_path(policy, G::Gaussian, x::AbstractVector, y::AbstractVector,
                       ∂y::AbstractVector, σ::Real, nsample::Int)
    C = G
    indices = zeros(Int, 0)
    rmse, r2 = zeros(nsample, 2), zeros(nsample, 2) # rmse for data and estimated gradients
    sstot = sum(abs2, y .- mean(y)) # total sum of squares
    ∂sstot = sum(abs2, ∂y .- mean(∂y))
    for i in 1:nsample # random sampling
        xi = policy(C, x)
        append!(indices, findfirst(==(xi), x))
        sort!(indices)
        C = G | @views (x[indices], y[indices], σ^2)
        # RMSE (root mean squared error) and R2
        μ = mean(C, x)
        δ = μ - y
        rmse[i, 1] = sqrt(mean(δ.^2))
        r2[i, 1] = 1 - (sum(abs2, δ) / sstot)
        # RMSE and R2 of gradient prediction
        D = GaussianDistributions.derivative(C) # derivative process
        ∂μ = mean(D, x)
        ∂δ = ∂μ - ∂y
        rmse[i, 2] = sqrt(mean(δ.^2))
        r2[i, 2] = 1 - (sum(abs2, ∂δ) / ∂sstot)
    end
    return rmse, r2
end

# on a single stripe
function acquisition_benchmark!(rmse, r2, policies, k, x, y, σ)
    nsample = size(rmse, 1)
    G = Gaussian(k)
    ∂y = get_gradient(G, x, y, σ) # computes gradient "ground truth" on all available data
    for (i, policy) in enumerate(policies)
        println("policy $i of $(length(policies))")
        @time rmse[:, :, i], r2[:, :, i] = sampling_path(policy, G, x, y, ∂y, σ, nsample)
    end
    return rmse, r2
end

function acquisition_benchmark_driver(position, optical, rescaling_parameters,
                                      σ::Real, y_normalization::Real)
    nstripes = length(position)
    # nrandom = 32 # uncomment these two lines to run multiple random samplings to report statistics
    # policies = [random_sampling for _ in 1:nrandom]
    policies = [random_sampling, uncertainty_sampling,
                integrated_uncertainty_sampling(σ), inner_sampling(σ)]
    npolicies = length(policies)
    nsample = 128 # number of samples the active learning algorithm collects
    nmetrics = 2 # number of error metrics to compute
    rmse = zeros(nsample, nmetrics, npolicies, nstripes)
    r2 = zeros(nsample, nmetrics, npolicies, nstripes)
    @threads for i in 1:nstripes
        println("stripe $i of $nstripes")
        θ = rescaling_parameters[i]
        # uncomment these two lines to run the acquisition benchmark with the specialized SARA kernel
        # kernel for inner loop acquisition (l cross-validated from kernel benchmarks)
        l = .2 # is best from kernel benchmarks
        k = iSARA(l, θ)
        # l = .05 # is best from kernel benchmarks for EQ
        # k = Kernel.Lengthscale(Kernel.EQ(), l)
        rmse_i = @view rmse[:, :, :, i]
        r2_i = @view r2[:, :, :, i]
        x = position[i]
        y = optical[i][1] # for simplicity, benchmark with first and most significant optical coefficient
        y = y .- mean(y) # normalized and centered
        y /= y_normalization
        acquisition_benchmark!(rmse_i, r2_i, policies, k, x, y, σ)
    end

    if true
        # f = h5open("inner_random_acquisition_benchmark.h5", "w")
        # f["policy names"] = ["random" for _ in 1:npolicies]
        f = h5open("inner_acquisition_benchmark.h5", "w")
        f["policy names"] = ["random", "uncertainty", "int. unc.", "int. grad. unc."]
        f["rmse"] = rmse
        f["r2"] = r2
        f["nsample"] = nsample
        # f["kernel"] = ["EQ"]
        f["kernel"] = ["iSARA+RGB"]
        f["nstripes"] = nstripes
        f["npolicies"] = npolicies
        f["y_normalization"] = y_normalization
        f["metrics"] = ["rmse", "r2", "rmse of gradient", "r2 of gradient"]
        close(f)
    end
    return rmse, r2
end

include("inner_load_data.jl")
using SARA: estimate_noise

path = "SARA/ScienceAdvances2021/inner/data/"
file = "Bi2O3_19F44_01_inner_loop_data.json"
f = JSON.parsefile(path * file)

position, optical, rescaling_parameters, T_max, log10_τ = load_data(f)

# estimating noise variance of normalized coefficient
x = position[1]
y = optical[1][1]
y = y .- mean(y)
y_normalization = std(y)
y /= y_normalization # std(y)
# θ = rescaling_parameters[1]
# k = iSARA(.2, θ)
l = .05
k = Kernel.Lengthscale(Kernel.EQ(), l)
σ = estimate_noise(x, y, k)
println(σ)
rmse, r2 = acquisition_benchmark_driver(position, optical, rescaling_parameters, σ, y_normalization)

# mrmse = reshape(mean(rmse, dims = (3, 4)), 64, 5)

# random analysis
m = mean(r2, dims = (3, 4))[:, 2]
s = std(r2, dims = (3, 4))[:, 2]

# standard errors
s ./= sqrt(prod(size(r2)[3:4]))
threshold(p) = findfirst(>(p), m)
println(threshold(.8))
println(threshold(.9))
println(s)


text:
using CSV
using DataFrames
using RCall

@rimport readr as readr

# load functions to fit data via least squares
println(snakemake)
include(string(snakemake.scriptdir, "/optim.jl"))

println("Read in data...")
println(snakemake.input["oesophagusmetadata"])
DFdonor = rcopy(readr.read_csv(snakemake.input["oesophagusmetadata"]))
DFcohort = rcopy(readr.read_csv(snakemake.input["oesophagusdnds"]))
DFcohort[:A] = 2 * DFcohort[:areacutoff];

println("Reading in per gene data...")
DF = rcopy(readr.read_csv(snakemake.input["oesophagusdndsgenes"]))
DF[:A] = 2 .* DF[:areacutoff];
println("Data read in")

myDFmiss = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, Float64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :nmutations], 0)

println("Analyse per patient data for missense mutations...")
for p in DFdonor[:patient]
    println("Analysing patient $p")
    DFpatient = filter(row -> row[:patient] == p, DFcohort);
    DFpatient = filter(row -> row[:name] == "wmis", DFpatient)
    if DFpatient[:mle][1] == 0.0
        continue
    end
    age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]
    x = LLoptimizationresults(DFpatient[:mle], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
    x.DF[:patient] = p
    x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
    x.DF[:Age2] = age
    x.DF[:nmutations] = DFpatient[:nmutations]
    append!(myDFmiss, x.DF)
end

myDFnon = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, Float64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :nmutations], 0)

println("Analyse per patient data for nonsense mutations...")
for p in DFdonor[:patient]
    println("Analysing patient $p")
    DFpatient = filter(row -> row[:patient] == p, DFcohort);
    DFpatient = filter(row -> row[:name] == "wnon", DFpatient)
    if DFpatient[:mle][1] == 0.0
        continue
    end
    age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]
    x = LLoptimizationresults(DFpatient[:mle], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
    x.DF[:patient] = p
    x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
    x.DF[:Age2] = age
    x.DF[:nmutations] = DFpatient[:nmutations]
    append!(myDFnon, x.DF)
end

myDFall = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, Float64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :nmutations], 0)

println("Analyse per patient data for all mutations...")
for p in DFdonor[:patient]
    println("Analysing patient $p")
    DFpatient = filter(row -> row[:patient] == p, DFcohort);
    DFpatient = filter(row -> row[:name] == "wall", DFpatient)
    if DFpatient[:mle][1] == 0.0
        continue
    end
    age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]
    x = LLoptimizationresults(DFpatient[:mle], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
    x.DF[:patient] = p
    x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
    x.DF[:Age2] = age
    x.DF[:nmutations] = DFpatient[:nmutations]
    append!(myDFall, x.DF)
end

println("Writing data to file....")
@rput myDFmiss
@rput myDFnon
@rput myDFall
R"""
library(readr)
write_csv(myDFmiss, $(snakemake.output["oesophagusfitmissense"]))
write_csv(myDFnon, $(snakemake.output["oesophagusfitnonsense"]))
write_csv(myDFall, $(snakemake.output["oesophagusfitall"]))
""";

println("Reading in per gene data...")
DF = rcopy(readr.read_csv(snakemake.input["oesophagusdndsgenes"]))
DF[:A] = 2 .* DF[:areacutoff];

delta = Float64[]
rlam = Float64[]
genevec = String[]
mutationtypevec = String[]
rsqvec = Float64[]
patientvec = String[]
myDFmiss = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, String, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :gene, :nmutations], 0)
myDFnon = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, String, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :gene, :nmutations], 0)

println("Analyse per patient per gene data for missense mutations...")
for gene in unique(DF[:gene_name])
#for gene in ["TP53", "NOTCH1"]
    println("Analysing gene: $gene")
    DFgene = filter(row -> row[:gene_name] == gene, DF);
    println("Analysing missense mutations")
    if gene == "SPHKAP"
        println("Errors with this gene, moving to next patient")
        continue
    end
    for p in DFdonor[:patient]
        DFpatient = filter(row -> row[:patient] == p, DFgene);
        age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]

        nmuts = DFpatient[:n_syn][end] + DFpatient[:n_mis][end]
        if nmuts < 5
            println("Patient $p, only $(nmuts) mutations, moving to next patient")
            continue
        end
        println("Analysing patient $p, $(nmuts) mutations")

        try
            global x = LLoptimizationresults(DFpatient[:wmis_cv], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
        catch e
            println("Patient $p failed, moving to next patient")
            continue
        end
        x.DF[:patient] = p
        x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
        x.DF[:Age2] = age
        x.DF[:gene] = gene
        x.DF[:nmutations] = nmuts
        append!(myDFmiss, x.DF)
        push!(delta, x.Δ)
        push!(rlam, x.rλ)
        push!(genevec, gene)
        push!(mutationtypevec, "Missense")
        push!(rsqvec, x.DF[:rsq][1])
        push!(patientvec, p)
    end

    println("Analysing nonsense mutations")
    for p in DFdonor[:patient]
        println("Analysing patient $p")
        DFpatient = filter(row -> row[:patient] == p, DFgene);
        age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]

        nmuts = DFpatient[:n_syn][end] + DFpatient[:n_non][end]
        if nmuts < 5
            println("Only $(nmuts) mutations, moving to next patient")
            continue
        end

        println("Analysing patient $p, $(nmuts) mutations")

        try
            global x = LLoptimizationresults(DFpatient[:wnon_cv], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
        catch e
            println("Patient $p failed, moving to next patient")
            continue
        end
        x.DF[:patient] = p
        x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
        x.DF[:Age2] = age
        x.DF[:gene] = gene
        x.DF[:nmutations] = nmuts
        append!(myDFnon, x.DF)
        push!(delta, x.Δ)
        push!(rlam, x.rλ)
        push!(genevec, gene)
        push!(mutationtypevec, "Nonsense")
        push!(rsqvec, x.DF[:rsq][1])
        push!(patientvec, p)
    end
end

println("Number of fits missense: $(length(unique(myDFmiss[:gene])))")
println("Number of fits nonsense: $(length(unique(myDFnon[:gene])))")

@rput myDFmiss
@rput myDFnon
R"""
library(readr)
write_csv(myDFmiss, $(snakemake.output["oesophagusfitmissensepergene"]))
write_csv(myDFnon, $(snakemake.output["oesophagusfitnonsensepergene"]))
""";

#DF = CSV.File(args["oesophagusdndsneutral"]) |> DataFrame
DF = rcopy(readr.read_csv((snakemake.input["oesophagusdndsneutral"])))
DF[:A] = 2 * DF[:areacutoff];

DFpatient = filter(row -> row[:name] == "wall", DF)
DFpatient = filter(row -> row[:A] .> 0.1, DFpatient)

x = LLoptimizationresults(DFpatient[:mle], DFpatient[:A]; t = 70.0, Amin = 0.14, ρ = 5000.0)
allmutations = x.DF
@rput allmutations
R"""
write_csv(allmutations, $(snakemake.output["oesophagusfitneutral"]))
""";

#########################################################
# Skin data
#########################################################
println("Analysing Skin data")
DFdonor = rcopy(readr.read_csv(snakemake.input["skinmetadata"]));
DFcohort = rcopy(readr.read_csv(snakemake.input["skindnds"]));
DFcohort[:A] = DFcohort[:areacutoff];
patient = unique(DFcohort[:patient])

myDFmiss = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, Float64, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :nmutations], 0)

for p in patient
    println(p)
    DFpatient = filter(row -> row[:patient] == p, DFcohort);
    DFpatient = filter(row -> row[:name] == "wmis", DFpatient)
    if DFpatient[:mle][1] == 0.0
        continue
    end
    age = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
    x = LLoptimizationresults(DFpatient[:mle], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
    x.DF[:patient] = p
    x.DF[:Age] = age
    x.DF[:nmutations] = DFpatient[:nmutations][1]
    append!(myDFmiss, x.DF)
end

myDFnon = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, Float64, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :nmutations], 0)

for p in patient
    println(p)
    DFpatient = filter(row -> row[:patient] == p, DFcohort);
    DFpatient = filter(row -> row[:name] == "wnon", DFpatient)
    if DFpatient[:mle][1] == 0.0
        continue
    end
    age = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
    x = LLoptimizationresults(DFpatient[:mle], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
    x.DF[:patient] = p
    x.DF[:Age] = age
    x.DF[:nmutations] = DFpatient[:nmutations][1]
    append!(myDFnon, x.DF)
end

myDFall = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, Float64, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :nmutations], 0)

for p in patient
    println(p)
    DFpatient = filter(row -> row[:patient] == p, DFcohort);
    DFpatient = filter(row -> row[:name] == "wall", DFpatient)
    if DFpatient[:mle][1] == 0.0
        continue
    end
    age = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
    x = LLoptimizationresults(DFpatient[:mle], DFpatient[:areacutoff]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
    x.DF[:patient] = p
    x.DF[:Age] = age
    x.DF[:nmutations] = DFpatient[:nmutations][1]
    append!(myDFall, x.DF)
end

@rput myDFmiss
@rput myDFnon
R"""
library(readr)
write_csv(myDFmiss, $(snakemake.output["skinfitmissense"]))
write_csv(myDFnon, $(snakemake.output["skinfitnonsense"]))
""";

#DF = CSV.File("data/skin/dnds_patient_genes_combined.csv") |> DataFrame
DF = rcopy(readr.read_csv(snakemake.input["skindndsgenes"], guess_max = 10^5))
DF[:A] = DF[:areacutoff];

delta = Float64[]
rlam = Float64[]
genevec = String[]
mutationtypevec = String[]
rsqvec = Float64[]
patientvec = String[]
myDFmiss= DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, Float64, String, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :gene, :nmutations], 0)
myDFnon = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, Float64, String, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :gene, :nmutations], 0)

for gene in unique(DF[:gene_name])
    println(gene)
    DFgene = filter(row -> row[:gene_name] == gene, DF);
    if gene == "DICER1"
        continue
        println("Skipping DICER1")
    end

    if gene == "KCNH5"
        continue
        println("Skipping KCNH5")
    end

    for p in DFdonor[:patient]
        println(p)
        DFpatient = filter(row -> row[:patient] == p, DFgene);
        age = filter(row -> row[:patient] == p, DFdonor)[:Age][1]

        nmuts = DFpatient[:n_syn][end] + DFpatient[:n_mis][end]
        if nmuts < 5
            println("Patient $p, only $(nmuts) mutations, moving to next patient")
            continue
        end
        println("Analysing patient $p, $(nmuts) mutations")

        try
            global x = LLoptimizationresults(DFpatient[:wmis_cv], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
        catch e
            println("Patient $p failed, moving to next patient")
            continue
        end
        x.DF[:patient] = p
        x.DF[:Age] = age
        x.DF[:gene] = gene
        x.DF[:nmutations] = DFpatient[:n_syn][end] + DFpatient[:n_mis][end] + DFpatient[:n_non][end]
        println(length(x.DF[:A]))
        append!(myDFmiss, x.DF)
        push!(delta, x.Δ)
        push!(rlam, x.rλ)
        push!(genevec, gene)
        push!(mutationtypevec, "Missense")
        push!(rsqvec, x.DF[:rsq][1])
        push!(patientvec, p)
    end

    for p in DFdonor[:patient]
        println(p)
        DFpatient = filter(row -> row[:patient] == p, DFgene);

        nmuts = DFpatient[:n_syn][end] + DFpatient[:n_mis][end]
        if nmuts < 5
            println("Patient $p, only $(nmuts) mutations, moving to next patient")
            continue
        end
        println("Analysing patient $p, $(nmuts) mutations")

        age = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
        try
            global x = LLoptimizationresults(DFpatient[:wnon_cv], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
        catch e
            println("Patient $p failed, moving to next patient")
            continue
        end
        x.DF[:patient] = p
        x.DF[:Age] = age
        x.DF[:gene] = gene
        x.DF[:nmutations] = DFpatient[:n_syn][1] + DFpatient[:n_mis][1] + DFpatient[:n_non][1]
        append!(myDFnon, x.DF)
        push!(delta, x.Δ)
        push!(rlam, x.rλ)
        push!(genevec, gene)
        push!(mutationtypevec, "Nonsense")
        push!(rsqvec, x.DF[:rsq][1])
        push!(patientvec, p)
    end
end

@rput myDFmiss
@rput myDFnon
R"""
write_csv(myDFmiss, $(snakemake.output["skinfitmissensepergene"]))
write_csv(myDFnon, $(snakemake.output["skinfitnonsensepergene"]))
"""


text:
n_values  = [200, 300, 400]
k_values  = [200, 300, 400]
m_values  = [200, 300, 400]

testset_name_suffix = "(main)"

include("_matmul.jl")


text:
function Base.show(io::IO, drv::Driver)
    drv.ptr == C_NULL && (return print(io, "NULL Driver"))
    print(io, "Driver: $(shortname(drv))/$(longname(drv))")
end

function Base.show(io::IO, dataset::AbstractDataset)
    dataset.ptr == C_NULL && (return print(io, "NULL Dataset"))
    println(io, "GDAL Dataset ($(getdriver(dataset)))")
    println(io, "File(s): ")
    for (i,filename) in enumerate(filelist(dataset))
        println(io, "  $filename")
        if i > 5
            println(io, "  ...")
            break
        end
    end
    nrasters = nraster(dataset)
    if nrasters > 0
        print(io, "\nDataset (width x height): ")
        println(io, "$(width(dataset)) x $(height(dataset)) (pixels)")
        println(io, "Number of raster bands: $nrasters")
        for i in 1:min(nrasters, 3)
            print(io, "  ")
            summarize(io, getband(dataset, i))
        end
        nrasters > 3 && println(io, "  ...")
    end

    nlayers = nlayer(dataset)
    if nlayers > 0
        println(io, "\nNumber of feature layers: $nlayers")
        ndisplay = min(nlayers, 5) # display up to 5 layers
        for i in 1:ndisplay
            layer = getlayer(dataset, i-1)
            layergeomtype = getgeomtype(layer)
            println(io, "  Layer $(i-1): $(getname(layer)) ($layergeomtype)")
        end
        if nlayers > 5
            print(io, "  Remaining layers: ")
            for i in 6:nlayers
                print(io, "$(getname(getlayer(dataset, i-1))) ")
                # display up to 5 layer names per line
                if i % 5 == 0 println() end
            end
        end
    end
end

function summarize(io::IO, rasterband::AbstractRasterBand)
    rasterband.ptr == C_NULL && (return print(io, "NULL RasterBand"))
    access = accessflag(rasterband)
    color = getname(getcolorinterp(rasterband))
    xsize = width(rasterband)
    ysize = height(rasterband)
    i = indexof(rasterband)
    pxtype = pixeltype(rasterband)
    println(io, "[$access] Band $i ($color): $xsize x $ysize ($pxtype)")
end

function Base.show(io::IO, rasterband::AbstractRasterBand)
    rasterband.ptr == C_NULL && (return print(io, "NULL RasterBand"))
    summarize(io, rasterband)
    (x,y) = blocksize(rasterband)
    sc = getscale(rasterband)
    ofs = getoffset(rasterband)
    norvw = noverview(rasterband)
    ut = getunittype(rasterband)
    nv = getnodatavalue(rasterband)
    print(io, "    blocksize: $(x)×$(y), nodata: $nv, ")
    println(io, "units: $(sc)px + $ofs$ut")
    print(io, "    overviews: ")
    for i in 1:norvw
        ovr_band = getoverview(rasterband, i-1)
        print(io, "($(i-1)) $(width(ovr_band))x$(height(ovr_band)) ")
        i % 3 == 0 && print(io, "\n               ")
    end
end

# assumes that the layer is reset, and will reset it after display
function Base.show(io::IO, layer::AbstractFeatureLayer)
    layer.ptr == C_NULL && (return println(io, "NULL Layer"))
    layergeomtype = getgeomtype(layer)
    println(io, "Layer: $(getname(layer))")
    featuredefn = layerdefn(layer)

    # Print Geometries
    n = ngeom(featuredefn)
    ngeomdisplay = min(n, 3)
    for i in 1:ngeomdisplay
        gfd = getgeomdefn(featuredefn, i-1)
        display = "  Geometry $(i-1) ($(getname(gfd))): [$(gettype(gfd))]"
        if length(display) > 75
            println(io, "$display[1:70]...")
            continue
        end
        if ngeomdisplay == 1 # only support printing of a single geom column
            for f in layer
                geomwkt = toWKT(getgeom(f))
                length(geomwkt) > 25 && (geomwkt = "$(geomwkt[1:20])...)")
                newdisplay = "$display, $geomwkt"
                if length(newdisplay) > 75
                    display = "$display, ..."
                    break
                else
                    display = newdisplay
                end
            end
        end
        println(io, display)
        resetreading!(layer)
    end
    n > 3 && println(io, "  ...\n  Number of Geometries: $n")

    # Print Features
    n = nfield(featuredefn)
    nfielddisplay = min(n, 5)
    for i in 1:nfielddisplay
        fd = getfielddefn(featuredefn, i-1)
        display = "     Field $(i-1) ($(getname(fd))): [$(gettype(fd))]"
        if length(display) > 75
            println(io, "$display[1:70]...")
            continue
        end
        for f in layer
            field = string(getfield(f, i-1))
            length(field) > 25 && (field = "$(field[1:20])...")
            newdisplay = "$display, $field"
            if length(newdisplay) > 75
                display = "$display, ..."
                break
            else
                display = newdisplay
            end
        end
        println(io, display)
        resetreading!(layer)
    end
    n > 5 && print(io, "...\n Number of Fields: $n")
end

function Base.show(io::IO, featuredefn::AbstractFeatureDefn)
    featuredefn.ptr == C_NULL && (return print(io, "NULL FeatureDefn"))
    n = ngeom(featuredefn)
    ngeomdisplay = min(n, 3)
    for i in 1:ngeomdisplay
        gfd = getgeomdefn(featuredefn, i-1)
        println(io, "  Geometry (index $(i-1)): $gfd")
    end
    n > 3 && println(io, "  ...\n  Number of Geometries: $n")

    n = nfield(featuredefn)
    nfielddisplay = min(n, 5)
    for i in 1:nfielddisplay
        fd = getfielddefn(featuredefn, i-1)
        println(io, "     Field (index $(i-1)): $fd")
    end
    n > 5 && print(io, "...\n Number of Fields: $n")
end

function Base.show(io::IO, fd::AbstractFieldDefn)
    fd.ptr == C_NULL && (return print(io, "NULL FieldDefn"))
    print(io, "$(getname(fd)) ($(gettype(fd)))")
end

function Base.show(io::IO, gfd::AbstractGeomFieldDefn)
    gfd.ptr == C_NULL && (return print(io, "NULL GeomFieldDefn"))
    print(io, "$(getname(gfd)) ($(gettype(gfd)))")
end

function Base.show(io::IO, feature::Feature)
    feature.ptr == C_NULL && (return println(io, "NULL Feature"))
    println(io, "Feature")
    n = ngeom(feature)
    for i in 1:min(n, 3)
        displayname = geomname(getgeom(feature, i-1))
        println(io, "  (index $(i-1)) geom => $displayname")
    end
    n > 3 && println(io, "...\n Number of geometries: $n")
    n = nfield(feature)
    for i in 1:min(n, 10)
        displayname = getname(getfielddefn(feature, i-1))
        print(io, "  (index $(i-1)) $displayname => ")
        println(io, "$(getfield(feature, i-1))")
    end
    n > 10 && print(io, "...\n Number of Fields: $n")
end

function Base.show(io::IO, spref::AbstractSpatialRef)
    spref.ptr == C_NULL && (return print(io, "NULL Spatial Reference System"))
    projstr = toPROJ4(spref)
    if length(projstr) > 45
        projstart = projstr[1:35]
        projend = projstr[end-4:end]
        print(io, "Spatial Reference System: $projstart ... $projend")
    else
        print(io, "Spatial Reference System: $projstr")
    end
end

function Base.show(io::IO, geom::AbstractGeometry)
    geom.ptr == C_NULL && (return print(io, "NULL Geometry"))
    print(io, "Geometry: ")
    geomwkt = toWKT(geom)
    if length(geomwkt) > 60
        print(io, "$(geomwkt[1:50]) ... $(geomwkt[end-4:end])")
    else
        print(io, "$geomwkt")
    end
end


text:
using ReachabilityAnalysis: _jacobian_function,
                            _jacobian_expression,
                            _hessian_function,
                            _hessian_expression

@testset "Symbolic-Numeric utils: Jacobian" begin
    # Brusselator model as a vector of symbolic expressions
    var = @variables x y
    f = [1 + x^2*y - 2.5x, 1.5x - x^2*y]
    Jexpr0 = [-2.5 + 2x*y x^2; 1.5-2x*y -x^2]

    Jexpr = _jacobian_expression(f, var)
    @test isequal(Jexpr, Jexpr0)
    J = _jacobian_function(f, var)
    v = [1.0, 2.0]
    Jv = [1.5 1; -2.5 -1]
    @test J(v) == Jv

    # same using static arrays
    f_st = SVector{2 ,Num}(f)
    Jexpr_st = _jacobian_function(f_st, var)
    J_st = _jacobian_function(f_st, var)
    v_st = SA[1.0, 2.0]
    Jv_st = SMatrix{2, 2, Float64}(Jv)
    aux = J_st(v_st)
    @test isa(aux, StaticArray) && aux == Jv_st

    # same using MathematicalSystem types
    sys = @system(u' = f(u), dim=2)
    J = _jacobian_function(sys, var)
    @test J(v) == Jv
    Jexpr = _jacobian_expression(sys, var)
    @test isequal(Jexpr0, Jexpr)

    U0 = VPolygon([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
    ivp = @ivp(u' = f(u), dim=2, u(0) ∈ U0)
    J = _jacobian_function(ivp, var)
    @test J(v) == Jv
    Jexpr = _jacobian_expression(ivp, var)
    @test isequal(Jexpr0, Jexpr)

    # test in-place version
    prob_iip = bruss(0.1)
    J = _jacobian_function(prob_iip, var)
    @test J(v) == Jv
    Jexpr = _jacobian_expression(prob_iip, var)
    @test isequal(Jexpr0, Jexpr)

    # linear systems
    A = [1 3; 2 -2.]
    sys = @system(x' = Ax)
    J = _jacobian_function(sys, var)
    aux = [1.0, 1.0] # can be any vector
    @test J(aux) == A
    Jexpr = _jacobian_expression(sys, var)
    @test isequal(Jexpr, A)
end

@testset "Symbolic-Numeric utils: Hessian" begin

end


text:
module RunLengthArrays

export RunLengthArray, runs, values

import Base

@doc raw"""
An array of objects of type `T` encoded using run-length compression.

This works like an array, but it is designed to be extremely efficient when the
array is long and consists of long repetitions («runs») of the same elements.
The length of each run is held using type `N`; this type should be chosen
according to the expected maximum length of a run.

Once a `RunLengthArray` has been created, new elements can be added only at the
end of it, using `push!` (add one element) or `append!` (add a sequence of
elements).

```julia
arr = RunLengthArray{Int, Int8}(Int8[3, 3, 3, 7, 7, 7, 7, 7, 7, 4])

longarr = collect(arr)

@assert length(arr) == 10
@assert sum(arr) == 30
@assert arr[1] == 3
@assert arr[2] == 3

println("$(sizeof(arr))")  # Prints 16 (bytes)
println("$(sizeof(longarr))")  # Prints 56 (bytes)
```

"""
mutable struct RunLengthArray{N <: Number,T} <: AbstractArray{T,1}
    runs::Array{N}
    values::Array{T}

    RunLengthArray{N,T}() where {N <: Number,T} = new(N[], T[])
    RunLengthArray{N,T}(runs, values) where {N <: Number,T} = new(collect(runs), collect(values))
end

@doc raw"""
    runs(arr::RunLengthArray{N,T})

Returns a `Array{N,1}` type containing the length of each run in the run-length
array.

```julia
x = RunLengthArray{Int,String}(["X", "X", "X", "O", "O"])
@assert runs(x) == Int[3, 2]
```

To get a list of the values, use `values`.
"""
runs(arr::RunLengthArray{N,T}) where {N,T} = arr.runs

Base.values(arr::RunLengthArray{N,T}) where {N,T} = arr.values

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

function RunLengthArray{N,T}(arr) where {N,T}
    isempty(arr) && return RunLengthArray{N,T}()

    runs = N[]
    values = T[]

    curval = arr[1]
    curcount = 1
    for idx = 2:length(arr)
        if arr[idx] == curval
            curcount += 1
        else
            push!(runs, curcount)
            push!(values, curval)
            curval = arr[idx]
            curcount = 1
        end
    end

    push!(runs, curcount)
    push!(values, curval)

    RunLengthArray{N,T}(runs, values)
end

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

function Base.show(io::IO, arr::RunLengthArray{N,T}) where {N,T}
    print(io, "RunLengthArray{$N,$T}([")
    for elem in eachindex(arr.runs)
        print(io, "$(arr.values[elem])×$(arr.runs[elem])")
        if elem < length(arr.runs)
            print(io, ", ")
        end
    end
    println(io, "])")
end

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

Base.length(arr::RunLengthArray{N,T}) where {N,T} = sum(arr.runs)
Base.size(arr::RunLengthArray{N,T}) where {N,T} = (length(arr),)
Base.IndexStyle(arr::RunLengthArray{N,T}) where {N,T} = IndexLinear()
Base.firstindex(arr::RunLengthArray{N,T}) where {N,T} = 1
Base.lastindex(arr::RunLengthArray{N,T}) where {N,T} = length(arr)

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

function Base.iterate(iter::RunLengthArray{N,T}) where {N,T}
    isempty(iter.runs) && return nothing

    (iter.values[1], (iter.runs[1] - 1, 1))
end

function Base.iterate(iter::RunLengthArray{N,T}, state) where {N,T}
    runsleft, curidx = state

    if runsleft > 0
        runsleft -= 1
        return (iter.values[curidx], (runsleft, curidx))
    else
        curidx == length(iter.runs) && return nothing

        curidx += 1
        runsleft = iter.runs[curidx] - 1
        return (iter.values[curidx], (runsleft, curidx))
    end
end

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

function Base.getindex(arr::RunLengthArray{N,T}, target_idx::Number) where {N,T}
    1 <= target_idx <= length(arr) || throw(BoundsError(arr, target_idx))

    runidx = 1
    elements_left = target_idx
    while elements_left > 0
        elements_left <= arr.runs[runidx] && return arr.values[runidx]
        elements_left -= arr.runs[runidx]
        runidx += 1
    end

    return arr.values[end]
end

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

function Base.sum(arr::RunLengthArray{N,T}) where {N,T}
    result = zero(T)

    for idx in eachindex(arr.runs)
        result += arr.runs[idx] * arr.values[idx]
    end

    result
end

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

Base.minimum(arr::RunLengthArray{N,T}) where {N,T} = minimum(arr.values)
Base.maximum(arr::RunLengthArray{N,T}) where {N,T} = maximum(arr.values)
Base.extrema(arr::RunLengthArray{N,T}) where {N,T} = extrema(arr.values)

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

function Base.sort!(arr::RunLengthArray{N,T}; kwargs...) where {N,T}
    perm = sortperm(arr.values, kwargs...)

    arr.runs = arr.runs[perm]
    arr.values = arr.values[perm]
end

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

function Base.sort(arr::RunLengthArray{N,T}; kwargs...) where {N,T}
    perm = sortperm(arr.values, kwargs...)

    RunLengthArray{N,T}(arr.runs[perm], arr.values[perm])
end

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

function Base.push!(arr::RunLengthArray{N,T}, newval::T) where {N,T}
    if newval == arr.values[end]
        arr.runs[end] += 1
    else
        push!(arr.runs, 1)
        push!(arr.values, newval)
    end
end

function Base.push!(arr::RunLengthArray{N,T}, newrun::Tuple{N,T}) where {N,T}
    push!(arr.runs, newrun[1])
    push!(arr.values, newrun[2])
end

@doc raw"""
    push!(arr::RunLengthArray{N,T}, newval::T) where {N <: Number,T}
    push!(arr::RunLengthArray{N,T}, newrun::Tuple{N,T}) where {N <: Number,T}

Append a element or a run of elements to the end of a `RunLengthArray`. In the
second case, the parameter `newrun` must be a tuple containing the number of
repetitions `n` and the value `v`: `(n, v)`.

```julia
arr = RunLengthArray{Int, Float64}([3, 2], [1.1, 6.5])
@assert length(arr) == 5

# Add 6 instances of the number "1.4" at the end of the array
push!(arr, (6, 1.4))

@assert length(arr) == 11
Base.push!
```

`"""
push!

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

@doc raw"""
    append!(arr::RunLengthArray{N,T}, otherarray)

Append an array to a `RunLengthArray`, modifying it in place.

"""
function Base.append!(arr::RunLengthArray{N,T}, otherarray) where {N,T}
    for elem in otherarray
        push!(arr, T(elem))
    end
end

end # module


text:
using Pkg
Pkg.activate(".")
verbose = true
using Revise
if verbose   println("# Loading RvSpecML")    end
using RvSpectML
if verbose   println("# Loading other packages")    end
using DataFrames, Query

# TODO: USER:  Either the paths that specify where datafiles are stored here or in examples/data_paths.jl
if verbose   println("# Creating manifest of files to process.")    end
if isdir("examples")
   cd("examples")
end
solar_data_path = "20190918"
ancilary_solar_data_path = "."
 if isfile("data_paths.jl")
    include("data_paths.jl")
 end
 solar_data_path

num_spectra_to_use = 20
 bjd_first_good = 2458745.1296134139
 bjd_last_good = 2458745.283
 df_files = NEID.make_manifest(solar_data_path)

df_files_use = df_files |>
  @filter( _.target == "Solar" ) |>
  @filter(bjd_first_good <= _.bjd < bjd_last_good) |>
#  @take(num_spectra_to_use) |>
  DataFrame

if verbose println("# Reading in FITS files.")  end
@time solar_data = map(NEID.read_solar_data,eachrow(df_files_use))

if verbose println("# Applying wavelength corrections.")  end
NEID.read_drift_corrections!(joinpath(ancilary_solar_data_path,"SolarRV20190918_JD_SciRV_CalRV.txt"), df_files_use)

NEID.read_barycentric_corrections!(joinpath(ancilary_solar_data_path,"SolarTelescope2019-09-18_inclGravRedshiftAirMassAltitude.csv"), df_files_use)

NEID.read_differential_extinctions!(joinpath(ancilary_solar_data_path,"20190918_diff_ex_full_fixed.txt"), df_files_use)

apply_doppler_boost!(solar_data,df_files_use)


text:
"""
      Simplified version of LowRankEigenAlign to only make use of the aligned
    edges, rather than considering non-matches, and punishing mis-matches.
"""
function lowRankEigenAlignEdgesOnly(A::SparseMatrixCSC{T, Int64}, B::SparseMatrixCSC{T, Int64}) where T

    m = size(A,1)
    n = size(B,1)

    tol=1e-8
    maxIter = 50 #likely to converge to desired tolerance

    v,_ = powerMethod(A, ones(m), tol, maxIter)
    u,_ = powerMethod(B, ones(n), tol, maxIter)
    mapping, _ = rank_one_matching(v,u)

    return mapping

end

#used for the LowRankEigenAlignEdgesOnly routine
function powerMethod(A::SparseMatrixCSC{T, Int64},x0::Array{T,1},tol::Float64,maxIter::Int;update_user =-1) where T

    y_k = copy(x0)
    y_k ./=norm(y_k)

    λ_k = Inf

    for i =1:maxIter

        y_k_1 = A*y_k
        λ_k_1 = y_k_1'*y_k

        y_k_1 ./=norm(y_k_1)

        if update_user != -1 && i % update_user == 0
            println("step $i : λ_k = $λ_k_1")
        end

        if abs(λ_k - λ_k_1) < tol || i == maxIter
            return y_k_1, λ_k_1
        else
            λ_k = λ_k_1
            y_k = y_k_1
        end

    end

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
using Manifolds, Manopt, Test, ManifoldsBase

@testset "Record Options" begin
    # helper to get debug as string
    io = IOBuffer()
    M = ManifoldsBase.DefaultManifold(2)
    x = [4.,2.]
    o = GradientDescentOptions(x, StopAfterIteration(20), ConstantStepsize(1.))
    f = y -> distance(M,y,x).^2
    ∇f = y -> -2*log(M,y,x)
    p = GradientProblem(M,f,∇f)
    a = RecordIteration()
    # constructors
    rO = RecordOptions(o,a)
    @test Manopt.dispatch_options_decorator(rO) === Val{true}()
    @test get_options(o) == o
    @test get_options(rO) == o
    @test_throws MethodError get_options(p)
    #
    @test get_initial_stepsize(p,rO) == 1.
    @test get_stepsize(p,rO,1) == 1.
    @test get_last_stepsize(p,rO,1) == 1.
    #
    @test rO.recordDictionary[:All] == a
    @test RecordOptions(o,[a]).recordDictionary[:All].group[1] == a
    @test RecordOptions(o,Dict(:A => a)).recordDictionary[:A] == a
    @test isa( RecordOptions(o, [:Iteration]).recordDictionary[:All].group[1], RecordIteration)
    @test !has_record(o)
    @test_throws ErrorException get_record(o)
    @test get_options(o) == o
    @test !has_record(DebugOptions(o,[]))
    @test has_record(rO)
    @test_throws ErrorException get_record(o)
    @test length( get_record(rO,:All) ) == 0
    @test length( get_record(rO) ) == 0
    @test length(  get_record(DebugOptions( rO,[]) )  ) == 0
    @test_throws ErrorException get_record(RecordOptions(o,Dict{Symbol,RecordAction}()))
    @test_throws ErrorException get_record(o)
    @test get_record(rO) == Array{Int64,1}()
    # RecordIteration
    @test a(p,o,0) == nothing # inactive
    @test a(p,o,1) == [1]
    @test a(p,o,2) == [1,2]
    @test a(p,o,9) == [1,2,9]
    @test a(p,o,-1) == []
    # RecordGroup
    @test length(RecordGroup().group) == 0
    b = RecordGroup([ RecordIteration(), RecordIteration() ])
    b(p,o,1)
    b(p,o,2)
    @test b.group[1].recordedValues == [1,2]
    @test b.group[2].recordedValues == [1,2]
    @test get_record(b) == [ (1,1), (2,2) ]
    #RecordEvery
    c = RecordEvery(a,10,true)
    @test c(p,o,0) == nothing
    @test c(p,o,1) == nothing
    @test c(p,o,10) == [10]
    @test c(p,o,20) == [10,20]
    @test c(p,o,-1) == []
    @test get_record(c) == []
    # RecordChange
    c = RecordChange()
    c(p,o,1)
    @test c.recordedValues == [0.] # no x0 -> assume x0 is the first iterate
    o.x = [3.,2.]
    c(p,o,2)
    @test c.recordedValues == [0.,1.] # no x0 -> assume x0 is the first iterate
    c = RecordChange([4.,2.])
    c(p,o,1)
    @test c.recordedValues == [1.] # no x0 -> assume x0 is the first iterate
    # RecordEntry
    o.x = x
    d = RecordEntry(x,:x)
    d(p,o,1)
    @test d.recordedValues == [x]
    d2 = RecordEntry(typeof(x),:x)
    d2(p,o,1)
    @test d2.recordedValues == [x]
    # RecordEntryChange
    o.x = x
    e = RecordEntryChange(:x, (p,o,x,y) -> distance(p.M,x,y))
    @test update_storage!(e.storage,o) == (:x,)
    e2 = RecordEntryChange(x,:x,(p,o,x,y) -> distance(p.M,x,y))
    @test e.field == e2.field
    e(p,o,1)
    @test e.recordedValues == [0.]
    o.x = [3.0,2.0]
    e(p,o,2)
    @test e.recordedValues == [0.,1.]
    # RecordIterate
    o.x = x
    f = RecordIterate(x)
    @test_throws ErrorException RecordIterate()
    f(p,o,1)
    @test f.recordedValues == [x]
    # RecordCost
    g = RecordCost()
    g(p,o,1)
    @test g.recordedValues == [0.]
    o.x = [3.,2.]
    g(p,o,2)
    @test g.recordedValues == [0.,1.]
    #RecordFactory
    o.∇ = [0.,0.]
    rf = RecordFactory(o,[:Cost,:∇])
    @test isa(rf[:All], RecordGroup)
    @test isa(rf[:All].group[1], RecordCost)
    @test isa(rf[:All].group[2], RecordEntry)
    @test isa( RecordFactory(o,[2])[:All], RecordEvery)
    @test rf[:All].group[2].field == :∇
    @test length(rf[:All].group) == 2
    @test all(  isa.( RecordFactory(o,[:Cost, :Iteration, :Change, :Iterate])[:All].group,
         [RecordCost, RecordIteration, RecordChange, RecordIterate]) )
    @test RecordActionFactory(o,g) == g
end


text:
using EchogramUtils
using Test

@test pow2db(1000) == 30
@test db2pow(30) == 1000

@test mag2db(10) == 20
@test db2mag(20) == 10


text:
using DynamicGrids, Test, Dates
import DynamicGrids: Remove, NoOpt

life = Life()
@test occursin("Life{:_default_,:_default_}", sprint((io, s) -> show(io, MIME"text/plain"(), s), life))

rs = Ruleset(;
    rules=(Life(),),
    timestep=Day(1),
    boundary=Remove(),
    opt=NoOpt(),
)
@test occursin("Ruleset =", sprint((io, s) -> show(io, MIME"text/plain"(), s), rs))
@test occursin("Life{:_default_,:_default_}", sprint((io, s) -> show(io, MIME"text/plain"(), s), rs))
@test occursin(r"opt = .*NoOpt()", "nopt = DynamicGrids.NoOpt()")
@test occursin(r"boundary = .*Remove()", sprint((io, s) -> show(io, MIME"text/plain"(), s), rs))

rule1 = Cell{:a,:b}() do a
    2a
end
@test occursin("Cell{:a,:b}", sprint((io, s) -> show(io, MIME"text/plain"(), s), rule1))

rule2 = Cell{Tuple{:b,:d},:c}() do b, d
    b + d
end
@test occursin("Cell{Tuple{:b,:d},:c}", sprint((io, s) -> show(io, MIME"text/plain"(), s), rule2))

chain = Chain(rule1, rule2)

@test occursin("Chain{Tuple{:a,:b,:d},Tuple{:b,:c}}", sprint((io, s) -> show(io, MIME"text/plain"(), s), chain))
@test occursin("    Cell{:a,:b}", sprint((io, s) -> show(io, MIME"text/plain"(), s), chain))
@test occursin("    Cell{Tuple{:b,:d},:c}", sprint((io, s) -> show(io, MIME"text/plain"(), s), chain))


text:
export CuTwoDGrid

struct CuTwoDGrid{T} <: AbstractTwoDGrid
  nx::Int
  ny::Int
  nk::Int
  nl::Int
  nkr::Int
  Lx::T
  Ly::T
  dx::T
  dy::T
  x::CuArray{T,2}
  y::CuArray{T,2}
  X::CuArray{T,2}
  Y::CuArray{T,2}
  k::CuArray{T,2}
  l::CuArray{T,2}
  kr::CuArray{T,2}
  K::CuArray{T,2}
  L::CuArray{T,2}
  Kr::CuArray{T,2}
  Lr::CuArray{T,2}
  KKsq::CuArray{T,2}      # K^2 + L^2
  invKKsq::CuArray{T,2}   # 1/KKsq, invKKsq[1, 1]=0
  KKrsq::CuArray{T,2}     # Kr^2 + Lr^2
  invKKrsq::CuArray{T,2}  # 1/KKrsq, invKKrsq[1, 1]=0

  fftplan::CuArrays.FFT.cCuFFTPlan{Complex{T},-1,false,2}
  rfftplan::CuArrays.FFT.rCuFFTPlan{T,-1,false,2}

  # Range objects that access the aliased part of the wavenumber range
  ialias::UnitRange{Int}
  iralias::UnitRange{Int}
  jalias::UnitRange{Int}
end

"""
    CuTwoDGrid(nx, Lx, ny=nx, Ly=Lx;  x0=-Lx/2, y0=-Ly/2)

Construct a CuTwoDGrid object. The two-dimensional domain has size `(Lx, Ly)`,
resolution `(nx, ny)` and bottom left corner at `(x0, y0)`. The type of `CuTwoDGrid` is
inferred from the type of `Lx`.
"""
function CuTwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-0.5*Lx, y0=-0.5*Ly)
  T = typeof(Lx)
  dx = Lx/nx
  dy = Ly/ny
  nk = nx
  nl = ny
  nkr = Int(nx/2+1)

  # Physical grid
  x = Array{T}(reshape(linspace(x0, x0+Lx-dx, nx), (nx, 1)))
  y = Array{T}(reshape(linspace(y0, y0+Ly-dy, ny), (1, ny)))
  X = [ x[i] for i = 1:nx, j = 1:ny]
  Y = [ y[j] for i = 1:nx, j = 1:ny]

  # Wavenubmer grid
  i1 = 0:Int(nx/2)
  i2 = Int(-nx/2+1):-1
  j1 = 0:Int(ny/2)
  j2 = Int(-ny/2+1):-1

  k  = reshape(2π/Lx*cat(1, i1, i2), (nk, 1))
  kr = reshape(2π/Lx*cat(1, i1), (nkr, 1))
  l  = reshape(2π/Ly*cat(1, j1, j2), (1, nl))

  K = [ k[i] for i = 1:nk, j = 1:nl]
  L = [ l[j] for i = 1:nk, j = 1:nl]
  Kr = [ kr[i] for i = 1:nkr, j = 1:nl]
  Lr = [ l[j]  for i = 1:nkr, j = 1:nl]

  KKsq  = @. K^2 + L^2
  invKKsq = 1./KKsq
  invKKsq[1, 1] = 0

  KKrsq = @. Kr^2 + Lr^2
  invKKrsq = 1./KKrsq
  invKKrsq[1, 1] = 0

  # Convert Arrays to CuArrays
  @cuconvertarrays x y X Y k kr l K L Kr Lr KKsq invKKsq KKrsq invKKrsq

  # FFT plans
    fftplan = plan_fft(CuArray{Complex{T},2}(nx, ny))
   rfftplan = plan_rfft(CuArray{T,2}(nx, ny))

  # Index endpoints for aliased i, j wavenumbers
  iaL, iaR = Int(floor(nk/3))+1, 2*Int(ceil(nk/3))-1
  jaL, jaR = Int(floor(nl/3))+1, 2*Int(ceil(nl/3))-1

  ialias  = iaL:iaR
  iralias = iaL:nkr
  jalias  = iaL:iaR

  CuTwoDGrid(nx, ny, nk, nl, nkr, Lx, Ly, dx, dy, x, y, X, Y,
             k, l, kr, K, L, Kr, Lr, KKsq, invKKsq, KKrsq, invKKrsq,
             fftplan, rfftplan, ialias, iralias, jalias)
end

"""
    CuTwoDGrid(T, args...; kwargs...)

Construct a CuTwoDGrid object with type T. The other args and kwargs are identical
to the constructor for CuTwoDGrid with no explicit type.
"""
CuTwoDGrid(T::DataType, nx, Lx, args...; kwargs...) = CuTwoDGrid(nx, T(Lx), args...; kwargs...)

"""
    CuTwoDGrid(g)

Construct a CuTwoDGrid object with the same type and size as the TwoDGrid object `g`.
"""
CuTwoDGrid(g::TwoDGrid) = CuTwoDGrid(g.nx, g.Lx, g.ny, g.Ly; x0=g.x[1], y0=g.y[1])
TwoDGrid(g::CuTwoDGrid) = TwoDGrid(g.nx, g.Lx, g.ny, g.Ly; x0=g.x[1], y0=g.y[1])

makefilter(g::CuTwoDGrid, T, sz; kwargs...) = CuArray(ones(T, sz)).*makefilter(g; realvars=sz[1]==g.nkr, kwargs...)


text:
"""
get_efficiency(nPassed, nTotal)
Description of ```get_efficiency```
------------------------------
Returns efficiency as ```ε = Nₚ/Nₜₒₜ ```
"""
function get_efficiency(nPassed, nTotal = 1e5)
    return nPassed/nTotal
end


text:
using Documenter
using GPMethylation

makedocs(
    sitename = "GPMethylation",
    format = Documenter.HTML(),
    modules = [GPMethylation]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "https://github.com/owensnick/GPMethylation.jl.git",
)


text:
function foo()
    a = 10
    b = a
    @assert(b == 10)
end

function fibonacci(n)::Int64
    if n == 0
        return 0
    elseif n == 1
        return 1
    else
        return fibonacci(n - 2) + fibonacci(n - 1)
    end
end

function mul_list()
    a::Vector = []
    for i = 0:4
        push!(a, i)
    end
    return repeat(a, 2)
end

function combinations(array)::Vector
    result = []
    for x in array
        for y in array
            push!(result, (x, y))
        end
    end
    return result
end

function mul_recvd_list(a::Vector)
    for i = 0:length(a)-1
        push!(a, i)
    end
    return repeat(a, 2)
end

function plus_test(x, y)::Any
    return x + y
end

function plus_test(x::String, y::String)::String
    return x * y
end

if abspath(PROGRAM_FILE) == @__FILE__
    foo()
    @assert(fibonacci(10) == 55)
    @assert((repeat("test", fibonacci(3))) == "testtest")
    a = []
    a_mul = mul_recvd_list(a)
    @assert(a_mul == [])
    x = "ss"
    y = "zz"
    res = plus_test(x, y)
    @assert(res == "sszz")
    println("OK")
end


text:
@inline function initialize!(integrator,cache::Rosenbrock23Cache,f=integrator.f)
  integrator.kshortsize = 2
  @unpack k₁,k₂,fsalfirst,fsallast = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = fsallast
  integrator.k = [k₁,k₂]
  f(integrator.t,integrator.uprev,integrator.fsalfirst)
end

@inline function perform_step!(integrator,cache::Rosenbrock23Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,uf,tf,linsolve_tmp,jac_config = cache
  jidx = eachindex(J)
  @unpack c₃₂,d = cache.tab

  # Setup Jacobian Calc
  sizeu  = size(u)
  tf.vf.sizeu = sizeu
  tf.uprev = uprev
  uf.vfr.sizeu = sizeu
  uf.t = t

  #if alg_autodiff(alg)
    ForwardDiff.derivative!(dT,tf,t) # Time derivative of each component
    ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
  #else
  #  Calculus.finite_difference!((t)->vecfreturn(t,uprev,du2),[t],dT)
  #  Calculus.finite_difference_jacobian!((du1,uprev)->vecf(t,uprev,du1),vec(uprev),vec(du1),J)
  #end
  a = -dt*d
  for i in 1:length(u), j in 1:length(u)
    W[i,j] = @muladd I[i,j]+a*J[i,j]
  end

  Wfact = integrator.alg.factorization(W)
  a = -a
  for i in uidx
    linsolve_tmp[i] = @muladd fsalfirst[i] + a*dT[i]
  end

  @into! vectmp = Wfact\linsolve_tmp
  recursivecopy!(k₁,reshape(vectmp,size(u)...))
  dto2 = dt/2
  for i in uidx
    u[i]= @muladd uprev[i]+dto2*k₁[i]
  end
  f(t+dto2,u,f₁)

  for i in uidx
    linsolve_tmp[i] = f₁[i]-k₁[i]
  end

  @into! vectmp2 = Wfact\linsolve_tmp
  for i in uidx
    k₂[i] = vectmp2[i] + k₁[i]
    u[i] = @muladd uprev[i] + dt*k₂[i]
  end
  if integrator.opts.adaptive
    f(t+dt,u,integrator.fsallast)

    for i in uidx
      linsolve_tmp[i] = @muladd integrator.fsallast[i] - c₃₂*(k₂[i]-f₁[i])-2(k₁[i]-fsalfirst[i])+dt*dT[i]
    end

    @into! vectmp3 = Wfact\linsolve_tmp
    k₃ = reshape(vectmp3,sizeu...)
    for i in uidx
      tmp[i] = (dt*(k₁[i] - 2k₂[i] + k₃[i])/6)./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
    end
    integrator.EEst = integrator.opts.internalnorm(tmp)
  end
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::Rosenbrock32Cache,f=integrator.f)
  integrator.kshortsize = 2
  @unpack k₁,k₂,fsalfirst,fsallast = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = fsallast
  integrator.k = [k₁,k₂]
  f(integrator.t,integrator.uprev,integrator.fsalfirst)
end

@inline function perform_step!(integrator,cache::Rosenbrock32Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,uf,tf,linsolve_tmp,jac_config = cache
  jidx = eachindex(J)
  @unpack c₃₂,d = cache.tab
  # Setup Jacobian Calc
  sizeu  = size(u)
  tf.vf.sizeu = sizeu
  tf.uprev = uprev
  uf.vfr.sizeu = sizeu
  uf.t = t

  #if alg_autodiff(alg)
    ForwardDiff.derivative!(dT,tf,t) # Time derivative of each component
    ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)

  a = -dt*d
  for i in 1:length(u), j in 1:length(u)
    W[i,j] = @muladd I[i,j]+a*J[i,j]
  end

  Wfact = integrator.alg.factorization(W)

  a = -a
  for i in uidx
    linsolve_tmp[i] = @muladd fsalfirst[i] + a*dT[i]
  end

  @into! vectmp = Wfact\linsolve_tmp

  recursivecopy!(k₁,reshape(vectmp,sizeu...))

  dto2 = dt/2
  for i in uidx
    u[i]= @muladd uprev[i]+dto2*k₁[i]
  end
  f(t+dto2,u,f₁)
  for i in uidx
    linsolve_tmp[i] = f₁[i]-k₁[i]
  end
  @into! vectmp2 = Wfact\linsolve_tmp
  tmp = reshape(vectmp2,sizeu...)
  for i in uidx
    k₂[i] = vectmp2[i] + k₁[i]
  end
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*k₂[i]
  end
  f(t+dt,tmp,integrator.fsallast)
  for i in uidx
    linsolve_tmp[i] = @muladd integrator.fsallast[i] - c₃₂*(k₂[i]-f₁[i])-2(k₁[i]-fsalfirst[i])+dt*dT[i]
  end
  @into! vectmp3 = Wfact\linsolve_tmp
  k₃ = reshape(vectmp3,sizeu...)
  for i in uidx
    u[i] = uprev[i] + dt*(k₁[i] + 4k₂[i] + k₃[i])/6
  end
  if integrator.opts.adaptive
    for i in uidx
      tmp[i] = (dt*(k₁[i] - 2k₂[i] + k₃[i])/6)./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
    end
    integrator.EEst = integrator.opts.internalnorm(tmp)
  end
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::Rosenbrock23ConstantCache,f=integrator.f)
  integrator.kshortsize = 2
  k = eltype(integrator.sol.k)(2)
  integrator.k = k
  integrator.fsalfirst = f(integrator.t,integrator.uprev)
end

@inline function perform_step!(integrator,cache::Rosenbrock23ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack c₃₂,d,tf,uf = cache
  # Time derivative
  tf.u = uprev
  uf.t = t
  dT = ForwardDiff.derivative(tf,t)
  J = ForwardDiff.derivative(uf,uprev)
  a = -dt*d
  W = @muladd 1+a*J
  a = -a
  k₁ = W\@muladd(integrator.fsalfirst + a*dT)
  dto2 = dt/2
  f₁ = f(t+dto2,@muladd uprev+dto2*k₁)
  k₂ = W\(f₁-k₁) + k₁
  u = @muladd uprev + dt*k₂
  if integrator.opts.adaptive
    integrator.fsallast = f(t+dt,u)
    k₃ = W\@muladd(integrator.fsallast - c₃₂*(k₂-f₁)-2(k₁-integrator.fsalfirst)+dt*dT)
    integrator.EEst = integrator.opts.internalnorm((dt*(k₁ - 2k₂ + k₃)/6)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
  end
  integrator.k[1] = k₁
  integrator.k[2] = k₂
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::Rosenbrock32ConstantCache,f=integrator.f)
  integrator.kshortsize = 2
  k = eltype(integrator.sol.k)(2)
  integrator.k = k
  integrator.fsalfirst = f(integrator.t,integrator.uprev)
end

@inline function perform_step!(integrator,cache::Rosenbrock32ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack c₃₂,d,tf,uf = cache
  tf.u = uprev
  uf.t = t
  # Time derivative
  dT = ForwardDiff.derivative(tf,t)
  J = ForwardDiff.derivative(uf,uprev)
  a = -dt*d
  W = 1+a*J
  #f₀ = f(t,uprev)
  a = -a
  k₁ = W\@muladd(integrator.fsalfirst + a*dT)
  dto2 = dt/2
  f₁ = f(t+dto2,@muladd(uprev+dto2*k₁))
  k₂ = W\(f₁-k₁) + k₁
  tmp = @muladd uprev + dt*k₂
  integrator.fsallast = f(t+dt,tmp)
  k₃ = W\@muladd(integrator.fsallast - c₃₂*(k₂-f₁)-2(k₁-integrator.fsalfirst)+dt*dT)
  u = uprev + dt*(k₁ + 4k₂ + k₃)/6
  if integrator.opts.adaptive
    integrator.EEst = integrator.opts.internalnorm((dt*(k₁ - 2k₂ + k₃)/6)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
  end
  integrator.k[1] = k₁
  integrator.k[2] = k₂
  @pack integrator = t,dt,u,k
end


text:
include("scopf_structjump.jl")

function main()
  if length(ARGS) < 4
    println("Usage: julia scopf_structjump_main.jl solver enable_profiling case_name lines_indexes")
    println("solver is one of PipsNlp, PipsNlpSerial or Ipopt")
    println("enable_profiling is true or false. If enable_profiling is true and solver is PipsNlp, PIPS should be compiled with NLP_TIMING option.")
    println("Cases are in 'data' directory: case9 case30 case57 case118 case300 case1354pegase case2383wp case2736sp case2737sop case2746wop case2869pegase case3012wp  case3120sp case3375wp case9241pegase")
    return
  end

  solver = ARGS[1]; shift!(ARGS)
  prof = eval(parse(ARGS[1])); shift!(ARGS)

  @timing prof tic() #t_user_model_loading

  @timing prof tic() #t_rawdata_loading

  raw = loadrawdata(ARGS[1])
  @timing prof t_rawdata_loading = toq()

  scopfmodel, scopfdata = scopf_model(raw)

  # Initial point - needed especially for pegase cases
  scopf_init_x(scopfmodel,scopfdata)

  @timing prof t_user_model_loading = toq()

  model,status = scopf_solve(scopfmodel,scopfdata, solver, prof)
  if status==:Optimal || status == :Solve_Succeeded
    scopf_structjump_outputAll(scopfmodel,scopfdata)
  end

  @timing prof begin
    @message @sprintf(" t_user_model_loading %s ", t_user_model_loading)
    @message @sprintf("t_rawdata_loading %s", t_rawdata_loading)
  end
end
#run the case 9 example
main()


text:
#=
@authors: Shunsuke Hori

Markov Perfect Equilibrium for the simple duopoly example.

See the lecture at http://quant-econ.net/jl/markov_perf.html for a
description of the model.
=#

using QuantEcon

# == Parameters == #
a0    = 10.0
a1    = 2.0
beta  = 0.96
gamma = 12.0

# == In LQ form == #

A  = eye(3)

B1 = [0.0, 1.0, 0.0]
B2 = [0.0, 0.0, 1.0]

R1 = [0.0      -a0/2.0 0.0;
      -a0/2.0  a1      a1/2.0;
      0.0      a1/2.0  0.0]

R2 = [0.0     0.0   -a0/2.0;
      0.0     0.0    a1/2.0;
      -a0/2.0 a1/2.0 a1]

Q1 = Q2 = gamma

S1 = S2 = W1 = W2 = M1 = M2 = 0.0

# == Solve using QE's nnash function == #
F1, F2, P1, P2 = nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2,
                          beta=beta)

# == Display policies == #
println("Computed policies for firm 1 and firm 2:")
println("F1 = $F1")
println("F2 = $F2")


text:
module StippleController

using InsuranceContractsController
using Genie.Renderer.Html, Genie.Renderer, SearchLight, Stipple, StippleUI
export Model, initModel, process, render

@reactive mutable struct Model <: ReactiveModel
    process::R{Bool} = false
    output::R{String} = ""
    input::R{String} = ""
    name::R{Vector{String}} = ["fliegner", "Burmeister"]
    given::R{Vector{String}} = ["michael", "irmgard"]
    age::R{Int} = 0
    warin::R{Bool} = true
end

function handlers(model)
    on(model.process) do _
        if (model.process[])
            model.output[] = model.input[] |> reverse
            model.given[1] = model.name[2] |> reverse
            model.given[2] = model.name[1] |> reverse
            println(
                "person" *
                model.name[1] *
                " " *
                model.name[2] *
                model.given[1] *
                " " *
                model.given[2],
            )
            model.process[] = false
        end
    end
    on(model.isready) do _
        push!(model)
    end
    model
end

function ui(model)
    pinputs = ""
    for i = 1:length(model.name[])
        pinputs =
            pinputs * p(
                [
                    "DynamicName Person $(i)"
                    input(
                        "",
                        @bind(Symbol("name[$(i-1)]")),
                        @on("keyup.enter", "process = true")
                    )
                ],
            )
        pinputs =
            pinputs * p(
                [
                    "Dynamic Given Person $(i)"
                    input(
                        "",
                        @bind(Symbol("given[$(i-1)]")),
                        @on("keyup.enter", "process = true")
                    )
                ],
            )
    end
    println("pinputs")
    println(pinputs)
    page(
        model,
        class = "container",
        [
            pinputs
            p([
                "Input "
                input("", @bind(:input), @on("keyup.enter", "process = true"))
            ])
            StippleUI.form(
                action = "/sub",
                method = "POST",
                [
                    textfield(
                        "What's your name *",
                        Symbol("name[0]"),
                        name = "name",
                        @iif(:warin),
                        :filled,
                        hint = "Name and surname",
                        "lazy-rules",
                        rules = "[val => val && val.length > 0 || 'Please type something']",
                    ),
                    numberfield(
                        "Your age *",
                        :age,
                        name = "age",
                        "filled",
                        :lazy__rules,
                        rules = "[val => val !== null && val !== '' || 'Please type your age',
                          val => val > 0 && val < 100 || 'Please type a real age']",
                    ),
                    btn("submit", type = "submit", color = "primary"),
                ],
            )
        ],
    )
end

function initModel(hid,vid)
    SearchLight.Configuration.load() |> SearchLight.connect
    csect=csection(hid,vid)
    global model = Stipple.init(Model, transport = Genie.WebThreads) |> handlers
    model.name = [ csect.contract_revision.description ]
end

function render(cmd::String)
    html(ui(model), context = ctx = @__MODULE__)
end

function process(cmd::String)
    age = params(:age, "99")
    name = params(:name, "Schnipperdey")
    println("PARAMS " * name * "  " * string(age))
    model.name = ["a", "b", "Schnipperdey"]
end

end

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp