text:
function secular(A::SymArrow{T}, λ::T) where T
    a, b, c = A.a, A.b, A.c
    ret = λ-c
    @inbounds for i = 1:length(a)
        ret += b[i]^2/(a[i]-λ)
    end
    ret
end

function secular_derivative(A::SymArrow{T}, λ::T) where T
    a, b, c = A.a, A.b, A.c
    ret = one(T)
    @inbounds for i = 1:length(a)
        ret += (b[i]/(a[i]-λ))^2
    end
    ret
end

function secular_second_derivative(A::SymArrow{T}, λ::T) where T
    a, b, c = A.a, A.b, A.c
    ret = zero(T)
    @inbounds for i = 1:length(a)
        ret += b[i]^2/(a[i]-λ)^3
    end
    2ret
end

function secular_over_derivative(A::SymArrow{T}, λ::T) where T
    a, b, c = A.a, A.b, A.c
    num = λ-c
    den = one(T)
    @inbounds for i = 1:length(a)
        bi2 = b[i]^2
        aiλ = a[i]-λ
        bi2daiλ = bi2/aiλ
        num += bi2daiλ
        den += bi2daiλ/aiλ
    end
    num/den
end

function pick_parameter(A::SymArrow{T}, j::Integer, λ::T) where T
    f = secular(A, λ)
    fp = secular_derivative(A, λ)
    fpp = secular_second_derivative(A, λ)
    c1 = inv(A.a[j-1]-λ)
    c2 = inv(A.a[j]-λ)
    # α, β, γ
    # γ = (fpp - 2c1*fp)/(2*(c2-c1)*c2^2)
    # β = (fp-c2^2*γ)/c1^2
    # α = f - c1*β - c2*γ
    c2γ = (fpp - 2c1*fp)/(2*(c2-c1)*c2)
    c1β = (fp-c2*c2γ)/c1
    return f - c1β - c2γ
end

function pick_zero_update(A::SymArrow{T}, j::Integer, α::T, λ::T) where T
    f = secular(A, λ)
    fp = secular_derivative(A, λ)
    a = α/((A.a[j-1]-λ)*(λ-A.a[j]))
    b = fp - (inv(A.a[j-1]-λ)+inv(A.a[j]-λ))*f
    c = -f
    return 2c/(b+sqrt(b^2-4a*c))
end

function pick_zero_update(A::SymArrow{T}, j::Integer, λ::T) where T
    f = secular(A, λ)
    fp = secular_derivative(A, λ)
    fpp = secular_second_derivative(A, λ)
    c1 = inv(A.a[j-1]-λ)
    c2 = inv(A.a[j]-λ)
    # α, β, γ
    # γ = (fpp - 2c1*fp)/(2*(c2-c1)*c2^2)
    # β = (fp-c2^2*γ)/c1^2
    # α = f - c1*β - c2*γ
    c2γ = (fpp - 2c1*fp)/(2*(c2-c1)*c2)
    c1β = (fp-c2*c2γ)/c1
    α = f - c1β - c2γ
    a = α/((A.a[j-1]-λ)*(λ-A.a[j]))
    b = fp - (inv(A.a[j-1]-λ)+inv(A.a[j]-λ))*f
    c = -f
    return 2c/(b+sqrt(b^2-4a*c))
end

function first_initial_guess(a1::T, nrmb2::T, c::T) where T
    ret = a1
    if c < a1
        amc2 = (a1-c)/2
        ret -= amc2+sqrt(amc2^2+nrmb2)
    else
        cma2 = (c-a1)/2
        ret -= nrmb2/(cma2+sqrt(cma2^2+nrmb2))
    end
    ret
end

function first_pick_zero_update(A::SymArrow{T}, λ::T) where T
    f = secular(A, λ)
    fp = secular_derivative(A, λ)
    n = length(A.a)
    a1 = A.a[1]
    α = one(T)
    @inbounds for i = 2:n
        α += (A.b[i]/(A.a[i]-λ))^2*(a1-A.a[i])/(λ-A.a[i])
    end
    a = -α/(λ-a1)
    b = fp + f/(λ-a1)
    c = -f
    return 2c/(b+sqrt(b^2-4a*c))
end

function last_initial_guess(an::T, nrmb2::T, c::T) where T
    ret = an
    if c < an
        amc2 = (an-c)/2
        ret += nrmb2/(amc2+sqrt(amc2^2+nrmb2))
    else
        cma2 = (c-an)/2
        ret += cma2+sqrt(cma2^2+nrmb2)
    end
    ret
end

function last_pick_zero_update(A::SymArrow{T}, λ::T) where T
    f = secular(A, λ)
    fp = secular_derivative(A, λ)
    n = length(A.a)
    an = A.a[n]
    α = one(T)
    @inbounds for i = 1:n-1
        α += (A.b[i]/(A.a[i]-λ))^2*(an-A.a[i])/(λ-A.a[i])
    end
    a = -α/(λ-an)
    b = fp + f/(λ-an)
    c = -f
    return 2c/(b+sqrt(b^2-4a*c))
end

# Fast Multipole Method used to accelerate the eigensolve

@inline cauchykernel{T}(x::T, y::T) = inv(x-y)
@inline coulombkernel{T}(x::T, y::T) = inv((x-y)^2)
@inline coulombprimekernel{T}(x::T, y::T) = inv((x-y)^3)
@inline logkernel{T}(x::T, y::T) = log(abs(x-y))

function secular(A::SymArrow{T}, λ::Vector{T}) where T
    a, b, c = A.a, A.b, A.c
    ret = λ-c
    for i = 1:length(λ)
        @inbounds for j = 1:length(a)
            ret[i] += b[j]^2/(a[j]-λ[i])
        end
    end
    ret
end

function secular_FMM(A::SymArrow{T}, λ::Vector{T}) where T
    a, b, c = A.a, A.b, A.c
    amin, amax = extrema(a)
    λmin, λmax = extrema(λ)
    K = KernelMatrix(cauchykernel, reverse(λ), reverse(a), λmax, λmin, amax, amin)
    #ret = λ - c + K*(b.^2)
    ret = λ - c - reverse!(K*(b.^2))
    ret
end

function secular_derivative(A::SymArrow{T}, λ::Vector{T}) where T
    a, b, c = A.a, A.b, A.c
    ret = ones(T, length(λ))
    for i = 1:length(λ)
        @inbounds for j = 1:length(a)
            ret[i] += (b[j]/(a[j]-λ[i]))^2
        end
    end
    ret
end

function secular_derivative_FMM(A::SymArrow{T}, λ::Vector{T}) where T
    a, b, c = A.a, A.b, A.c
    amin, amax = extrema(a)
    λmin, λmax = extrema(λ)
    K = KernelMatrix(coulombkernel, reverse(λ), reverse(a), λmax, λmin, amax, amin)
    ret = ones(T, length(λ)) + K*(b.^2)
    reverse!(ret)
end

function secular_second_derivative(A::SymArrow{T}, λ::Vector{T}) where T
    a, b, c = A.a, A.b, A.c
    ret = zeros(T, length(λ))
    for i = 1:length(λ)
        @inbounds for j = 1:length(a)
            ret[i] += 2b[j]^2/(a[j]-λ[i])^3
        end
    end
    ret
end

function secular_second_derivative_FMM(A::SymArrow{T}, λ::Vector{T}) where T
    a, b, c = A.a, A.b, A.c
    amin, amax = extrema(a)
    λmin, λmax = extrema(λ)
    K = KernelMatrix(coulombprimekernel, reverse(λ), reverse(a), λmax, λmin, amax, amin)
    ret = K*(-2b.^2)
    reverse!(ret)
end

# For all 2 ≤ j ≤ n-1
function pick_zero_update(A::SymArrow{T}, λ::Vector{T}) where T
    n = size(A, 1)
    f = secular(A, λ)
    fp = secular_derivative(A, λ)
    fpp = secular_second_derivative(A, λ)
    c1 = T[inv(A.a[j-1]-λ[j-1]) for j in 2:n-1]
    c2 = T[inv(A.a[j]-λ[j-1]) for j in 2:n-1]
    c2γ = (fpp - 2c1.*fp)./(2*(c2-c1).*c2)
    c1β = (fp-c2.*c2γ)./c1
    α = f - c1β - c2γ
    a = T[α[j-1]/((A.a[j-1]-λ[j-1])*(λ[j-1]-A.a[j])) for j in 2:n-1]
    b = T[fp[j-1] - (inv(A.a[j-1]-λ[j-1])+inv(A.a[j]-λ[j-1]))*f[j-1] for j in 2:n-1]
    c = -f
    return 2c./(b+sqrt.(b.^2-4a.*c))
end

function pick_zero_update_FMM(A::SymArrow{T}, λ::Vector{T}) where T
    n = size(A, 1)
    f = secular_FMM(A, λ)
    fp = secular_derivative_FMM(A, λ)
    fpp = secular_second_derivative_FMM(A, λ)
    c1 = T[inv(A.a[j-1]-λ[j-1]) for j in 2:n-1]
    c2 = T[inv(A.a[j]-λ[j-1]) for j in 2:n-1]
    c2γ = (fpp - 2c1.*fp)./(2*(c2-c1).*c2)
    c1β = (fp-c2.*c2γ)./c1
    α = f - c1β - c2γ
    a = T[α[j-1]/((A.a[j-1]-λ[j-1])*(λ[j-1]-A.a[j])) for j in 2:n-1]
    b = T[fp[j-1] - (inv(A.a[j-1]-λ[j-1])+inv(A.a[j]-λ[j-1]))*f[j-1] for j in 2:n-1]
    c = -f
    return 2c./(b+sqrt.(b.^2-4a.*c))
end


text:
ENV["JULIA_DEBUG"] = "HierarchicalTemporalMemory"
plot_enabled= false

using HierarchicalTemporalMemory, BenchmarkTools, CSV, Printf, Lazy, DataFrames, Test
using Plots; gr()
import Random.seed!
seed!(0)

include("utils/utils.jl")

display_evaluation(t,sp,sp_activity,spDims)= begin
  @info("t=$t")
  sparsity= count(sp_activity)/prod(spDims)*100

  # Evaluation metrics
  #@info "Sparsity: $(sparsity)"
  #histogram(sp.b.b)|> display
  #histogram(sp.synapses.Dₚ[sp.synapses.Dₚ.>0]|>Vector, yaxis=(:log10))|> display
  #heatmap(@> sp_activity reshape(64,32))|> display
  #heatmap(@> sp.b.a_Tmean reshape(64,32))|> display
  #heatmap(sp.synapses.Dₚ)|> display
end
process_data!(encHistory,spHistory,encANDspHistory,tN,data,encParams,sp)=
  for t in 1:tN
    z,a= _process_sp(t,tN,data,encParams,sp,display_evaluation)
    encHistory[:,t]= z; spHistory[:,t]= a
    # when were the most similar SDRs to z,a in history? These indices correspond to times
    _,similarEncIdx= top_similar_sdr(encHistory[:,1:t],10)
    _,similarSpIdx= top_similar_sdr(spHistory[:,1:t],10)
    encOnlyIdx= setdiff(similarEncIdx, similarSpIdx)
    spOnlyIdx= setdiff(similarSpIdx, similarEncIdx)
    encANDspHistory[t]= (encANDsp= intersect(similarEncIdx,similarSpIdx), Nenc= length(similarEncIdx))
    t%60==0 && plot_enabled &&
        plot_ts_similarEncSp(t,data.power_hourly_kw,
                             encOnlyIdx,spOnlyIdx,encANDspHistory)
  end

# Define Spatial Pooler
inputDims= ((14,6,4).*25,)
spDims= (2048,)
#inputDims= (8,8)
#spDims= (12,12)
@info "creating Spatial Pooler"
sp= SpatialPooler(SPParams(
      szᵢₙ= (inputDims.|> sum), szₛₚ=spDims,
      γ=600,
      s=0.02,
      prob_synapse=0.15,
      θ_stimulus_activate=4,
      p⁺_01= 0.07,
      p⁻_01= 0.12,
      β=5,
      Tboost=400,
      enable_local_inhibit=false,
      enable_boosting=true))
# Define input data
data,tN= read_gympower("test_data/gym_power_benchmark.csv")
encParams= initenc_powerDay(data.power_hourly_kw, data.hour, data.is_weekend,
                 encoder_size=inputDims[1], w=(23,27,27))
# Histories
encHistory= falses(map(sum,inputDims)|>prod,tN)
spHistory= falses(spDims|>prod,tN)
encANDspHistory= Vector{NamedTuple{(:encANDsp,:Nenc),Tuple{Vector{Int},Int}}}(undef,tN)

@info("processing data")
process_data!(encHistory,spHistory,encANDspHistory, tN,data,encParams,sp)
total_overlap= [1; map(x->length(x.encANDsp), encANDspHistory[2:end]) ./
                   map(x->x.Nenc, encANDspHistory[2:end])]
#plot(total_overlap)|>display
early_totalOverlap= mean(total_overlap[50:100])
late_totalOverlap= mean(total_overlap[500:end])
@info @sprintf("Mean SP performance: [%.2f,%.2f]\n", early_totalOverlap, late_totalOverlap)

# If this isn't true, something's quite wrong with the model
@test early_totalOverlap >= 0.65 && late_totalOverlap >= 0.78 && late_totalOverlap > early_totalOverlap

#sparsity= (i-> count(spHistory[:,i])/prod(spDims)).(1:tN)
#mean(sparsity.*100)|> display
#using Statistics
#std(sparsity.*100)|> display

ext: jl


text:
export WindingGeometry, Windings, windings, winding_geometry

"""
    WindingGeometry

2D representation of one layer of a winding.

**Fields**
- `width`   -- width of trace (m)
- `length`  -- length of trace (m)
- `turns`   -- number of turns

`WindingGeometry` objects are typicaly created with `winding_geometry`
"""
struct WindingGeometry
  width :: Float64
  length :: Float64
  turns :: Int
end
"""
    winding_geometry(pcb, core, turns)

Create a `WindingGeometry` object.
"""
function winding_geometry(pcb :: PCB_Specification,
                          core :: CoreGeometry,
                          turns :: Int)
  trace_width = (core.winding_aperture-2*pcb.trace_edge_gap-
    (turns-1)*pcb.trace_trace_gap)/turns
  trace_length = 0.0
  for i in 0:turns-1
    r = core.half_center_width+pcb.trace_edge_gap+
        0.5*trace_width+i*(trace_width+pcb.trace_trace_gap)
    trace_length += 2π*r
  end
  trace_length += 2*core.center_length
  if trace_width<=0.0
    warn("$(turns) on $(name(core)) will result in negative trace width.  Setting to NaN.")
    trace_width = NaN
  end
  WindingGeometry(trace_width, trace_length, turns)
end

"""
    Windings

**Fields**
- `pcb`                     -- `PCB_Specification`
- `core`                    -- `CoreGeometry`
- `primarywindinggeometry`  -- `WindingGeometry` for all layers of primary
- `secondarywindinggeometry`-- `WindingGeometry` for all layers of secondary
- `isprimary`               -- tuple of `Bool` for each conductor layer
- `sides`                   -- for internal use
- `isprimaryseries`         -- true if primary is connected in series
- `issecondaryseries`       -- true if secondary is connected in series

Create `Windings` objects with `windings`.
"""
struct Windings
  pcb :: PCB_Specification
  core :: CoreGeometry
  primarywindinggeometry :: WindingGeometry
  secondarywindinggeometry :: WindingGeometry
  isprimary :: Tuple
  sides :: Tuple
  isprimaryseries :: Bool
  issecondaryseries :: Bool
end
"""
    windings(pcb, core,
            primaryturnsperlayer, secondaryturnsperlayer,
            isprimary,
            isprimaryeries, issecondaryseries)

Create a `Windings` object.
"""
function windings(pcb::PCB_Specification, core::CoreGeometry,
                  primaryturnsperlayer::Int, secondaryturnsperlayer::Int,
                  isprimary,
                  isprimaryeries::Bool, issecondaryseries::Bool)
  if length(isprimary)!=div(length(pcb.stackup.material),2)+1
    throw(ArgumentError("length of isprimary must be the same as number of conductor layers"))
  end
  x = 0
  for ip in isprimary
    x += ip ? 1 : -1
    if abs(x)>1
      throw(ArgumentError("cannot analyze due to layer connections to primay and secondary"))
    end
  end
  Windings(pcb,core,
           winding_geometry(pcb,core,primaryturnsperlayer),
           winding_geometry(pcb,core,secondaryturnsperlayer),
           isprimary,
           ntuple(i->sides(isprimary,i),length(isprimary)),
           isprimaryeries, issecondaryseries
          )
end
"""
    sides(isprimary,index)

For a given layer (index), returns the number (1 or 2) of sides which are part
of the same winding.

This function is for internal use calculating effective_resistance and
 leakage_inductance.
"""
function sides(isprimary,i)::Float64
  i==1 && return 1.0
  i==length(isprimary) && return 1.0
  s = 0.0
  isprimary[i] != isprimary[i-1] && (s+=1.0)
  isprimary[i] != isprimary[i+1] && (s+=1.0)
  return s
end

effective_volume(w::Windings) = effective_volume(core(w))
effective_area(w::Windings) = effective_area(core(w))
effective_length(w::Windings) = effective_length(core(w))
winding_aperture_height(w::Windings) = winding_aperture_height(core(w))
winding_aperture(w::Windings) = winding_aperture(core(w))
core_name(w::Windings) = name(core(w))


text:
module NNConv1d

include("impl.jl")
include("flux.jl")

end


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("CovFuncs.jl")
include("Utils.jl")
using LinearAlgebra
using Distributions
using Distances
using Random
using Match
using JLD
Random.seed!(123)

N = 10000
d = 2
l = 0.5
σList = [1e3, 1e2, 1e1, 1e0, 1e-1, 1e-2, 1e-3, 1e-4]
ρ = 3.0
λ = 1.5

# Creating uniform random data
x = rand(d, N)

# reordering points and forming skeletons
colptr, rowval, P, revP, distances = sortSparseRev(x, ρ, 1)
xOrd = x[:, P]
skeletons = construct_skeletons(colptr, rowval, distances, λ)

KLTrueNaive = zeros(3, length(σList))
KLNaiveTrue = zeros(3, length(σList))
KLTrueApp = zeros(3, length(σList))
KLAppTrue = zeros(3, length(σList))
KLTrueLarge = zeros(3, length(σList))
KLLargeTrue = zeros(3, length(σList))
KLTrueExact = zeros(3, length(σList))
KLExactTrue = zeros(3, length(σList))

σBase = exp.(randn(N))
σBase .= max.(σBase, eps(Float64))

covfunc12!(r) = matern12!(r,l)
covfunc32!(r) = matern32!(r,l)
covfunc52!(r) = matern52!(r,l)
for (iOrder, maternOrder) in enumerate(["12", "32", "52"])
  # Make sure to avoid NANs
  @match maternOrder begin
    "12" => (covfunc! = covfunc12!)
    "32" => (covfunc! = covfunc32!)
    "52" => (covfunc! = covfunc52!)
  end

  for (iScale, σScale) in enumerate(σList)
    σ = σScale * σBase
    # Creating reference solution
    Θ = covfunc!(pairwise(Euclidean(), xOrd, dims=2)) + Diagonal(σ)

    # Creating naive vecchia approximation
    LNaive = assembleL(skeletons, xOrd, covfunc!, σ)
    ΘNaive = inv(Matrix(LNaive)') * inv(Matrix(LNaive))

    # Creating "smart" vecchia approximation
    L = assembleL(skeletons, xOrd, covfunc!)
    U = squareSparse(L) + spdiagm(0 => 1 ./ σ)
    icholU_high_level!(U)
    noisecov = NoiseCov(L, U, σ)
    ΘExact = exactMatrix(noisecov)
    ΘApp = appMatrix(noisecov)

    # Creating smart vecchia approximation with larger stencil for IC0
    ULarge = triu(L * L' + spdiagm(0 => 1 ./ σ))
    icholU_high_level!(ULarge)
    noisecovLarge = NoiseCov(L, ULarge, σ)
    ΘLarge = appMatrix(noisecovLarge)

    KLTrueNaive[iOrder, iScale] = KL(Θ, ΘNaive)
    KLNaiveTrue[iOrder, iScale] = KL(ΘNaive, Θ)

    KLTrueApp[iOrder, iScale] = KL(Θ, ΘApp)
    KLAppTrue[iOrder, iScale] = KL(ΘApp, Θ)

    KLTrueLarge[iOrder, iScale] = KL(Θ, ΘLarge)
    KLLargeTrue[iOrder, iScale] = KL(ΘLarge, Θ)

    KLTrueExact[iOrder, iScale] = KL(Θ, ΘExact)
    KLExactTrue[iOrder, iScale] = KL(ΘExact, Θ)
  end
end

save("./out/jld/KLVarySigma.jld",
     "σList", σList,
     "KLTrueNaive", KLTrueNaive,
     "KLNaiveTrue", KLNaiveTrue,
     "KLTrueApp", KLTrueApp,
     "KLAppTrue", KLAppTrue,
     "KLTrueLarge", KLTrueLarge,
     "KLLargeTrue", KLLargeTrue,
     "KLTrueExact", KLTrueExact,
     "KLExactTrue", KLExactTrue)


text:
using POMDPs, StatsBase, POMDPToolbox, SARSOP, QMDP, PyCall
using AutomotiveDrivingModels, Parameters, ArgParse, AutoViz, Reel

s = ArgParseSettings()
@add_arg_table s begin
    "--weights"
        help = "the path to the NN weights to load"
        arg_type = String
        default = ""
end

parsed_args  = parse_args(ARGS, s)
weight_path = parsed_args["weights"]

include("../admenv/crosswalk/occluded_crosswalk_env.jl") # env config
include("eval-env/constant_pedestrian.jl") # pedestrian model
include("eval-env/pedestrian_flow.jl")
include("eval-env/simple_sensor.jl")
include("eval-env/ego_control.jl")
include("simulation.jl")
include("eval_config.jl") # eval config
include("baseline_policy.jl")
include("render_helpers.jl")

config = EvalConfig(n_episodes=100, time_out = 200, n_ped=1) # use default
params = EnvParams(ped_rate = 0.9)
env = CrosswalkEnv(params)
cam = FitToContentCamera(0.0)

include("../admenv/crosswalk/pomdp_types.jl")
include("../admenv/crosswalk/spaces.jl")
include("../admenv/crosswalk/distributions.jl")
include("../admenv/crosswalk/belief.jl")
include("../adm-workspace/adm_helpers.jl")
include("../admenv/crosswalk/decomposition.jl")
include("eval-env/qmdp_eval_model.jl")
include("eval-env/nn_eval_model.jl")
include("eval-env/continuous_observations.jl")

@pyimport tensorflow as tf
nn_wrapper = pyimport("a3c.nn_wrapper")

sess = tf.Session();
rng = MersenneTwister(3);
pomdp = OCPOMDP();
pomdp.pos_res = 0.5
pomdp.vel_res = 0.5
pomdp.p_birth = 0.2
up = OCUpdater(pomdp);
nn = nn_wrapper[:NNPolicy](sess, weight_path, n_actions(pomdp))
policy = NNPolicy(pomdp, nn);
sensor = POMDPSensor(config.n_ped, pomdp, SimpleSensor(0.1, 0.1))
d0 = initialstate_distribution(pomdp)
crash = 0
success = 0
time_out = 0
avg_step = 0
crash_flag = true
for i=1:config.n_episodes
    s0 = rand(config.rng, d0)
    ego_initialstate = s0.ego
    ego = Vehicle(ego_initialstate, VehicleDef(), 1)
    #TODO function initialize belief
    b0 = Dict{Int64, NNBelief}()
    for j=1:config.n_ped
        b0[j+1] = Float64[]
    end
    models = Dict{Int, DriverModel}()
    up = NNUpdater(pomdp);
    reset_policy!(policy)
    models[1] = CrosswalkDriver(5, 0, env, 0., b0, Dict{Int64, OCState}(), ego, sensor, policy, up);
    # todo function init_scene
    scene = Scene()
    push!(scene, ego)
    nticks = 200
    rec = SceneRecord(nticks+1, config.sim_dt)

    # execute the simulation
    simulate!(rec, scene, env, models, nticks, config.rng, config.n_ped, config.callbacks)
    avg_step += nframes(rec)

    # analyse output
    if is_crash(rec[0]) && crash_flag
        crash += 1
        duration, fps, render_rec = animate_record(rec)
        speed_factor = 1
        film = roll(render_rec, fps = fps*speed_factor, duration = duration/speed_factor)
        write("crash.mp4", film)
        crash_flag = false
        # println("Episode $i CRASH")
        # break
    elseif is_terminal(rec[0], env)
        success += 1
        # println("Episode $i SUCCESS")
    elseif nframes(rec) >= config.time_out
        time_out += 1
        # println("Episode $i TIME OUT")
    else
        # println("Episode $i nothing happened, frames: $(nframes(rec))")
    end
end
println("End of the simulation: Crashes $crash , Successes $success, Time out $time_out, Avg steps $(avg_step/config.n_episodes)")

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
#Defines the Guidance system type, including its discrete time behaviour
#include("DiscreteTimeSystem.jl")

#Define ControlSystem as a subtype of the DiscreteTimeSystem
abstract type AbstractGuidanceSystem <: DiscreteTimeSystem end

########### Template function ###########################
mutable struct GuidanceSystem{T<:Number} <: AbstractGuidanceSystem
  Active::Bool
  Ts::T #sample time
  GuidanceRelation::Function
  Tprev::T #storage for previous time it was updated, hence must be a mutable type
end
Guidance() = GuidanceSystem(false, 0.0, x -> nothing, 0.0)
GuidanceSystem(Ts::Number) = GuidanceSystem(false, Float64(Ts), x -> nothing, 0.0)
GuidanceSystem(Ts::Number, csfunc::Function) = GuidanceSystem(true, Float64(Ts), csfunc, 0.0)
GuidanceSystem(Active::Bool, Ts::Number, csfunc::Function) = GuidanceSystem(Active, Float64(Ts), csfunc, 0.0)
#Set the callback affect condition of the ControlSystem
#placeholder!
function affect!(int, dts::Type{<:AbstractGuidanceSystem})
    for cu in user_cache(int)
        cu.u = [dts.GuidanceRelation(int.u)...]
    end
  dts.Tprev = int.t
end


text:
#' ---
#' title: Finance
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' Optimization models play an increasingly important role in financial decisions.
#' Many computational finance problems can be solved efficiently using modern optimization techniques.
#' In this tutorial we will discuss 3 such examples.
#' The first two are taken from the book Optimization Methods in Finance[[1]](#c1),
#' while the third is the famous Markowitz Portfolio Optimization problem with data from
#' lecture notes from a course taught at Georgia Tech[[2]](#c2).

using JuMP
using GLPK

#' ## Short Term Financing

#' Corporations routinely face the problem of financing short term cash commitments such as the following:

#' **Month**        |Jan    |Feb    |Mar    |Apr    |May    |Jun
#' :-----:          |:-----:|:-----:|:-----:|:-----:|:-----:|:-----:
#' **Net Cash Flow**|-150   |-100   |200    |-200   |50     |300

#' Net cash flow requirements are given in thousands of dollars. The company has the following sources of funds:

#' - A line of credit of up to $100K at an interest rate of 1% per month,
#' - In any one of the first three months, it can issue 90-day commercial paper bearing a total interest of 2% for the 3-month period,
#' - Excess funds can be invested at an interest rate of 0.3% per month.

#' Our task is to find out the most economical way to use these 3 sources such that
#' we end up with the most amount of money at the end of June.
#' We model this problem in the following manner:

#' We will use the following decision variables:
#' - the amount $u_{i}$ drawn from the line of credit in month $i$
#' - the amount $v_{i}$ of commercial paper issued in month $i$
#' - the excess funds $w_{i}$ in month $i$

#' Here we have three types of constraints:
#' 1. for every month, cash inflow = cash outflow for each month
#' 2. upper bounds on $u_{i}$
#' 3. nonnegativity of the decision variables $u_{i}$, $v_{i}$ and $w_{i}$.

#' Our objective will be to simply maximimse the company's wealth in June, which say we represent with the variable $m$.

financing = Model(GLPK.Optimizer)

@variables(financing, begin
    0 <= u[1:5] <= 100
    0 <= v[1:3]
    0 <= w[1:5]
    m
end)

@objective(financing, Max, m) # Money at the end of June

@constraints(financing, begin
    u[1] + v[1] - w[1] == 150 # January
    u[2] + v[2] - w[2] - 1.01u[1] + 1.003w[1] == 100 # February
    u[3] + v[3] - w[3] - 1.01u[2] + 1.003w[2] == -200 # March
    u[4] - w[4] - 1.02v[1] - 1.01u[3] + 1.003w[3] == 200 # April
    u[5] - w[5] - 1.02v[2] - 1.01u[4] + 1.003w[4] == -50 # May
    -m - 1.02v[3] - 1.01u[5] + 1.003w[5] == -300 # June
end)

optimize!(financing)
@show objective_value(financing);

#' ## Combinatorial Auctions
#' In many auctions, the value that a bidder has for a set of items may not be
#' the sum of the values that he has for individual items.
#' It may be more or it may be less.
#' Examples are equity trading, electricity markets, pollution right auctions and auctions for airport landing slots.
#' To take this into account, combinatorial auctions allow the bidders to submit bids on combinations of items.

#' Let $M=\{1,2, \ldots, m\}$ be the set of items that the auctioneer has to sell.
#' A bid is a pair $B_{j}=\left(S_{j}, p_{j}\right)$ where $S_{j} \subseteq M$ is a nonempty set of items and
#' $p_{j}$ is the price offer for this set.
#' Suppose that the auctioneer has received $n$ bids $B_{1}, B_{2}, \ldots, B_{n} .$
#' The goal of this problem is to help an auctioneer determine the winners in order to maximize his revenue.

#' We model this problem by taking a decision variable $y_{j}$ for every bid.
#' We add a constraint that each item $i$ is sold at most once. This gives us the following model:

#' $$
#' \begin{align*}
#' \max && \sum_{i=1}^{n} p_{j} y_{j} \\
#' \text { s.t. }  && \sum_{j : i \in S_{j}} y_{j} \leq 1 && \forall i=\{1,2 \ldots m\} \\
#' && y_{j} \in\{0,1\} && \forall j \in\{1,2 \ldots n\}
#' \end{align*}
#' $$

bid_values = [6 3 12 12 8 16]
bid_items = [[1], [2], [3 4], [1 3], [2 4], [1 3 4]]

auction = Model(GLPK.Optimizer)
@variable(auction, y[1:6], Bin)
@objective(auction, Max, sum(y' .* bid_values))
for i in 1:6
    @constraint(auction, sum(y[j] for j in 1:6 if i in bid_items[j]) <= 1)
end

optimize!(auction)

@show objective_value(auction);
@show value.(y);

#' ## Portfolio Optimization
#' Suppose we are considering investing 1000 dollars in three non-dividend paying stocks,
#' IBM (IBM), Walmart (WMT), and Southern Electric (SEHI), for a onemonth period.
#' This means we will use the money to buy shares of the three stocks at the current market prices,
#' hold these for one month, and sell the shares off at the prevailing market prices at the end of the month.
#' As a rational investor, we hope to make some profit out of this endeavor, i.e.,
#' the return on our investment should be positive.
#' Suppose we bought a stock at $p$ dollars per share in the beginning of the month,
#' and sold it off at $s$ dollars per share at the end of the month.
#' Then the one-month return on a share of the stock is $ \frac{s-p}{p} $.

#' Since the stock prices are quite uncertain, so is the end-of-month return on our investment.
#' Our goal is to invest in such a way that the expected end-of-month return is at least $50.00 or 5%.
#' Furthermore, we want to make sure that the “risk” of not achieving our desired return is minimum.

#' Note that we are solving the problem under the following assumptions:
#' 1. We can trade any continuum of shares.
#' 2. No short-selling is allowed.
#' 3. There are no transaction costs.

#' We model this problem by taking decision variables $x_{i}, i=1,2,3,$ denoting the dollars invested in each of the 3 stocks.
#' Let us denote by $\tilde{r}_{i}$ the random variable corresponding to the monthly return
#' (increase in the stock price) per dollar for stock $i .$

#' Then, the return (or profit) on $x_{i}$ dollars invested in stock $i$ is $\tilde{r}_{i} x_{i},$ and
#' the total (random) return on our investment is $\sum_{i=1}^{3} \tilde{r}_{i} x_{i} .$
#' The expected return on our investment is then
#' $\mathbb{E}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right]=\sum_{i=1}^{3} \overline{r}_{i} x_{i},$
#' where $\overline{r}_{i}$ is the expected value of the $\tilde{r}_{i} .$

#' Now we need to quantify the notion of “risk” in our investment.
#' Markowitz, in his Nobel prize winning work, showed that a rational investor’s notion of minimizing risk
#' can be closely approximated by minimizing the variance of the return of the investment portfolio.
#' This variance is given by:

#' $$
#' \operatorname{Var}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right] = \sum_{i=1}^{3} \sum_{j=1}^{3} x_{i} x_{j} \sigma_{i j}
#' $$

#' where $\sigma_{i j}$ is the covariance of the return of stock $i$ with stock $j$.

#' Note that the right hand side of the equation is the most reduced form of the expression and
#' we have not shown the intermediate steps involved in getting to this form.
#' We can also write this equation as:

#' $$
#' \operatorname{Var}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right] =x^{T} Q x
#' $$

#' Where $Q$ is the covariance matrix for the random vector $\tilde{r}$.

#' Finally, we can write the model as:

#' $$
#' \begin{align*}
#' \min x^{T} Q x \\
#' \text { s.t. } \sum_{i=1}^{3} x_{i} \leq 1000.00 \\
#' \overline{r}^{T} x \geq 50.00 \\
#' x \geq 0
#' \end{align*}
#' $$

#' After that long discussion, lets now use JuMP to solve the portfolio optimization problem for the data given below.

#' | Month        |  IBM     |  WMT    |  SEHI  |
#' |--------------|----------|---------|--------|
#' | November-00  |  93.043  |  51.826 |  1.063 |
#' | December-00  |  84.585  |  52.823 |  0.938 |
#' | January-01   |  111.453 |  56.477 |  1.000 |
#' | February-01  |  99.525  |  49.805 |  0.938 |
#' | March-01     |  95.819  |  50.287 |  1.438 |
#' | April-01     |  114.708 |  51.521 |  1.700 |
#' | May-01       |  111.515 |  51.531 |  2.540 |
#' | June-01      |  113.211 |  48.664 |  2.390 |
#' | July-01      |  104.942 |  55.744 |  3.120 |
#' | August-01    |  99.827  |  47.916 |  2.980 |
#' | September-01 |  91.607  |  49.438 |  1.900 |
#' | October-01   |  107.937 |  51.336 |  1.750 |
#' | November-01  |  115.590 |  55.081 |  1.800 |

using Statistics # Useful for calculations
using Ipopt      # Ipopt since our objective is quadratic

stock_data = [
93.043 51.826 1.063;
84.585 52.823 0.938;
111.453 56.477 1.000;
99.525 49.805 0.938;
95.819 50.287 1.438;
114.708 51.521 1.700;
111.515 51.531 2.540;
113.211 48.664 2.390;
104.942 55.744 3.120;
99.827 47.916 2.980;
91.607 49.438 1.900;
107.937 51.336 1.750;
115.590 55.081 1.800;
]

# Calculating stock returns

stock_returns = Array{Float64}(undef, 12, 3)

for i in 1:12
    stock_returns[i, :] = (stock_data[i + 1, :] .- stock_data[i, :]) ./ stock_data[i, :]
end

# Calculating the expected value of monthly return

r = [Statistics.mean(stock_returns[:,1]) Statistics.mean(stock_returns[:,2]) Statistics.mean(stock_returns[:,3])]'

# Calculating the covariance matrix Q

Q = Statistics.cov(stock_returns)

#+

# JuMP Model

portfolio = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(portfolio, x[1:3] >= 0)
@objective(portfolio, Min, x' * Q * x)
@constraint(portfolio, sum(x) <= 1000)
@constraint(portfolio, sum(r .* x) >= 50)

optimize!(portfolio)

@show objective_value(portfolio);
@show value.(x);

#' ### References
#'
#' 1. Cornuéjols, G., Peña, J., & Tütüncü, R. (2018). Optimization Methods in Finance. Cambridge: Cambridge University Press. doi:10.1017/9781107297340
#'
#' 2. Ahmed, S. ISyE6669: Deterministic Optimization. Available at: https://www2.isye.gatech.edu/~sahmed/isye6669/


text:
# --- PRIVATE METHODS ---------------------------------------------------------- #
function is_kegg_organism_code_ok(kegg_organism_code::String)
end

function is_cobra_dictionary_ok(cobra_dictionary::Dict{String,Any})
end

function lookup_gene_ec_mapping_record(kegg_gene_id::String)::Union{Set{String},Nothing}

    # initialize -
    record_set = Set{String}()

    # call the KEGG API to get the ec number -
    list_of_ec_numbers = get_ec_number_for_gene(kegg_gene_id)
    if list_of_ec_numbers != nothing

        # create a record, add to the mapping buffer -
        for ec_number in list_of_ec_numbers
            push!(record_set,ec_number)
        end

        # return -
        return record_set
    end

    return nothing
end
# ------------------------------------------------------------------------------ #

# --- PUBLIC METHODS ----------------------------------------------------------- #
function load_cobra_model_file(path_to_cobra_mat_file::String, model_name::String)::Dict{String,Any}

    # check file ok?
    is_file_path_ok(path_to_cobra_mat_file)

    # open the cobra file -
    file = matopen(path_to_cobra_mat_file)
    cobra_dictionary = read(file,model_name)
    close(file)

    # return -
    return cobra_dictionary
end

function extract_files_from_cobra_model(cobra_dictionary::Dict{String,Any}, kegg_organism_code::String, path_to_model_directory::String)

    # check - does this path exist?
    is_dir_path_ok(path_to_model_directory)

    # check the cobra dictionary - is it legit?
    is_cobra_dictionary_ok(cobra_dictionary)

    # ok, looks ok if we get here - write out a bunch of files that we'll need

    # reaction mapping -
    path_to_reaction_mapping_file = joinpath(path_to_model_directory,"ReactionNameMap.dat")
    export_reaction_mapping_file(cobra_dictionary,path_to_reaction_mapping_file)

    # gene order -
    path_to_gene_order_file = joinpath(path_to_model_directory,"GeneOrder.dat")
    export_gene_order_file(cobra_dictionary, kegg_organism_code, path_to_gene_order_file)

    # reaction tags to gene mapping -
    path_reaction_gene_mapping_file = joinpath(path_to_model_directory,"ReactionGeneMap.dat")
    export_reaction_tag_to_gene_mapping_file(cobra_dictionary, kegg_organism_code, path_reaction_gene_mapping_file)

    # reaction tag to ec mapping -
    path_reaction_to_ec_number_mapping = joinpath(path_to_model_directory, "ReactionECNumberMap.dat")
    export_reaction_tag_to_ec_mapping_file(cobra_dictionary, kegg_organism_code, path_reaction_to_ec_number_mapping)
end

function write_cobra_model_file(path_to_cobra_mat_file::String,cobra_dictionary::Dict{String,Any})

    # check the dir -
    is_dir_path_ok(path_to_cobra_mat_file)

    # get the model name -
    model_name = splitdir(splitext(path_to_cobra_mat_file)[1])[2]

    # ok, the dir is legit, so let's write out the cobra file -
    file = matopen(path_to_cobra_mat_file, "w")
    write(file, "$(model_name)", cobra_dictionary)
    close(file)
end

function generate_cobra_dictionary_from_vff(path_to_vff_file::String, path_to_mapping_file_dir::String)::Dict{String,Any}

    # check the path -
    is_file_path_ok(path_to_vff_file)

    # does the mapping file dir exist?
    is_dir_path_ok(path_to_mapping_file_dir)

    # initialize -
    cobra_dictionary = Dict{String,Any}()

    # parse the sections of the vff -
    problem_dictionary = load_vff_model_file(path_to_vff_file)

    # get some stuff from the problem dictionary -
    list_of_reactions = problem_dictionary["list_of_reactions"]
    reaction_order_array = problem_dictionary["reaction_order_array"]
    gene_order_array = problem_dictionary["gene_order_array"]
    list_of_genes = problem_dictionary["list_of_genes"]

    # let's fill in the fields for the cobra dictionary -
    # rxns: reaction tags -
    cobra_dictionary["rxns"] = reaction_order_array

    # rxnNames: build reaction name array -
    reaction_name_array = Array{String,1}()
    path_reaction_name_mapping_file = joinpath(path_to_mapping_file_dir,"ReactionNameMap.dat")
    is_file_path_ok(path_reaction_name_mapping_file)
    reaction_tag_name_map = build_mapping_dictionary(path_reaction_name_mapping_file)
    for reaction_tag in reaction_order_array

        # check, do we have this tag in the map?
        if (haskey(reaction_tag_name_map,reaction_tag) == true)
            reaction_name_mapping_wrapper = reaction_tag_name_map[reaction_tag]
            value_set = reaction_name_mapping_wrapper.value
            push!(reaction_name_array, pop!(value_set))
        else
            push!(reaction_name_array,reaction_tag)
        end
    end
    cobra_dictionary["rxnNames"] = reaction_name_array

    # genes: list of genes -
    cobra_dictionary["genes"] = gene_order_array

    # rules: rules -
    rule_text_array = Array{String,1}()
    list_of_rules = problem_dictionary["list_of_rules"]
    for reaction_tag in reaction_order_array

        # do we have this reaction_tag in our rule list?
        if (haskey(list_of_rules,reaction_tag) == true)
            rule_text = list_of_rules[reaction_tag].rule_text
            if (rule_text != "[]")
                push!(rule_text_array,rule_text)
            else
                push!(rule_text_array,"")
            end
        else
            push!(rule_text_array,"")
        end
    end
    cobra_dictionary["rules"] = rule_text_array

    # rxnECNumbers: Add the ec numbers of the cobra dictionary -
    ec_number_array = Array{String,1}()
    for reaction_tag in reaction_order_array

        # check - do we have this reaction tag?
        if (haskey(list_of_reactions,reaction_tag) == true)

            # get the wrapper -
            reaction_wrapper = list_of_reactions[reaction_tag]

            # get the ec_number field -
            ec_number_record = reaction_wrapper.ec_number

            # cache -
            push!(ec_number_array, ec_number_record)
        else
            # for some reason, we don't have this reaction tag?
            # if this happens, we have some crazy sh*tnitz happening, better throw an exception
            throw(ErrorException("reaction tag not in the list of reactions."))
        end
    end
    cobra_dictionary["rxnECNumbers"] = ec_number_array

    # mets: add metabolite symbols -
    # get list of reaction wrappers in the proper order -
    reaction_wrapper_array = Array{CBMetabolicReaction,1}()
    for reaction_tag in reaction_order_array

        # do we have this reaction?
        reaction_wrapper = list_of_reactions[reaction_tag]
        push!(reaction_wrapper_array, reaction_wrapper)
    end

    # Get the metabolites from the reaction list -
    metabolite_wrapper_array = build_metabolite_symbol_array(reaction_wrapper_array)

    # lb: lower flux bound (from the vff) -
    lb_array = Array{Float64,1}()
    rev_array = Float64[]
    for reaction_wrapper in reaction_wrapper_array

        # get lower bound -
        reverse_bound_value = reaction_wrapper.reverse
        if (reverse_bound_value == "-inf")
            push!(lb_array,-1000.0)
        else
            push!(lb_array,parse(Float64,reverse_bound_value))
        end

        if (parse(Float64,reverse_bound_value)<0.0)
            push!(rev_array,1.0)
        else
            push!(rev_array,0.0)
        end
    end
    cobra_dictionary["lb"] = lb_array
    cobra_dictionary["rev"] = rev_array

    # ub: upper flux bound (from the vff) -
    ub_array = Array{Float64,1}()
    for reaction_wrapper in reaction_wrapper_array

        # get upper bound -
        forward_bound_value = reaction_wrapper.forward
        if (forward_bound_value == "inf")
            push!(ub_array,1000.0)
        else
            push!(ub_array,parse(Float64,forward_bound_value))
        end
    end
    cobra_dictionary["ub"] = ub_array

    # build metabolite name array -
    metabolite_symbol_order_array = problem_dictionary["metabolite_symbol_order_array"]
    list_of_metabolite_wrappers = problem_dictionary["list_of_metabolites"]
    met_names_array = Array{String,1}()
    met_kegg_id_array = Array{String,1}()
    met_symbol_array = Array{String,1}()
    for metabolite_symbol in metabolite_symbol_order_array

        # do we have this metabolite wrapper?
        if (haskey(list_of_metabolite_wrappers,metabolite_symbol) == true)

            # get the met name -
            met_wrapper = list_of_metabolite_wrappers[metabolite_symbol]
            met_symbol = met_wrapper.symbol
            met_name = met_wrapper.name
            met_kegg_id = met_wrapper.kegg_id

            # cache -
            push!(met_symbol_array, met_symbol)
            push!(met_names_array,met_name)
            push!(met_kegg_id_array,met_kegg_id)
        else

            # cache -
            push!(met_symbol_array,"")
            push!(met_names_array,"")
            push!(met_kegg_id_array,"")
        end
    end
    cobra_dictionary["metNames"] = met_names_array
    cobra_dictionary["metKEGGID"] = met_kegg_id_array
    cobra_dictionary["mets"] = met_symbol_array

    # S: generate the stoichiometric matrix -
    S = build_stoichiometric_matrix(met_symbol_array,reaction_wrapper_array)
    cobra_dictionary["S"] = S

    # c: objective vector - default is 0's
    # how many reactions do we have?
    (number_of_species,number_of_reactions) = size(S)
    c_vector = zeros(number_of_reactions)
    cobra_dictionary["c"] = c_vector

    # rxnGeneMat - reaction to gene mapping array -
    number_of_genes = length(gene_order_array)
    rxnGeneMat = zeros(number_of_reactions, number_of_genes)
    for (row_index, reaction_tag) in enumerate(reaction_order_array)
        for (col_index, gene_symbol) in enumerate(gene_order_array)

            # ok, get the reation-gene map for this reaction -
            if (haskey(list_of_genes, reaction_tag) == true)

                # grab the genes -
                gene_set = list_of_genes[reaction_tag].genes

                # does this set contain the gene?
                if (in(gene_symbol,gene_set) == true)
                    rxnGeneMat[row_index,col_index] = 1.0
                end
            end
        end
    end
    cobra_dictionary["rxnGeneMat"] = rxnGeneMat

    # return -
    return cobra_dictionary
end

function export_reaction_mapping_file(cobra_dictionary::Dict{String,Any}, path_to_mapping_file::String)

    # Check the dictionary -
    is_cobra_dictionary_ok(cobra_dictionary)

    # initalize -
    mapping_buffer = String[]

    # add header text to mapping file -
    filename = splitdir(path_to_mapping_file)[2]
    +(mapping_buffer,"// ------------------------------------------------------------- //")
    +(mapping_buffer,"// $(filename)")
    +(mapping_buffer,"// GENERATED BY: CBModelTools")
    +(mapping_buffer,"// GENERATED ON: $(Dates.now())")
    +(mapping_buffer,"// SOURCE: https://github.com/varnerlab/CBModelTools")
    +(mapping_buffer,"//")
    +(mapping_buffer,"// Reaction->reaction name map")
    +(mapping_buffer,"// record: reacton_tag=reaction_name")
    +(mapping_buffer,"// reaction_tag: from the rxns field of the COBRA mat file")
    +(mapping_buffer,"// reaction_name: from the rxnNames field of the COBRA mat file")
    +(mapping_buffer,"//")
    +(mapping_buffer,"// The order of the reactions in this file is used to order the")
    +(mapping_buffer,"// columns of the stoichiometric matrix.")
    +(mapping_buffer,"// ------------------------------------------------------------ //")

    list_of_reactions = cobra_dictionary["rxns"]
    list_of_reaction_names = cobra_dictionary["rxnNames"]

    # process each gene -
    for (index,reaction_tag) in enumerate(list_of_reactions)

        # init -
        line = ""

        # get content -
        reaction_name = list_of_reaction_names[index]

        # build line -
        line *= "$(reaction_tag)=$(reaction_name)"

        # cache -
        +(mapping_buffer, line)
    end

    # Write out the mapping file
    write_file_to_path(path_to_mapping_file,mapping_buffer)
end

function export_gene_order_file(cobra_dictionary::Dict{String,Any}, kegg_organism_code::String, path_to_mapping_file::String)

    # Check organism code -
    is_kegg_organism_code_ok(kegg_organism_code)

    # Check the dictionary -
    is_cobra_dictionary_ok(cobra_dictionary)

    # initalize -
    mapping_buffer = String[]
    list_of_genes = cobra_dictionary["genes"]

    # process each gene -
    for gene_id in list_of_genes

        # make the KEGG gene id -
        kegg_gene_id = "$(kegg_organism_code):$(gene_id)"

        # cache -
        +(mapping_buffer, kegg_gene_id)
    end

    # Write out the vff file -
    open("$(path_to_mapping_file)", "w") do f

        for line_item in mapping_buffer
            write(f,"$(line_item)\n")
        end
    end
end

function export_reaction_tag_to_gene_mapping_file(cobra_dictionary::Dict{String,Any}, kegg_organism_code::String, path_to_mapping_file::String)

    # Check organism code -
    is_kegg_organism_code_ok(kegg_organism_code)

    # Check the dictionary -
    is_cobra_dictionary_ok(cobra_dictionary)

    # initalize -
    mapping_buffer = String[]

    # add header text to mapping file -
    filename = splitdir(path_to_mapping_file)[2]
    +(mapping_buffer,"// --------------------------------------------------------- //")
    +(mapping_buffer,"// $(filename)")
    +(mapping_buffer,"// GENERATED BY: CBModelTools")
    +(mapping_buffer,"// GENERATED ON: $(Dates.now())")
    +(mapping_buffer,"// SOURCE: https://github.com/varnerlab/CBModelTools")
    +(mapping_buffer,"")
    +(mapping_buffer,"// Reaction->gene mapping - ")
    +(mapping_buffer,"// record: reacton_tag={gene_id}")
    +(mapping_buffer,"// reaction_tag: from the rxns field of the COBRA mat file")
    +(mapping_buffer,"// gene_id: (kegg organism code):(gene location code)")
    +(mapping_buffer,"// --------------------------------------------------------- //")

    # get the list of reaction tags, and genes from the cobra dictionary -
    list_of_reaction_tags = cobra_dictionary["rxns"]
    list_of_genes = cobra_dictionary["genes"]

    # get the rxn gene mapping matrix -
    RGM = Matrix(cobra_dictionary["rxnGeneMat"])

    # What is the size of the system?
    (number_of_reactions, number_of_genes) = size(RGM)

    # main loop -
    for reaction_index = 1:number_of_reactions

        # what is the tag for this reaction?
        reaction_tag = list_of_reaction_tags[reaction_index]

        # initialize -
        record = ""
        record *= "$(reaction_tag)="

        # does this reaction have associated genes?
        idx_genes = findall(x->x==1,RGM[reaction_index,:])
        local_gene_set = Set{String}()
        if (!isempty(idx_genes))

            # ok, we have genes, grab them -
            gene_id_array = list_of_genes[idx_genes]

            # process each gene -
            for gene_id in gene_id_array

                # make the KEGG gene id -
                kegg_gene_id = gene_id
                if (occursin(".",gene_id) == true)
                    # need to cutoff the trailing *.1
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id[1:end-2])"
                else
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id)"
                end

                # cache -
                push!(local_gene_set,kegg_gene_id)
            end
        end

        if (isempty(local_gene_set) == false)

            # make an ec record -
            gene_record = ""
            for gene_number in local_gene_set
                gene_record*="$(gene_number),"
            end

            # remove trailing ,
            gene_record = gene_record[1:end-1]

            # add the record -
            record*="$(gene_record)"

            # push into buffer and go around again -
            +(mapping_buffer,record)
        end
    end

    # Write out the file -
    write_file_to_path(path_to_mapping_file, mapping_buffer)
end

function export_reaction_tag_to_ec_mapping_file(cobra_dictionary::Dict{String,Any}, kegg_organism_code::String, path_to_mapping_file::String)

    # Check organism code -
    is_kegg_organism_code_ok(kegg_organism_code)

    # Check the dictionary -
    is_cobra_dictionary_ok(cobra_dictionary)

    # initalize -
    mapping_buffer = String[]

    # add header text to mapping file -
    filename = splitdir(path_to_mapping_file)[2]
    +(mapping_buffer,"// --------------------------------------------------------- //")
    +(mapping_buffer,"// $(filename)")
    +(mapping_buffer,"// GENERATED BY: CBModelTools")
    +(mapping_buffer,"// GENERATED ON: $(Dates.now())")
    +(mapping_buffer,"// SOURCE: https://github.com/varnerlab/CBModelTools")
    +(mapping_buffer,"//")
    +(mapping_buffer,"// Reaction->ecnumber mapping - ")
    +(mapping_buffer,"// record: reacton_tag={ecnumber}")
    +(mapping_buffer,"// reaction_tag: from the rxns field of the COBRA mat file")
    +(mapping_buffer,"// ec_number: possible ecnumber estimated from KEGG")
    +(mapping_buffer,"// --------------------------------------------------------- //")

    # get the list of reaction tags, and genes from the cobra dictionary -
    list_of_reaction_tags = cobra_dictionary["rxns"]
    list_of_genes = cobra_dictionary["genes"]

    # get the rxn gene mapping matrix -
    RGM = Matrix(cobra_dictionary["rxnGeneMat"])

    # What is the size of the system?
    (number_of_reactions, number_of_genes) = size(RGM)

    # Declare a progress meter for user feedback -
    p = Progress(number_of_reactions,color=:yellow)

    # main loop -
    for reaction_index = 1:number_of_reactions

        # what is the tag for this reaction?
        reaction_tag = list_of_reaction_tags[reaction_index]

        # user message -
        msg = "Starting $(reaction_tag) ($(reaction_index) of $(number_of_reactions))"

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:status,msg)])

        # initialize -
        record = ""
        record *= "$(reaction_tag)="

        # does this reaction have associated genes?
        idx_genes = findall(x->x==1,RGM[reaction_index,:])

        # init -
        ec_record_set = Set{String}()
        if (!isempty(idx_genes))

            # ok, we have genes, grab them -
            gene_id_array = list_of_genes[idx_genes]

            # process each gene -
            for gene_id in gene_id_array

                # make the KEGG gene id -
                kegg_gene_id = gene_id
                if (occursin(".",gene_id) == true)
                    # need to cutoff the trailing *.1
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id[1:end-2])"
                else
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id)"
                end

                # lookup -
                local_ec_record_set = lookup_gene_ec_mapping_record(kegg_gene_id)
                if local_ec_record_set != nothing

                    # push into ec_record_set -
                    for ec_number in local_ec_record_set
                        push!(ec_record_set,ec_number)
                    end
                end
            end
        end

        if (isempty(ec_record_set) == false)

            # make an ec record -
            ec_record = ""
            for ec_number in ec_record_set
                ec_record*="$(ec_number),"
            end

            # remove trailing ,
            ec_record = ec_record[1:end-1]

            # add the record -
            record*="$(ec_record)"

            # push into buffer and go around again -
            +(mapping_buffer,record)
        end
    end

    # Write out the file -
    write_file_to_path(path_to_mapping_file,mapping_buffer)
end
# ------------------------------------------------------------------------------ #

max_stars_repo_name: varnerlab/CBModelTools

max_issues_repo_name: varnerlab/CBModelTools

max_forks_repo_name: varnerlab/CBModelTools


text:
module BenchLinAlg

using BenchmarkTools
using SparseArrays, LinearAlgebra

# include("../src/linalg.jl")
# include(joinpath(@__DIR__, "../test/setup_dqmc.jl"))
using QMC

suite = BenchmarkGroup()

mc = mc_from_inxml(joinpath(@__DIR__, "../test/parameters/O3_generic_small_system.in.xml"))
dtype = ComplexF64
n = size(mc.s.greens,1)
A, B, C = rand(dtype, n,n), rand(dtype, n,n), rand(dtype, n,n)
S = sparse(Diagonal(rand(dtype, n)))

U1,D1,T1 = decompose_udt(mc.s.greens)
U2,D2,T2 = decompose_udt(mc.s.greens)
U3,D3,T3 = decompose_udt(mc.s.greens)
res = copy(mc.s.greens)

suite["decompose_udt"] = @benchmarkable decompose_udt($A)
suite["decompose_udt!"] = @benchmarkable decompose_udt!($A, $D1)

suite["Matrix*Sparse (mul!)"] = @benchmarkable mul!($C, $A, $S)
suite["Sparse*Matrix (mul!)"] = @benchmarkable mul!($C, $S, $A)

suite["inv_udt"] = @benchmarkable inv_udt($U1, $D1, $T1)
suite["inv_udt!"] = @benchmarkable inv_udt!($C, $U1, $D1, $T1)

suite["inv_one_plus_two_udts!"] = @benchmarkable inv_one_plus_two_udts!($mc,$U1,$D1,$T1,$U2,$D2,$T2,$U3,$D3,$T3)
suite["inv_one_plus_two_udts!"] = @benchmarkable inv_one_plus_two_udts!($mc,$res,$U2,$D2,$T2,$U3,$D3,$T3)
suite["inv_one_plus_udt_scalettar!"] = @benchmarkable inv_one_plus_udt_scalettar!($mc,$res,$U3,$D3,$T3)
suite["inv_sum_udts_scalettar!"] = @benchmarkable inv_sum_udts_scalettar!($mc,$res,$U1,$D1,$T1,$U2,$D2,$T2)

end  # module
BenchLinAlg.suite

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
#! /usr/bin/env julia

# Copyright 2016 John Hanley.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# The software is provided "AS IS", without warranty of any kind, express or
# implied, including but not limited to the warranties of merchantability,
# fitness for a particular purpose and noninfringement. In no event shall the
# authors or copyright holders be liable for any claim, damages or other
# liability, whether in an action of contract, tort or otherwise, arising from,
# out of or in connection with the software or the use or other dealings in the
# software.

# usage:
#     $ ./class_deps.jl *.java

using PyCall
@pyimport pydot

import_re = r"^import(\s+static|)\s+([\w\.]*\w)"

function get_edges(src_files)
    Task() do
        for fspec in src_files
            fspec1 = replace(basename(fspec), ".java", "")
            open(fspec) do fin
                for line in readlines(fin)
                    m = match(import_re, line)
                    if m != nothing
                        produce(fspec1, m[2])
                    end
                end
            end
        end
    end
end

macro p_str(s) s end

ignore_re = Regex(join([
    p"^(javax?|android|org\.slf4j|com\.google)\.",
    p"^org.apache.commons.lang3.(StringUtils|builder.HashCodeBuilder|time.FastDateFormat)",
    p"\.util\.(Util|CardValue|Clock|LoggerBg|STException)$",
    p"APP_VERSION|\.R$"
], "|"))

function get_filtered_edges(src_files)
    Task() do
        for (src, dst) in get_edges(src_files)
            if !ismatch(ignore_re, dst)
                produce(src, dst)
            end
        end
    end
end

function abbrev(class)
    replace(class, r".*\.", "")
end

function report(out_file, src_files, ignores)
    g = pydot.Dot(rankdir="LR")

    for (src, dst) in get_filtered_edges(src_files)
        println(src, " ", dst)
        e = pydot.Edge(src, abbrev(dst))
        g[:add_edge](e)
    end
    g[:write](out_file)
end

report("/tmp/deps.dot", ARGS, Set())


text:
module FourierFlows

export
  # Helper variables and macros for determining if machine is CUDA-enabled.
  @has_cuda,

  Device,
  CPU,
  GPU,
  ArrayType,

  cxtype,
  fltype,
  innereltype,

  AbstractVars,
  AbstractParams,

  AbstractGrid,
  ZeroDGrid,
  OneDGrid,
  TwoDGrid,
  ThreeDGrid,
  dealias!,
  gridpoints,

  Diagnostic,
  resize!,
  update!,
  increment!,

  Output,
  saveoutput,
  saveproblem,
  groupsize,
  savediagnostic,

  @zeros,
  @devzeros,
  @createarrays,
  @superzeros,
  devzeros,
  superzeros,
  supersize,

  TimeStepper,
  ForwardEulerTimeStepper,
  FilteredForwardEulerTimeStepper,
  RK4TimeStepper,
  FilteredRK4TimeStepper,
  ETDRK4TimeStepper,
  FilteredETDRK4TimeStepper,
  AB3TimeStepper,
  FilteredAB3TimeStepper,
  stepforward!

using
  FFTW,
  JLD2,
  Statistics,
  Interpolations,
  CUDAapi,
  Requires

import Base: resize!, getindex, setindex!, lastindex, push!, append!

using Base: fieldnames
using AbstractFFTs: fftfreq, rfftfreq

abstract type AbstractGrid{T, Ta} end
abstract type AbstractTimeStepper{T} end
abstract type AbstractParams end
abstract type AbstractVars end
abstract type AbstractDiagnostic end

abstract type Device end
struct CPU <: Device end
struct GPU <: Device end

# The main show
include("problem.jl")
include("domains.jl")
include("utils.jl")
include("diagnostics.jl")
include("output.jl")
include("timesteppers.jl")

# Physics
include("diffusion.jl")
include("diffusion2D.jl")
# Import CUDA utilities if cuda is detected.
if has_cuda()
  try
    using CuArrays # we have CUDA, so this should not fail
  catch ex
    # something is wrong with the user's set-up (or there's a bug in CuArrays)
    @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex, catch_backtrace())
  end
end

macro has_cuda(ex)
  return has_cuda() ? :($(esc(ex))) : :(nothing)
end

function __init__()
  @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("CuFourierFlows.jl")
end

end # module


text:
###############################################################################
#
#   fq_rel_series.jl : Power series over flint finite fields
#
###############################################################################

export fq_rel_series, FqRelSeriesRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::fq_rel_series)
   val = pol_length(a) + valuation(a) - 1
   val < 0 && throw(DomainError())
   z = fq_rel_series(base_ring(a), Array{fq}(0), 0, val, val)
   z.parent = parent(a)
   return z
end

elem_type(::Type{FqRelSeriesRing}) = fq_rel_series

parent_type(::Type{fq_rel_series}) = FqRelSeriesRing

base_ring(R::FqRelSeriesRing) = R.base_ring

var(a::FqRelSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::FqRelSeriesRing) = R.prec_max

function normalise(a::fq_rel_series, len::Int)
   ctx = base_ring(a)
   if len > 0
      c = base_ring(a)()
      ccall((:fq_poly_get_coeff, :libflint), Void,
         (Ref{fq}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
                                                         c, a, len - 1, ctx)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         ccall((:fq_poly_get_coeff, :libflint), Void,
            (Ref{fq}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
                                                         c, a, len - 1, ctx)
      end
   end
   return len
end

function pol_length(x::fq_rel_series)
   return ccall((:fq_poly_length, :libflint), Int,
                (Ref{fq_rel_series}, Ref{FqFiniteField}), x, base_ring(x))
end

precision(x::fq_rel_series) = x.prec

function polcoeff(x::fq_rel_series, n::Int)
   z = base_ring(x)()
   if n < 0
      return z
   end
   ccall((:fq_poly_get_coeff, :libflint), Void,
         (Ref{fq}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
                                                      z, x, n, base_ring(x))
   return z
end

zero(R::FqRelSeriesRing) = R(0)

one(R::FqRelSeriesRing) = R(1)

function gen(R::FqRelSeriesRing)
   z = fq_rel_series(base_ring(R), [base_ring(R)(1)], 1, max_precision(R) + 1, 1)
   z.parent = R
   return z
end

function deepcopy_internal(a::fq_rel_series, dict::ObjectIdDict)
   z = fq_rel_series(base_ring(a), a)
   z.prec = a.prec
   z.val = a.val
   z.parent = parent(a)
   return z
end

function renormalize!(z::fq_rel_series)
   i = 0
   zlen = pol_length(z)
   zval = valuation(z)
   zprec = precision(z)
   while i < zlen && iszero(polcoeff(z, i))
      i += 1
   end
   z.prec = zprec
   if i == zlen
      z.val = zprec
   else
      z.val = zval + i
      ccall((:fq_poly_shift_right, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
                                                      z, z, i, base_ring(z))
   end
   return nothing
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::FqRelSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{fq_rel_series}) = show_minus_one(fq)

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::fq_rel_series)
   z = parent(x)()
   ccall((:fq_poly_neg, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series}, Ref{FqFiniteField}),
               z, x, base_ring(x))
   z.prec = x.prec
   z.val = x.val
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fq_rel_series, b::fq_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   z = parent(a)()
   ctx = base_ring(a)
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fq_poly_set_trunc, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, b, max(0, lenz - b.val + a.val), ctx)
      ccall((:fq_poly_shift_left, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, z, b.val - a.val, ctx)
      ccall((:fq_poly_add_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, z, a, lenz, ctx)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fq_poly_set_trunc, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, a, max(0, lenz - a.val + b.val), ctx)
      ccall((:fq_poly_shift_left, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, z, a.val - b.val, ctx)
      ccall((:fq_poly_add_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, z, b, lenz, ctx)
   else
      lenz = max(lena, lenb)
      ccall((:fq_poly_add_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, ctx)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function -(a::fq_rel_series, b::fq_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   lenz = max(lena, lenb)
   z = parent(a)()
   ctx = base_ring(a)
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fq_poly_set_trunc, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, b, max(0, lenz - b.val + a.val), ctx)
      ccall((:fq_poly_shift_left, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, z, b.val - a.val, ctx)
      ccall((:fq_poly_neg, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Ref{FqFiniteField}),
            z, z, ctx)
      ccall((:fq_poly_add_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, z, a, lenz, ctx)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fq_poly_set_trunc, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, a, max(0, lenz - a.val + b.val), ctx)
      ccall((:fq_poly_shift_left, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, z, a.val - b.val, ctx)
      ccall((:fq_poly_sub_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, z, b, lenz, ctx)
   else
      lenz = max(lena, lenb)
      ccall((:fq_poly_sub_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, ctx)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function *(a::fq_rel_series, b::fq_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z = parent(a)()
   z.val = a.val + b.val
   z.prec = prec + z.val
   if lena == 0 || lenb == 0
      return z
   end
   lenz = min(lena + lenb - 1, prec)
   ccall((:fq_poly_mullow, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, base_ring(a))
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::fq, y::fq_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:fq_poly_scalar_mul_fq, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq}, Ref{FqFiniteField}),
               z, y, x, base_ring(y))
   return z
end

*(x::fq_rel_series, y::fq) = y * x

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fq_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   z = fq_rel_series(base_ring(x), x)
   z.prec = x.prec + len
   z.val = x.val + len
   z.parent = parent(x)
   return z
end

function shift_right(x::fq_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   xval = valuation(x)
   z = parent(x)()
   if len >= xlen + xval
      z.prec = max(0, x.prec - len)
      z.val = max(0, x.prec - len)
   else
      z.prec = max(0, x.prec - len)
      z.val = max(0, xval - len)
      zlen = min(xlen + xval - len, xlen)
      ccall((:fq_poly_shift_right, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Int, Ref{FqFiniteField}),
               z, x, xlen - zlen, base_ring(x))
      renormalize!(z)
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::fq_rel_series, prec::Int)
   prec < 0 && throw(DomainError())
   xlen = pol_length(x)
   xprec = precision(x)
   xval = valuation(x)
   if xprec + xval <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   if prec <= xval
      z = parent(x)()
      z.val = prec
      z.prec = prec
   else
      z.val = xval
      ccall((:fq_poly_set_trunc, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Int, Ref{FqFiniteField}),
               z, x, min(prec - xval, xlen), base_ring(x))
   end
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fq_rel_series, b::Int)
   b < 0 && throw(DomainError())
   if isgen(a)
      z = parent(a)()
      z = setcoeff!(z, 0, base_ring(a)(1))
      z.prec = a.prec + b - 1
      z.val = b
   elseif pol_length(a) == 0
      z = parent(a)()
      z.prec = b*valuation(a)
      z.val = b*valuation(a)
   elseif pol_length(a) == 1
      return parent(a)([polcoeff(a, 0)^b], 1,
                           (b - 1)*valuation(a) + precision(a), b*valuation(a))
   elseif b == 0
      return one(parent(a))
   else
      bit = ~((~UInt(0)) >> 1)
      while (UInt(bit) & b) == 0
         bit >>= 1
      end
      z = a
      bit >>= 1
      while bit != 0
         z = z*z
         if (UInt(bit) & b) != 0
            z *= a
         end
         bit >>= 1
      end
   end
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fq_rel_series, y::fq_rel_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)
   if prec <= x.val && prec <= y.val
      return true
   end
   if x.val != y.val
      return false
   end
   xlen = normalise(x, min(pol_length(x), prec - x.val))
   ylen = normalise(y, min(pol_length(y), prec - y.val))
   if xlen != ylen
      return false
   end
   return Bool(ccall((:fq_poly_equal_trunc, :libflint), Cint,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Int, Ref{FqFiniteField}),
               x, y, xlen, base_ring(x)))
end

function isequal(x::fq_rel_series, y::fq_rel_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || x.val != y.val || pol_length(x) != pol_length(y)
      return false
   end
   return Bool(ccall((:fq_poly_equal, :libflint), Cint,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Int, Ref{FqFiniteField}),
               x, y, pol_length(x), base_ring(x)))
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fq_rel_series, y::fq_rel_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   yval = valuation(y)
   xval = valuation(x)
   if yval != 0
      if xval >= yval
         x = shift_right(x, yval)
         y = shift_right(y, yval)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec - x.val, y.prec - y.val)
   z = parent(x)()
   z.val = xval - yval
   z.prec = prec + z.val
   if prec != 0
      ccall((:fq_poly_div_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, x, y, prec, base_ring(x))
   end
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fq_rel_series, y::fq)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.prec = x.prec
   z.val = x.val
   ccall((:fq_poly_scalar_div_fq, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq}, Ref{FqFiniteField}),
               z, x, y, base_ring(x))
   return z
end

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fq_rel_series)
   iszero(a) && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ainv.val = 0
   ccall((:fq_poly_inv_series, :libflint), Void,
         (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               ainv, a, a.prec, base_ring(a))
   return ainv
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(x::fq_rel_series)
  ccall((:fq_poly_zero, :libflint), Void,
                   (Ref{fq_rel_series}, Ref{FqFiniteField}), x, base_ring(x))
  x.prec = parent(x).prec_max
  return x
end

function fit!(z::fq_rel_series, n::Int)
   ccall((:fq_poly_fit_length, :libflint), Void,
         (Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
         z, n, base_ring(z))
   return nothing
end

function setcoeff!(z::fq_rel_series, n::Int, x::fmpz)
   ccall((:fq_poly_set_coeff_fmpz, :libflint), Void,
                (Ref{fq_rel_series}, Int, Ref{fmpz}, Ref{FqFiniteField}),
               z, n, x, base_ring(z))
   return z
end

function setcoeff!(z::fq_rel_series, n::Int, x::fq)
   ccall((:fq_poly_set_coeff, :libflint), Void,
                (Ref{fq_rel_series}, Int, Ref{fq}, Ref{FqFiniteField}),
               z, n, x, base_ring(z))
   return z
end

function mul!(z::fq_rel_series, a::fq_rel_series, b::fq_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z.val = a.val + b.val
   z.prec = prec + c.val
   lenz = min(lena + lenb - 1, prec)
   if lena <= 0 || lenb <= 0
      lenz = 0
   end
   ccall((:fq_poly_mullow, :libflint), Void,
         (Ref{fq_rel_series}, Ref{fq_rel_series},
          Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, base_ring(z))
   return z
end

function addeq!(a::fq_rel_series, b::fq_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   ctx = base_ring(a)
   if a.val < b.val
      z = fq_rel_series(base_ring(a))
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fq_poly_set_trunc, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, b, max(0, lenz - b.val + a.val), ctx)
      ccall((:fq_poly_shift_left, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            z, z, b.val - a.val, ctx)
      ccall((:fq_poly_add_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               a, a, z, lenz, ctx)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fq_poly_truncate, :libflint), Void,
            (Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            a, max(0, lenz - a.val + b.val), ctx)
      ccall((:fq_poly_shift_left, :libflint), Void,
            (Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
            a, a, a.val - b.val, ctx)
      ccall((:fq_poly_add_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               a, a, b, lenz, ctx)
   else
      lenz = max(lena, lenb)
      ccall((:fq_poly_add_series, :libflint), Void,
                (Ref{fq_rel_series}, Ref{fq_rel_series},
                 Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               a, a, b, lenz, ctx)
   end
   a.prec = prec
   a.val = val
   renormalize!(a)
   return a
end

function add!(c::fq_rel_series, a::fq_rel_series, b::fq_rel_series)
   ctx = base_ring(a)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenc = max(lena, lenb)
   c.prec = prec
   ccall((:fq_poly_add_series, :libflint), Void,
     (Ref{fq_rel_series}, Ref{fq_rel_series}, Ref{fq_rel_series}, Int, Ref{FqFiniteField}),
               c, a, b, lenc, ctx)
   return c
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{fq_rel_series}, ::Type{T}) where {T <: Integer} = fq_rel_series

promote_rule(::Type{fq_rel_series}, ::Type{fmpz}) = fq_rel_series

promote_rule(::Type{fq_rel_series}, ::Type{fq}) = fq_rel_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FqRelSeriesRing)()
   ctx = base_ring(a)
   z = fq_rel_series(ctx)
   z.prec = a.prec_max
   z.val = a.prec_max
   z.parent = a
   return z
end

function (a::FqRelSeriesRing)(b::Integer)
   ctx = base_ring(a)
   if b == 0
      z = fq_rel_series(ctx)
      z.prec = a.prec_max
   else
      z = fq_rel_series(ctx, [ctx(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FqRelSeriesRing)(b::fmpz)
   ctx = base_ring(a)
   if iszero(b)
      z = fq_rel_series(ctx)
      z.prec = a.prec_max
   else
      z = fq_rel_series(ctx, [ctx(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FqRelSeriesRing)(b::fq)
   ctx = base_ring(a)
   if iszero(b)
      z = fq_rel_series(ctx)
      z.prec = a.prec_max
   else
      z = fq_rel_series(ctx, [b], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FqRelSeriesRing)(b::fq_rel_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::FqRelSeriesRing)(b::Array{fq, 1}, len::Int, prec::Int, val::Int)
   ctx = base_ring(a)
   z = fq_rel_series(ctx, b, len, prec, val)
   z.parent = a
   return z
end


text:
# have Petsc report errors to Julia
"""
  Error handler registered with PETSc.  This function prints the information
  Petsc supplies about the error, a Julia stack trace, and then invokes the
  `PetscTraceBackErrorHandler` to print the Petsc stack trace.
"""
function error_handler(comm::comm_type, line::Cint, func::Ptr{UInt8}, file::Ptr{UInt8}, n::PetscErrorCode, p::PetscErrorType, mess::Ptr{UInt8}, ctx::Ptr{Void})
# receives the error call from Petsc

  func_string = unsafe_string(func)
  file_string = unsafe_string(file)

  if p == PETSC_ERROR_INITIAL
    p_string = "Initial Petsc Error"
  elseif p == PETSC_ERROR_REPEAT
    p_string = "Repeat Petsc Error"
  elseif  p == PETSC_ERROR_IN_CXX
    p_string == "CXX Petsc Error"
  else
    pstring = "Unknown PetscErrorType"
  end

  mess_string = unsafe_string(mess)
  print_with_color(:red, STDERR, string("\n### ERROR: PETSc Internal Error ###\n"))
  print_with_color(:red, STDERR, string("Error in function ", func_string, ", file ", file_string, "\n"))
  print_with_color(:red, STDERR, string("Error number: ", n, ", Error type: ", p_string, "\n"))
  print_with_color(:red, STDERR, string("Error message: \n", mess_string, "\n"))
  print_with_color(:red, STDERR, string("### Error Message Finished ###\n"))

  bt = backtrace()
  s = sprint(io->Base.show_backtrace(io, bt))
  print_with_color(:red, STDERR, string("Julia backtrace: ", s))
  print(STDERR, "\n\n")

  # invoking Petsc error handler
  PetscTraceBackErrorHandler(comm, line, func, file, n, p, mess, ctx)

  println(STDERR, "finished invoking Petsc error handler")

  return PetscErrorCode(0)

end

function PetscTraceBackErrorHandler(arg1::comm_type,arg2::Cint,arg3::Ptr{UInt8},
           arg4::Ptr{UInt8},arg5::PetscErrorCode,arg6::PetscErrorType,
           arg7::Ptr{UInt8},arg8::Ptr{Void})

    ccall((:PetscTraceBackErrorHandler,petsc),PetscErrorCode,
      (comm_type,Cint,Ptr{UInt8},Ptr{UInt8},PetscErrorCode,PetscErrorType,Ptr{UInt8},Ptr{Void}),
      arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8)
end

# tell Petsc abut the Error handler
cfunc = cfunction(error_handler, PetscErrorCode, (comm_type, Int32, Ptr{UInt8}, Ptr{UInt8}, PetscErrorCode, PetscErrorType, Ptr{UInt8}, Ptr{Void}) )
ctx = C_NULL
ierr = ccall( (:PetscPushErrorHandler, petsc), PetscErrorCode, (Ptr{Void}, Ptr{Void}), cfunc, ctx)

function chkerrq(i::PetscErrorCode)
  if i != 0
    error("Petsc Error Code: $i")
  end

  return nothing
end


text:
"""
    ScalarCartesianSet{N, TS}

Represents the Cartesian product `S` of several scalar sets.
`S = S_1 × S_2 × S_3 ... × S_N`, such that a vector `x ∈ S ⇔ x[i] ∈ S[i] ∀ i`.
The dimension of the set is the number of elementary scalar sets.
"""
struct ScalarCartesianSet{N, TS <: NTuple{N, MOI.AbstractScalarSet}} <: MOI.AbstractVectorSet
    sets::TS
end

ScalarCartesianSet(sets...) = ScalarCartesianSet(sets)

MOI.dimension(::ScalarCartesianSet{N}) where {N} = N

# leverage bridges here?
function MOI.add_constraint(model::MOI.ModelLike, vec::MOI.VectorOfVariables, set::ScalarCartesianSet)
    cons_indices = [
        MOI.add_constraint(model, x[i], set.sets[i]) for i in eachindex(set.sets)
    ]
    # TODO index type?
end


text:
if nworkers()<2
	addprocs(3-nworkers())
end

using jInv.Mesh
using jInv.LinearSolvers
using jInv.ForwardShare
using ForwardHelmholtz
using FWI
using Base.Test

# get a small mesh
domain = [0 2.0 0 2.2]
n      = [12;8;]
M      = getRegularMesh(domain,n)
omega  = 0.1;

# sources / receivers on top edge
idx    = reshape(collect(1:prod(n+1)),tuple(n+1...))
ib     = idx[:,1];
I      = speye(prod(n+1))
Q      = I[:,vec(ib)]
R      = copy(Q)

# get param without parallelizaion
ABLpad = 2;
gamma = getABL(M,true,ones(Int64,M.dim)*ABLpad,8*pi);
attenuation = 0.01*4*pi;
gamma += attenuation; # adding Attenuation.
gamma = gamma[:];
Ainv = getJuliaSolver();

pFor = getFWIparam(omega, one(Complex128),gamma,Q,R,M,Ainv,[workers()[1]])[1];

m0   = rand(tuple(n+1...))+1.0
dobs, = getData(vec(m0),fetch(pFor[1]));

# parallelize over sources
pForp,continuationDivision,SourcesSubInd = getFWIparam(omega, one(Complex128),gamma,Q,R,M,Ainv,workers())
dobsRF, = getData(vec(m0),pForp)
Dobs = zeros(eltype(dobs),size(dobs));
for k=1:length(dobsRF)
	dobsk = fetch(dobsRF[k]);
	Dobs[:,SourcesSubInd[k]] = dobsk;
end

@test vecnorm(Dobs-dobs)/vecnorm(dobs) < 1e-8


text:
export read_from, file_import, write_to

"""
    read_from(infile::String; kwargs...)
Wraps `delimited()`, `genepop()`, `bcf()`, and `vcf()` to read a file in as a `PopObj`. File type is
inferred from the file extension (case insensitive): \n

| File Format         | Extensions             | Docstring     |
| :------------------ | :--------------------- | :------------ |
| delimited           | `.csv`, `.txt`, `.tsv` | `?delimited`  |
| genepop             | `.gen`, `.genepop`     | `?genepop`    |
| variant call format | `.vcf`, `.bcf`         | `?vcf, ?bcf`  |

This function uses the same keyword arguments (and defaults) as the file importing
functions it wraps; please see their respective docstrings in the Julia help console.
(e.g. `?genepop`) for specific usage details. Use the alias function `file_import`
interchangeably if you prefer the explicit name instead.

## Examples
```
read_from("cavernous_assfish.gen", digits = 3)

file_import("bos_tauros.csv", silent = true)

read_from("juglans_nigra.vcf")
```
"""
function read_from(infile::String; kwargs...)
    ext = split(infile, ".")[end] |> lowercase
    if ext in ["gen", "genepop"]
        return genepop(infile;kwargs...)

    elseif ext in ["csv", "txt", "tsv"]
        return delimited(infile; kwargs...)

    elseif ext in ["vcf", "bcf"]
        ext == "vcf" && return vcf(infile)
        ext == "bcf" && return bcf(infile)

    else
        @error "File type not recognized by filename extension \n delimited: .csv | .tsv | .txt \n genepop: .gen | .genepop \n variant call formant: .bcf | .vcf"
    end
end

const file_import = read_from

"""
    write_to(data::PopData; filename::String, kwargs...)
Writes `PopData` to a specified file type inferred from the extension of `filename = ` (case insensitive). Additional keyword
arguments `kwargs...` are specific to the intended file type, and are listed in the docstrings of the specific
file writer with the format `?popdata2filetype`. For example, to find the appropriate keywords for a conversion
to Genepop format, call up the docstring to `popdata2genepop` with `?popdata2genepop`.

| File Format | Extensions             | Docstring          |
| :---------- | :--------------------- | :----------------- |
| genepop     | '.gen', '.genepop'     | ?popdata2genepop   |
| JLD2        | '.jld2'                | ?popdata2jld2      |
| delimited   | '.csv', '.txt', '.tsv' | ?popdata2delimited |

### Example
```
cats = nancycats();
fewer_cats = omit(cats, name = samples(cats)[1:10]);
write_to(fewer_cats, filename = "filtered_nancycats.gen", digits = 3, format = "h")
```
"""
function write_to(data::PopData; filename::String, kwargs...)
    ext = split(filename, ".")[end] |> lowercase
    if ext in ["gen", "genepop"]
        popdata2genepop(data, filename = filename; kwargs...)
    elseif ext == "jld2"
        popdata2jdl2(data, filename = filename)
    elseif ext in ["csv", "txt", "tsv"]
        popdata2delimited(data, filename = filename; kwargs...)
    else
        @error "File type not recognized by filename extension. Please see the docstring"
    end
end

"""
    popdata2jdl2(data::PopData; filename::String)
Write PopData to a `JLD2` file.
"""
function popdata2jdl2(data::PopData; filename::String)
    JLD2.@save filename data
end


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))

using StaticArrays
using Distributions
using Random
ImproperPrior()
DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))

# Let's load-in the data
# -----------------------
#obs = open(joinpath(OUT_DIR, "autoreg50fo.dat")) do f
#    [map(x->parse(Float64, x), split(l, ' '))[[2,3,4,1]] for l in eachline(f)]
#end
#obs_time = collect(range(0.0, length(obs)-1, step=1))

# Let's load-in the partially observed data
# -----------------------
_obs = open(joinpath(OUT_DIR, "prokaryote_custom.dat")) do f
    [map(x->parse(Float64, x), split(l, ' '))[[2,3]] for (i,l) in enumerate(eachline(f)) if i != 1]
end
_obs_time = collect(range(1.0, length(_obs), step=1))

obs = vcat([[0.0,0.0]], _obs)
obs_time = vcat(0.0, _obs_time)

# let's start from the true values that generated the data
θ_init = [0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1]
K = 10.0
P˟ = Prokaryote(θ_init..., K)
x0 = ℝ{4}([8.0,8.0,8.0,5.0])

auxFlag = Val{:custom}()
start_v = @SVector[5.0, 5.0, 5.0, 5.0]
P̃ = [ProkaryoteAux(θ_init..., K, t₀, u, T, v, auxFlag, start_v) for (t₀, T, u, v)
     in zip(obs_time[1:end-1], obs_time[2:end], obs[1:end-1], obs[2:end])]

Σ = @SMatrix[1.0 0.0; 0.0 2.0]
L = @SMatrix[1.0 0.0 0.0 0.0;
             0.0 1.0 2.0 0.0]

model_setup = DiffusionSetup(Pˣ, P̃, PartObs())
set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs_vals,
                  obs_time)
set_imputation_grid!(model_setup, 1/100)
set_x0_prior!(model_setup,
              GsnStartingPt(x0, SArray{Tuple{4,4}, Float64}(I)), # prior over starting point
              x0)
set_auxiliary!(model_setup; skip_for_save=thin_path,
               adaptive_prop=Adaptation(x0, [0.7, 0.4, 0.2, 0.2, 0.2],
                                            [500, 500, 500, 500, 500],
                                        1))
initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(num_steps_for_change_pt))

set_transition_kernels!(setup,
                        [ RandomWalk(fill(1.0, 8), collect(1:8)) for i in 1:7],
                        0.9, true, [[1], [2], [3], [4], [5], [7], [8]],
                        (MetropolisHastingsUpdt(),
                         MetropolisHastingsUpdt(),
                         MetropolisHastingsUpdt(),
                         MetropolisHastingsUpdt(),
                         MetropolisHastingsUpdt(),
                         MetropolisHastingsUpdt(),
                         MetropolisHastingsUpdt(),
                         ),
                         Adaptation(x0,
                                    fill(0.0, 10),
                                    fill(100, 10),
                                    1)
                        )
set_priors!(setup,
            Priors((ImproperPosPrior{Val{1}}(),ImproperPosPrior{Val{2}}(),
            ImproperPosPrior{Val{3}}(),
                    ImproperPosPrior{Val{4}}(),
                    ImproperPosPrior{Val{5}}(),
                    ImproperPosPrior{Val{7}}(),
                    ImproperPosPrior{Val{8}}())),
            KnownStartingPt(x0)
            )
set_mcmc_params!(setup, 1*10^4, 1*10^2, 10^1, 10^0, 0,
                 (50, 0.1, 0.00001, 0.99999, 0.234, 50),
                 (50, 0.1, -999, 999, 0.234, 50, (1,2,3,4,5,6,7), (1,2,3,4,5,7,8)))
set_blocking!(setup)
set_solver!(setup, Vern7(), NoChangePt())
initialise!(eltype(x0), setup)

Random.seed!(4)
out, elapsed = @timeit mcmc(setup)
display(out.accpt_tracker)

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_paths(out; ylims=[(0,15), nothing, nothing, nothing])
plot_chains(out; truth=[0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1])


text:
# # Saving and loading a mesh

# In this tutorial, we will learn
#
#    -  How to save a mesh to a file.
#    -  How to load a mesh from a file.
#    -  How to find out the number of elements and their type.
#    -  How to find out the characteristics of the shape (shape descriptor).

# We will generate the triangular mesh inside a rectangular block.
# The block will have the dimensions shown below:
a, b = 1.0, 1.0
# The triangles will be generated in a regular pattern, with the number of
# edges per side of the block given as
na, nb = 2, 3

# The mesh will be generated by the package `MeshSteward`. We choose an
# orientation of the diagonals (in other words the pattern in which the block is
# tiled) as labeled `:b`.
using MeshSteward: T3block
conn = T3block(a, b, na, nb, :b);

# The variable `conn` is an incidence relation. This incidence relation will
# become the base relation of the mesh when it gets attached to the mesh.
using MeshSteward: Mesh, attach!
m = Mesh()
attach!(m, conn);

# The embedding space in which the mesh lives is two-dimensional. We can verify
# that by
using MeshSteward: nspacedims
@show nspacedims(m)

# We can export the mesh into a set of files that are referred to as the MESH
# format. This will actually store the mesh as three separate files. The file
# `"Unit-square-mesh.mesh"` is the directory of the bundle of the mesh files,
# so to speak.
using MeshSteward: save
save(m, "Unit-square-mesh")
# Note the trio of files `"Unit-square-mesh.mesh"` (mesh  "directory"), and
# `"Unit-square-mesh-conn.dat"` (connectivity) and
# `"Unit-square-mesh-xyz.dat"` (locations of vertices).
readdir(".")

# We will now load the mesh we saved just now into another mesh that we create
# specifically for the purpose of comparing the two meshes.
using MeshSteward: load
m2 = load(Mesh(), "Unit-square-mesh");

# Now we compare the number of  elements  stored in those two meshes. First the
# vertices (those are the shapes on the right of the incidence relation).
using MeshCore: nshapes
using MeshSteward: baseincrel
@show nshapes(baseincrel(m).right) == nshapes(baseincrel(m2).right)

# Now  we compare the number of the triangular elements in those two meshes:
@show nshapes(baseincrel(m).left) == nshapes(baseincrel(m2).left)

# What is the name of the shape descriptor of the shape collection on the left
# (that is what we usually consider the "elements" in a finite element mesh) of
# the base incidence relation of the mesh?
using MeshCore: shapedesc
@show shapedesc(baseincrel(m).left).name

# We would expect it to be the name of a three-node triangle, "T3".

# Are the shape descriptors of the shape collection on the left the same for the
# two meshes?
@show shapedesc(baseincrel(m).left) == shapedesc(baseincrel(m2).left)

# What is the name of the shape descriptor of the shape collection on the right
# (that is what we usually consider the "nodes" in a finite element mesh) of
# the base incidence relation of the mesh?
using MeshCore: shapedesc
@show shapedesc(baseincrel(m).right).name

# P1 is the shape descriptor of a point shape. So vertices (nodes) are points.

# Are the shape descriptors of the shape collection on the right the same for the
# two meshes?
@show shapedesc(baseincrel(m).right) == shapedesc(baseincrel(m2).right)


text:
@testset "Functions" begin
    @polyvar x
    α = MOI.VariableIndex(1)
    @test α * x isa AbstractTerm{MOI.VariableIndex}
    @test x * α isa AbstractTerm{MOI.VariableIndex}
    @test α * x^2 isa AbstractTerm{MOI.VariableIndex}
    @test x^2 * α isa AbstractTerm{MOI.VariableIndex}
    @test (1α) * x^2 isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test (1x)^2 * α isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test α * (1x) isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test (1x) * α isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test α * (1.0x) isa AbstractTerm{MOI.ScalarAffineFunction{Float64}}
    @test (1.0x) * α isa AbstractTerm{MOI.ScalarAffineFunction{Float64}}
    @test α * (1x + 2) isa AbstractPolynomial{MOI.ScalarAffineFunction{Int}}
    @test (1x + 2) * α isa AbstractPolynomial{MOI.ScalarAffineFunction{Int}}
    @test α * (1.0x + 2.0) isa AbstractPolynomial{MOI.ScalarAffineFunction{Float64}}
    @test (1.0x + 2.0) * α isa AbstractPolynomial{MOI.ScalarAffineFunction{Float64}}
end


text:
using GLVisualize, GeometryTypes, GLAbstraction
using Colors, Reactive, FileIO
if !isdefined(:runtests)
	window = glscreen()
	timesignal = bounce(linspace(0f0, 1f0, 360))
end

cat    = GLNormalMesh(loadasset("cat.obj"))
sphere = GLNormalMesh(Sphere{Float32}(Vec3f0(0), 1f0), 12)

function scale_gen(v0, nv)
	l = length(v0)
	@inbounds for i=1:l
		v0[i] = Vec3f0(1,1,sin((nv*l)/i))/2
	end
	v0
end
function color_gen(v0, t)
	l = length(v0)
	@inbounds for x=1:l
		v0[x] = RGBA{U8}(x/l,(cos(t)+1)/2,(sin(x/l/3)+1)/2.,1.)
	end
	v0
end

t      		 = const_lift(x->x+0.1, timesignal)
ps 			 = sphere.vertices
scale_start  = Vec3f0[Vec3f0(1,1,rand()) for i=1:length(ps)]
scale 		 = foldp(scale_gen, scale_start, t)
colorstart 	 = color_gen(zeros(RGBA{U8}, length(ps)), value(t))
color  		 = foldp(color_gen, colorstart, t)
rotation 	 = -sphere.normals

cats = visualize((cat, ps), scale=scale, color=color, rotation=rotation)

_view(cats, window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
module TestCUDA
import CUDA
using Test
using ONNXRunTime
const OX = ONNXRunTime
using ONNXRunTime: SessionOptionsAppendExecutionProvider_CUDA

#using Libdl
#Libdl.dlopen("/home/jan/.julia/artifacts/e2fd6cdf04b830a1d802fb35a6193788d0a3811a/lib/libcudart.so.11.0")

@testset "CUDA high level" begin
    @testset "increment2x3.onnx" begin
        path = OX.testdatapath("increment2x3.onnx")
        model = OX.load_inference(path, execution_provider=:cuda)
        @test OX.input_names(model) == ["input"]
        @test OX.output_names(model) == ["output"]
        input = randn(Float32, 2,3)
        y = model((;input=input,), ["output"])
        @test y == (output=input .+ 1f0,)
    end
end

using ONNXRunTime.CAPI
@testset "CUDA low level" begin
    api = GetApi(execution_provider=:cuda)
    env = CreateEnv(api, name="myenv")
    path = ONNXRunTime.testdatapath("increment2x3.onnx")
    session_options = CreateSessionOptions(api)
    cuda_options = OrtCUDAProviderOptions()
    SessionOptionsAppendExecutionProvider_CUDA(api, session_options, cuda_options)
    session = CreateSession(api, env, path, session_options)
    @test SessionGetInputCount(api, session) == 1
    @test SessionGetOutputCount(api, session) == 1
    mem = CreateCpuMemoryInfo(api)
    allocator = CreateAllocator(api, session, mem)
    @test SessionGetInputName(api, session, 0, allocator) == "input"
    @test_throws OX.OrtException SessionGetInputName(api, session, 1, allocator)
    @test SessionGetOutputName(api, session, 0, allocator) == "output"
    @test_throws OX.OrtException SessionGetOutputName(api, session, 1, allocator)
    input_vec = randn(Float32, 6)
    input_array = [
        input_vec[1] input_vec[2] input_vec[3];
        input_vec[4] input_vec[5] input_vec[6];
    ]
    input_tensor = CreateTensorWithDataAsOrtValue(api, mem, input_vec, (2,3))
    run_options = CreateRunOptions(api)
    input_names = ["input"]
    output_names = ["output"]
    inputs = [input_tensor]
    outputs = Run(api, session, run_options, input_names, inputs, output_names)
    @test length(outputs) == 1
    output_tensor = first(outputs)
    @test output_tensor isa OrtValue
    output_array = GetTensorMutableData(api, output_tensor)
    @test typeof(output_array) <: AbstractMatrix{Float32}
    @test size(output_array) == (2,3)
    @test output_array ≈ 1 .+ input_array
end
end#module


text:
function get_affine_map(ds::ArchGDAL.Dataset)
    # ArchGDAL fails hard on datasets without
    # an affinemap. GDAL documents that on fail
    # a default affinemap should be returned.
    try
        global gt = ArchGDAL.getgeotransform(ds)
    catch y
        @warn y.msg
        global gt = [0.,1.,0.,0.,0.,1.]
    end
    geotransform_to_affine(SVector{6, Float64}(gt))
end

function geotransform_to_affine(gt::SVector{6, Float64})
    # See https://lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html
    # for an explanation of the geotransform format
    AffineMap(SMatrix{2,2}([gt[2] gt[3]; gt[5] gt[6]]), SVector{2}([gt[1], gt[4]]))
end
geotransform_to_affine(A::Vector{Float64}) = geotransform_to_affine(SVector{6}(A))

function affine_to_geotransform(am::AffineMap{SArray{Tuple{2,2},Float64,2,4},SArray{Tuple{2},Float64,1,2}})
    l = am.linear
    t = am.translation
    [t[1], l[1], l[2], t[2], l[3], l[4]]
end

function affine_to_geotransform(am::AffineMap{Array{Float64,2},Array{Float64,1}})
    l = am.linear
    t = am.translation
    (length(l) != 4 || length(t) != 2) || error("AffineMap has wrong dimensions.")
    [t[1], l[1], l[2], t[2], l[3], l[4]]
end

"""Check wether the AffineMap of a GeoArray contains rotations."""
function is_rotated(ga::GeoArray)
    ga.f.linear[2] != 0. || ga.f.linear[3] != 0.
end

function bbox(ga::GeoArray)
    ax, ay = ga.f(SVector(0,0))
    bx, by = ga.f(SVector(size(ga)[1:2]))
    (min_x=min(ax, bx), min_y=min(ay, by), max_x=max(ax, bx), max_y=max(ay, by))
end

# Placeholder for setting bbox
function bbox!(ga::GeoArray, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y),Tuple{Float64, Float64, Float64, Float64}})
    nothing
end

"Generate bounding boxes for GeoArray cells."
function bboxes(ga::GeoArray)
    c = coords(ga)::Array{StaticArrays.SArray{Tuple{2},Float64,1,2},2}
    m, n = size(c)
    cellbounds = Matrix{NamedTuple}(undef, (m-1, n-1))
    for j = 1:n-1, i = 1:m-1
        v = c[i:i+1, j:j+1]
        minx, maxx = extrema(first.(v))::Tuple{Float64, Float64}
        miny, maxy = extrema(last.(v))::Tuple{Float64, Float64}
        cellbounds[i, j] = (min_x=minx, max_x=maxx, min_y=miny, max_y=maxy)
    end
    cellbounds
end

# Extend CoordinateTransformations
CoordinateTransformations.compose(ga::GeoArray, t2::AffineMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::LinearMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::Translation) = CoordinateTransformations.compose(ga.f, t2)

"""Transform an GeoArray by applying a Transformation."""
function compose!(ga::GeoArray, t2::AffineMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::LinearMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::Translation)
    ga.f = compose(ga, t2)
end

"""Function to flip GeoArray upside down to adjust to GDAL ecosystem."""
function flipud!(ga::GeoArray)
    # Flip data upside down
    ga.A = reverse(ga.A, dims=2)

    # Find new corner coordinates
    _, uy = ga.f(SVector{2}([0, size(ga)[2]]))
    _, ux = ga.f(SVector{2}([size(ga)[1], 0]))

    # Define y mirror and compose
    lm = LinearMap(SMatrix{2,2}([1.0 0.0; 0.0 -1.0]))
    am = compose(lm, ga.f)  # AffineMap
    translate = SVector{2}([ux, uy])
    f = AffineMap(am.linear, translate)

    ga.f = f
    ga
end