text:
@testset "light" begin
    M.set_interface_mode(M.LightInterface())
    @test M.get_interface_mode() isa M.LightInterface

    # For matrix objects (:other) we don't need `FullInterface`
    # to run `matrix` method as they are already matrices.
    X   = zeros(3, 4)
    mX  = matrix(X)
    mtX = matrix(X; transpose=true)

    @test mX === X
    @test mtX == permutedims(X)

    # for other objects `:table` or `:other` we need
    # `FullInterface` to get the corresponding trait
    X = (x=[1, 2, 3], y=[1, 2, 3])
    @test_throws M.InterfaceError M.vtrait(X)
    @test_throws M.InterfaceError matrix(X)
    X = (1, 2, 3, 4)
    @test_throws M.InterfaceError M.vtrait(X)
    @test_throws M.InterfaceError matrix(X)
end

@testset "full" begin
    M.set_interface_mode(M.FullInterface())
    @test M.get_interface_mode() isa M.FullInterface

    M.matrix(::M.FullInterface, ::Val{:table}, X; kw...) = Tables.matrix(X; kw...)
    M.vtrait(::FI, X, s) = Val{ifelse(Tables.istable(X), :table, :other)}()

    X = (x=[1, 2, 3], y=[1, 2, 3])
    mX = matrix(X)
    @test mX isa Matrix
    @test mX == hcat(X.x, X.y)
end


text:
#Fermats Theorem
#https://en.wikipedia.org/wiki/Fermat%27s_little_theorem

Fermat's little theorem is the basis for the Fermat primality test and is one of the fundamental results of elementary number theory. The theorem is named after Pierre de Fermat, who stated it in 1640. It is called the "little theorem" to distinguish it from Fermat's Last Theorem.[3]

function checkfermat(a::Int64,b::Int64,c::Int64,n::Int64)
    if (a^n) + (b^n) == (c^n)
        return (0)
    elseif (a^n)+(b^n) != (c^n)
        return (1)
    end
end

a=Array(1:50)
function check_(x::Array{Int64,1})
    o=0
    for a in x
        for b in x
            for c in x
                for n in x
                    q=checkfermat(a,b,c,n)
                    o=o+1
                    if q==0
                        if n>2
                            println("Theorem Proved WRONG")
                            println(a,'@',b,'@',c,'@',n)
                        end
                    end
                end
            end
        end
    end
    println(o)
end

check_(a)

ext: jl


text:
export
    NextEvent,
    LinearBound,
    nextevent_bps,
    nextevent_bps_q,
    nextevent_zz

abstract type IPPSamplingMethod end
abstract type Thinning <: IPPSamplingMethod end

"""
    LinearBound

Indicating you have access to a global bound on the eigenvalues of the Hessian.
"""
struct LinearBound <: Thinning
    xstar::Vector{Real}
    gllstar::Vector{Real}  # grad log likelihood around xstar
    b::Real                # N*L where L is Lipschitz constant
    a::Function
    function LinearBound(xstar, gllstar, b)
        new(xstar, gllstar, b,
            (x, v) -> max(dot(-gllstar, v), 0.) + norm(x-xstar) * b)
    end
end

"""
    NextEvent

Object returned when calling a `nextevent_*` type function. The bouncing time
`tau` as well as whether it should be accepted or rejected `dobounce` (set to
true always for analytical sampling). `flipindex` is to support ZZ sampling.
"""
struct NextEvent
    tau::Real # candidate first arrival time
    dobounce::Function # do bounce (if accept reject step)
    flipindex::Int     # flipindex (if ZZ style)
    function NextEvent(tau, dobounce, flipindex)
        new(tau, dobounce, flipindex)
    end
end
NextEvent(tau; dobounce=(g,v)->true, flipindex=-1) =
    NextEvent(tau, dobounce, flipindex)

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

"""
    nextevent_bps(g::MvGaussian, x, v)

Return a bouncing time and corresponding intensity in the specific case of a
Gaussian density (for which the simulation of the first arrival time of the
corresponding IPP can be done exactly).
"""
function nextevent_bps(g::MvGaussian,
                       x::Vector{<:Real}, v::Vector{<:Real})
    # precision(g)*x - precision(g)*mu(g) --> precmult, precmu in Gaussian.jl
    a = dot(mvg_precmult(g, x) - mvg_precmu(g), v)
    b = dot(mvg_precmult(g, v), v)
    c = a / b
    e = max(0.0, c) * c # so either 0 or c^2
    tau = -c + sqrt(e + 2randexp() / b)
    return NextEvent(tau)
end

function nextevent_bps(g::PMFGaussian,
                       x::Vector{<:Real}, w::Vector{<:Real})
    # unmasking x, w
    xu, xv = x[1:g.d], x[g.d+1:end]
    wu, wv = w[1:g.d], w[g.d+1:end]
    # precomputing useful dot products
    xuwv, xvwu = dot(xu, wv), dot(xv, wu)
    wuwv, dxw = dot(wu, wv), xuwv + xvwu
    # real root
    t0 = -0.5dxw / wuwv
    # e(x) := ⟨x_u, x_v⟩ - r
    ex = dot(xu, xv) - g.r
    # (quadratic in t) e(x+tw) = e(x)+(⟨wu,xv⟩+⟨wv,xu⟩)t+⟨wu,wv⟩t^2
    p1 = Poly([ex, dxw, wuwv])
    # (linear in t) ⟨xv,wu⟩+⟨xu,wv⟩+2⟨wu,wv⟩t
    p2 = Poly([dxw, 2.0wuwv])
    # ⟨∇E(x+tw),w⟩ is a cubic in t (and the intensity)
    p  = p1 * p2

    rexp = randexp()
    tau = 0.0

    ### CASES (cf. models/pmf.jl)

    Δ = t0^2 - ex / wuwv # discriminant of p1
    if Δ <= 0
        # only single real root (t0)
        # two imaginary roots
        if t0 <= 0
            tau = pmf_caseA(rexp, p)
        else
            tau = pmf_caseB(rexp, p, t0)
        end
    else
        # three distinct real roots,
        tm, tp = t0 .+ sqrt(Δ) * [-1.0,1.0]
        if tp <= 0 # case: |/
            tau = pmf_caseA(rexp, p)
        elseif t0 <= 0 # case: |_/
            tau = pmf_caseB(rexp, p, tp)
        elseif tm <= 0 # case: |\_/
            tau = pmf_caseC(rexp, p, t0, tp)
        else # case: |_/\_/
            tau = pmf_caseD(rexp, p, tm, t0, tp)
        end
    end
    return NextEvent(tau)
end

"""
    nextevent_zz(g::MvGaussian, x, v)

Same as nextevent but for the Zig Zag sampler.
"""
function nextevent_zz(g::MvGaussian,
                      x::Vector{<:Real}, v::Vector{<:Real})
    # # precision(g)*x - precision(g)*mu(g) --> precmult, precmu in Gaussian.jl
    u1 = mvg_precmult(g, x)-mvg_precmu(g)
    u2 = mvg_precmult(g, v)
    taus = zeros(g.p)
    for i in 1:g.p
        ai = u1[i] * v[i]
        bi = u2[i] * v[i]
        ci = ai ./ bi
        ei = max(0.0, ci) * ci
        taus[i] = -ci + sqrt(ei + 2.0randexp() / abs(bi))
    end
    tau, flipindex = findmin(taus)
    return NextEvent(tau, flipindex=flipindex)
end

"""
    nextevent_bps(lb::LinearBound, x, v)

Return a bouncing time and corresponding intensity corresponding to a linear
upperbound described in `lb`.
"""
function nextevent_bps(lb::LinearBound,
                       x::Vector{<:Real}, v::Vector{<:Real})
    a = lb.a(x, v)
    b = lb.b
    @assert a >= 0.0 && b > 0.0 "<ippsampler/nextevent_bps/linearbound>"
    tau = -a / b + sqrt((a / b)^2 + 2randexp() / b)
    lambdabar = a + b * tau
    return NextEvent(tau, dobounce=(g,v)->(rand()<-dot(g, v)/lambdabar))
end

function nextevent_bps_q(gll::Function, x::Vector{<:Real}, v::Vector{<:Real},
                         tref::Real; n=100)

    chi(t) = max(0.0, dot(-gll(x + t * v), v))
    S      = Chebyshev(0.0..tref)
    p      = points(S, n)
    v      = chi.(p)
    f      = Fun(S, ApproxFun.transform(S,v))
    If     = cumsum(f) # integral from 0 to t with t < = tref
    tau    = ApproxFun.roots(If - randexp())
    tau    = (length(tau)>0) ? minimum(tau) : Inf
    return NextEvent(tau)
end


text:
# Testing:
# @assert 1==42 "Shouldn't this be so?"
#> ERROR: assertion failed: Shouldn't this be so?
using Base.Test
# @test 1==3 #> ERROR: test failed: 1 == 3
# @test_approx_eq 1 1.1
#> ERROR: assertion failed: |1 - 1.1| <= 2.220446049250313e-12
@test_approx_eq_eps 1 1.1 0.2

# Debugging:
arr = [1, 2] #> 2-element Array{Int64,1}: 1 2
@which sort(arr) #> sort(v::AbstractArray{T,1}) at sort.jl:334
456 * 789 + (@show 2 + 3) #>  2 + 3 => 5     359789

# Benchmarking:
@time [x^2 for x in 1:1000] #> elapsed time: 3.911e-6 seconds (8064 bytes allocated) # 1000-element Array{Int64,1}: …
@timed [x^2 for x in 1:1000]
#> ([1,4,9,16,25,36,49,64,81,100  …  982081,984064,986049,988036,990025,992016,9940
# 09,996004,998001,1000000],3.911e-6,8064,0.0)
@elapsed [x^2 for x in 1:1000] #> 3.422e-6
@allocated [x^2 for x in 1:1000] #> 8064
# timing:
tic() # start timer
[x^2 for x in 1:1000]
toc() # stop timer and print time: elapsed time: 0.024395069 seconds
# toq() # stop timer and return time

# starting a task:
a = @async 1 + 2
# Task (done) @0x000000002d70faf0
consume(a) # 3


text:
"""

```
writePDB(atoms::Vector{Atom}, filename, selection)
```

Function that writes a PDB file given the vector of atoms, with
optional definition of a selection to be print.

"""
function writePDB(atoms::AbstractVector{Atom}, filename, selection)
  query = parse_query(selection)
  writePDB(atoms,filename,only=atom->apply_query(query,atom))
end

function writePDB(atoms::AbstractVector{Atom}, filename; only = all)
  file = open(filename,"w")
  for atom in atoms
    if only(atom)
      println(file,write_atom(atom))
    end
  end
  close(file)
end


text:
function foo()
  return 123
end


text:
module DiagnosisClassification

import CSV
import Downloads
import Pkg
import ProgressMeter
import Scratch
import TOML

include("types.jl")

include("dict_utils.jl")
include("download.jl")

include("convenience/icd9.jl")
include("convenience/icd10.jl")

include("dots/icd9.jl")
include("dots/icd10.jl")

include("init.jl")

include("TODO.jl")

end # module


text:
using ITensors

import ITensors: op

function op(::OpName"expτSS", ::SiteType"S=1/2", s1::Index, s2::Index; τ)
  h = 1/2 * op("S+",s1) * op("S-",s2) +
      1/2 * op("S-",s1) * op("S+",s2) +
            op("Sz",s1) * op("Sz",s2)
  return exp(τ * h)
end

function main(; N = 10, cutoff = 1E-8,
                δt = 0.1, ttotal = 5.0)
  # Compute the number of steps to do
  Nsteps = Int(ttotal / δt)

  # Make an array of 'site' indices
  s = siteinds("S=1/2", N; conserve_qns = true)

  # Make gates (1,2),(2,3),(3,4),...
  gates = ops([("expτSS", (n, n+1), (τ = -δt*im/2,)) for n in 1:N-1], s)

  # Include gates in reverse order too
  # (N,N-1),(N-1,N-2),...
  append!(gates, reverse(gates))

  # Function that measures <Sz> on site n
  function measure_Sz(psi::MPS, n)
    psi = orthogonalize(psi, n)
    sn = siteind(psi,n)
    Sz = scalar(dag(prime(psi[n],"Site"))*op("Sz",sn)*psi[n])
    return real(Sz)
  end

  # Initialize psi to be a product state (alternating up and down)
  psi0 = productMPS(s, n -> isodd(n) ? "Up" : "Dn")

  c = div(N,2)

  # Compute and print initial <Sz> value
  t = 0.0
  Sz = measure_Sz(psi0,c)
  println("$t $Sz")

  # Do the time evolution by applying the gates
  # for Nsteps steps
  psi = psi0
  for step in 1:Nsteps
    psi = apply(gates, psi; cutoff=cutoff)
    t += δt
    Sz = measure_Sz(psi,c)
    println("$t $Sz")
  end

  # Now do the same evolution with an MPO
  rho0 = MPO(psi0)
  rho = rho0
  for step in 1:Nsteps
    rho = apply(gates, rho; cutoff=cutoff, apply_dag = true)
    t += δt
  end
  @show inner(psi, rho, psi)
  @show inner(psi, psi)
  @show tr(rho)

  return nothing
end


text:
using ClenshawCurtisBessel
using Test
using DoubleFloats
using LinearAlgebra

# """
# Consider an example recurrence relation, which we will call `Example61161`
# ```math
# 6 y(s) - 11 y(s+1) + 6 y(s+2) - y(s+3) = 0.
# ```
# with initial conditions ``y(1) = 2``, ``y(2)=5``, and ``y(3) = 15``.
# """
@testset "Oliver's Method" begin
    y61161(s) = 1 - 2^(s-1) + 2 * 3^(s-1)
    struct Example61161{T,N,M} <: ClenshawCurtisBessel.OliverProblem{T,N,M} end
    function ClenshawCurtisBessel.OliverP(OP::Example61161, s, i)
        if s == 0
            return -1
        elseif s == 1
            return 6
        elseif s == 2
            return -11
        elseif s == 3
            return 6
        end
    end
    ClenshawCurtisBessel.OliverR(OP::Example61161, i) = 0

    # loop over a broad range of BC and sizes
    for a in -2:10
        for b in (a+3):(a+400)
            for m in 0:3
                a, b, n = 1, 9, 3

                YBC = Dict{Int,Int}()
                for i in (a+1):(a+n-m)
                    YBC[i] = y61161(i)
                end
                for i in (b-m+1):b
                    YBC[i] = y61161(i)
                end

                ex = Example61161{Float64,n,m}()
                P = ClenshawCurtisBessel.assembleP(ex, a, b)
                ρ = ClenshawCurtisBessel.assembleρ(ex, a, b, YBC)

                # """$(P * y.((a+n-m+1):(b-m))) $(ρ)"""
                sol = P \ ρ
                ref = y61161.((a+n-m+1):(b-m))
                @test all( ((sol .- ref) ./ ref) .< 1e-11 )
            end
        end
    end
end

# MATHEMATICA INPUT FOR THESE INTEGRALS
# NumberForm[NIntegrate[1/2 ChebyshevT[k, x] BesselJ[\[Nu], (x + 1) a/2] /.
# {a -> 10, \[Nu] -> 3 , k -> 3 + 11} , {x, -1, 1}, AccuracyGoal -> 40,
# PrecisionGoal -> 40, WorkingPrecision -> 40], 40]
@testset "Moment Boundary Conditions" begin
    @test ClenshawCurtisBessel.bessel_moment_G(10, 15/2, 2) ≈ 0.007089503058490749
    @test ClenshawCurtisBessel.bessel_moment_G(100, 15/2, 2) ≈ 1.482273470904133e24
    @test ClenshawCurtisBessel.bessel_moment_G(70, 59/2, 1) ≈ 2.695679919047616e-43
    @test ClenshawCurtisBessel.bessel_moment_G(100, 5000, 1) ≈ 0.0
    @test ClenshawCurtisBessel.bessel_moment_G(1, 5, 1) ≈ 0.00003601420867202275

    @test ClenshawCurtisBessel.Mk_asymptotic(2.0, 3.0, 21.0) ≈ -0.0001477774680739161 rtol=1e-6
    @test ClenshawCurtisBessel.Mk_asymptotic(3.0, 5.0, 21.0) ≈ -0.00004961332181983246 rtol=1e-6
    @test ClenshawCurtisBessel.Mk_asymptotic(2.0, 6.0, 24.0) ≈ -1.06156879586713e-6 rtol=1e-6
    @test ClenshawCurtisBessel.Mk_asymptotic(3., 7.0, 28.0) ≈ -1.647048323970303e-6 rtol=1e-6
end

@testset "Double64 Raw Moments" begin
	CCB = ClenshawCurtisBessel
    bpp = CCB.BrandersPiessensProblem{Double64,8,2}(10.0, 3 + 1/2)
	indexBC, MBC = CCB.generate_BC(bpp, 256)
	ρ = CCB.assembleρ(bpp, -3, indexBC, MBC)
	𝐏 = CCB.assembleP(bpp, -3, indexBC)
	sol = qr(𝐏) \ ρ
	ref = big"0.0002069511037367724863632484164263304887654"
	@test (sol[11] .- ref) ./ ref < 1e-15

    bpp = CCB.BrandersPiessensProblem{Double64,8,2}(10.0, 3)
	indexBC, MBC = CCB.generate_BC(bpp, 256)
	ρ = CCB.assembleρ(bpp, -3, indexBC, MBC)
	𝐏 = CCB.assembleP(bpp, -3, indexBC)
	sol = qr(𝐏) \ ρ
	ref = big"-0.0001955581634668262545942178152189286110836"
	@test (sol[11] .- ref) ./ ref < 1e-15
	ref = big"-0.00009457154862533141412967956613252113399756"
	@test (sol[16] .- ref) ./ ref < 1e-15
	ref = big"-0.00001631452917562497345315040327964024786504"
	@test (sol[40] .- ref) ./ ref < 1e-15
end

@testset "Float64 Raw Moments" begin
	CCB = ClenshawCurtisBessel
    bpp = CCB.BrandersPiessensProblem{Float64,8,2}(10.0, 3 + 1/2)
	indexBC, MBC = CCB.generate_BC(bpp, 256)
	ρ = CCB.assembleρ(bpp, -3, indexBC, MBC)
	𝐏 = CCB.assembleP(bpp, -3, indexBC)
	sol = qr(𝐏) \ ρ
	ref = big"0.0002069511037367724863632484164263304887654"
	@test (sol[11] .- ref) ./ ref < 1e-6

    bpp = CCB.BrandersPiessensProblem{Float64,8,2}(10.0, 3)
	indexBC, MBC = CCB.generate_BC(bpp, 256)
	ρ = CCB.assembleρ(bpp, -3, indexBC, MBC)
	𝐏 = CCB.assembleP(bpp, -3, indexBC)
	sol = qr(𝐏) \ ρ
	ref = big"-0.0001955581634668262545942178152189286110836"
	@test (sol[11] .- ref) ./ ref < 1e-6
end

@testset "Double64 momentM" begin

    ref = [big"0.0007290114582988161923809468346226923322190",
        big"0.0005938052788264036365639278651356395999367",
        big"0.0002786292946361793775883850681999983130955",
        big"-0.00002805260432115378687369846500283467077757",
        big"-0.0001950344048967695886304586424720522836879",
        big"-0.0002181086699578411818278638336785129268495",
        big"-0.0001688344453897982834345306352230680868134",
        big"-0.0001124694826392090914811134620250485009294",
        big"-0.00007437704414078541354435912139888596042476",
        big"-0.00005298335002672124481290725819716655026908",
        big"-0.00004066037955955940010097890154609528857357",
        big"-0.00003259729530260456078582391639304253599224",
        big"-0.00002680486467601182798029350986861852376451",
        big"-0.00002246247645744408637896618989749724244207",
        big"-0.00001911995263563955270424343152423340974346",
        big"-0.00001648755827827632619531486527620890550541",
        big"-0.00001437360220319041810429287944102381603564",
        big"-0.00001264808824960577566084971381515466084373",
        big"-0.00001121999433319545122043778233304758581356",
        big"-0.00001002382210812218822542518716028642086529",
        big"-9.011366672874117063012554033327179984774e-6"]
    # note that res starts at 0
    res = momentM(Float64, DoubleFloat(6.0), DoubleFloat(10.0), 20)
    for i in 0:19
        @test abs(res[i] .- ref[i+1]) < 1e-16
    end

    # now test the half-integer case
    ref = [big"0.0003784769456217184219402128017246464332000",
        big"0.0003116266706253777562730807652926920232194",
        big"0.0001539041845785743007034880912831886326126",
        big"-3.947735731437409097894699973916750717989e-6",
        big"-0.00009548112619274286495691961673995490768193",
        big"-0.0001140573571560400900107740517746052585751",
        big"-0.00009226229577888985713288717817974834079267",
        big"-0.00006329298661640707880778541216716573215075",
        big"-0.00004220681212615568225179049644556672266289",
        big"-0.00002983588500153932889283750688859037769246",
        big"-0.00002268346477382750545447550291183506504071"]
    # note that res starts at 0
    res = momentM(Float64, DoubleFloat(6.0), DoubleFloat(10.5), 20)
    for i in 0:10
        @test abs(res[i] .- ref[i+1]) < 1e-16
    end
end


text:
function FVort2VecDSS(v1,v2,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
dXdxIC = Param.cache.dXdxIC
v1CG=reshape(v1[reshape(CG.Glob,OP*OP*NF,1),:]
  ,OP,OP,NF,nz);
v2CG=reshape(v2[reshape(CG.Glob,OP*OP*NF,1),:]
  ,OP,OP,NF,nz);
vC1=reshape(
  CG.DS*reshape(dXdxIC[:,:,:,:,1,1].*v2CG -
  dXdxIC[:,:,:,:,1,2].*v1CG
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz) -
  permute(
  reshape(
  CG.DS*reshape(
  permute(-v2CG.*dXdxIC[:,:,:,:,2,1] +
  v1CG.*dXdxIC[:,:,:,:,2,2]
  ,[2 1 3 4])
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz)
  ,[2 1 3 4]);

Vort=zeros(CG.NumG,nz);
for iM=1:size(CG.FaceGlob,1)
  i = CG.FaceGlob[iM].Ind
  arr = reshape(CG.Glob[:,i,:],OP*OP*size(i,1))
  Vort[arr,:] = Vort[arr,:] .+
    reshape(vC1[:,:,i,:]
    ,OP*OP*size(i,1),nz);
end
Vort=Vort./CG.M;
return Vort
end


text:
using Test
using DiffEqBase
using Sundials

@testset "Test Reactors" begin
@testset "Test liquid phase reactor simulation" begin
#Constant T and V Ideal Dilute Liquid
phaseDict = readinput("../src/testing/liquid_phase.rms")
spcs = phaseDict["phase"]["Species"]; #mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
rxns = phaseDict["phase"]["Reactions"];
solv = phaseDict["Solvents"][1];
liq = IdealDiluteSolution(spcs,rxns,solv;name="phase",diffusionlimited=true) #Define the phase (how species thermodynamic and kinetic properties calculated)
initialconds = Dict(["T"=>450.0,"V"=>1.0e-6*1e6,"octane"=>6.154e-3*1e6,"oxygen"=>4.953e-6*1e6]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantTVDomain(phase=liq,initialconds=initialconds,constantspecies=["oxygen"]) #Define the domain (encodes how system thermodynamic properties calculated)
react = Reactor(domain,y0,(0.0,140000.01);p=p) #Create the reactor object

sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-8); #solve the ode associated with the reactor

spcnames = getfield.(liq.species,:name)
octaneind = findfirst(isequal("octane"),spcnames)
y = sol(32977.61568)
@test y[octaneind]/sum(y) ≈ 0.461599061 rtol=3e-2 #from RMG simulator I believe the slight difference is due to better calculation of diffusion limits in RMS
end;

#Use superminimal example to test
phaseDict = readinput("../src/testing/superminimal.rms") #load mechanism dictionary
spcs = phaseDict["phase"]["Species"]; #mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
rxns = phaseDict["phase"]["Reactions"];
ig = IdealGas(spcs,rxns,name="phase")

#Constant T and P Ideal Gas
@testset "Test constant T and P reactor simulation" begin
 #Define the phase (how species thermodynamic and kinetic properties calculated)
initialconds = Dict(["T"=>1000.0,"P"=>1e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantTPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,150.11094);p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor
sim = Simulation(sol,domain);

spcnames = getfield.(ig.species,:name)
h2ind = findfirst(isequal("H2"),spcnames)
o2ind = findfirst(isequal("O2"),spcnames)
h2oind = findfirst(isequal("H2O"),spcnames)
y = sol(20.44002454)
N = sim.N(20.44002454)
@test y[h2ind]/N ≈ 0.412883111 rtol=1e-4 #from RMG simulator
@test y[o2ind]/N ≈ 0.200419093 rtol=1e-4
@test y[h2oind]/N ≈ 0.386618602 rtol=1e-4

#sensitivities
dps = getadjointsensitivities(sim,"H2",CVODE_BDF();sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),abstol=1e-16,reltol=1e-6)
react2 = Reactor(domain,y0,(0.0,150.11094);p=p,forwardsensitivities=true)
sol2 = solve(react2.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-6); #solve the ode associated with the reactor
sim2 = Simulation(sol2,domain)

x,dp = extract_local_sensitivities(sol2,150.11094);
ind = findfirst(isequal("H2"),sim2.names)
dpvs = [v[ind] for v in dp]
dpvs[length(domain.phase.species)+1:end] .*= domain.p[length(domain.phase.species)+1:end]
dpvs ./= sol2(150.11094)[ind]
rerr = (dpvs .- dps')./dpvs
rerr = [isinf(x) ? 0.0 : x for x in rerr]
@test all((abs.(rerr) .> 1e-2).==false)
end;

#Constant V adiabatic Ideal Gas
#uses superminimal.yml mechanism
@testset "Constant volume adiabatic reactor simulation" begin
initialconds = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantVDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,0.101),p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor

ts = exp.(range(log(1e-15),length=10000,stop=log(0.1)))
IDT = ts[argmax(diff([sol(t)[end] for t in ts]))] #Ignition Delay Time based on argmax(dTdt(t))

@test IDT ≈ 0.038384723436228063 rtol=1e-5 #from Cantera simulation
end;

#Constant P adiabatic Ideal Gas
#uses ethane.rms mechanism
@testset "Constant pressure adiabatic reactor simulation" begin

phaseDict = readinput("../src/testing/ethane.rms")
spcs = phaseDict["phase"]["Species"]
rxns = phaseDict["phase"]["Reactions"]
ig = IdealGas(spcs,rxns,name="phase")

initialconds = Dict(["T"=>1000.0,"P"=>2.0e5,"ethane"=>1.0,"Ar"=>1.0,"O2"=>3.5]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,1.0);p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6); #solve the ode associated with the reactor

ts = exp.(range(log(1e-15),length=10000,stop=log(0.2)))
IDT = ts[argmax(diff([sol(t)[end] for t in ts]))] #Ignition Delay Time based on argmax(dTdt(t))

@test IDT ≈ 0.07324954954380769 rtol=1e-5
end;

end;


text:
module VehicleCalculate
include("VehicleStruct.jl")
include("VehicleFunction.jl")

export vehiclespeed,tableplot
export PrimTransmission,DrivingTransmission,DrivenTransmission,Flywheel,Transmission
# Write your package code here.

end

max_stars_repo_name: powerzf/VehicleCalculate

max_issues_repo_name: powerzf/VehicleCalculate

max_forks_repo_name: powerzf/VehicleCalculate


text:
using Test
using FiniteDifferences
using Random: seed!
using Zygote, Flux, ForwardDiff, LinearAlgebra
using AugmentedGaussianProcesses; const AGP = AugmentedGaussianProcesses
@info "Packages loaded"
seed!(42)
M = 10
N = 100
D= 2
T = 3000
fullgp = false
X = rand(N,D)*2 .-1.0
y = norm.(eachrow(X)).-0.5#randn(N)
xgrid = range(-1,1,length=100)
Xtest = vcat(transpose.(collect.(Iterators.product(xgrid,xgrid)))...)
kernel = SqExponentialKernel(3.0)
# kernel = SqExponentialKernel([3.0,4.0])
Z_init = rand(M,D)*2 .-1.0

##
AGP.setadbackend(:reverse_diff)
grads_reverse = []
ELBO_reverse = []
l_reverse = []
σ_reverse = []
function cb_reverse(model,iter)
    if iter > 3
        push!(grads_reverse,AGP.grads[first(Flux.params(model.f[1].kernel))])
    end
    push!(ELBO_reverse,ELBO(model))
    push!(l_reverse,first(model.f[1].kernel.transform.s))
end
model = fullgp ? VGP(X,sign.(y),deepcopy(kernel),LogisticLikelihood(),AnalyticVI()) : SVGP(X,sign.(y),deepcopy(kernel),LogisticLikelihood(),AnalyticVI(),M)

if !fullgp
    model.f[1].Z.Z .= copy(Z_init)
end
@time train!(model,T,callback=cb_reverse)

##
AGP.setadbackend(:forward_diff)
grads_forward = []
gradσ_forward = []
ELBO_forward = []
l_forward = []
σ_forward = []
function cb_forward(model,iter)
    if iter > 3
        push!(grads_forward,AGP.grads[first(Flux.params(model.f[1].kernel))])
    end
    push!(ELBO_forward,ELBO(model))
    push!(l_forward,first(model.f[1].kernel.transform.s))
end
# model = OnlineSVGP(deepcopy(kernel),LogisticLikelihood(),AnalyticVI(),)
model = fullgp ? VGP(X,sign.(y),deepcopy(kernel),LogisticLikelihood(),AnalyticVI()) : SVGP(X,sign.(y),deepcopy(kernel),LogisticLikelihood(),AnalyticVI(),M)

if !fullgp
    model.f[1].Z.Z .= copy(Z_init)
end
@time train!(model,T,callback=cb_forward)
# @time train!(model,X,sign.(y),iterations=20)

##
using Plots;
pyplot()
default(lw=3.0)
y_pred,_ = proba_y(model,Xtest)

p = plot(layout=(3,2))
for j in 1:length(grads_forward[1])
    plot!(1:length(grads_forward),getindex.(grads_forward,j),lab="forward",subplot=1)
    plot!(1:length(grads_reverse),getindex.(grads_reverse,j),title="D=$j",linestyle=:dash,lab="reverse",subplot=1)
    plot!(1:length(ELBO_forward),ELBO_forward,lab="forward",subplot=2)
    plot!(1:length(ELBO_reverse),ls=:dash,ELBO_reverse,lab="reverse",title="ELBO",subplot=2)
    plot!(1:length(l_forward),l_forward,lab="forward",subplot=3)
    plot!(1:length(l_reverse),ls=:dash,l_reverse,lab="reverse",title="l",subplot=3,yaxis=:log)
    Plots.contourf!(collect(xgrid),collect(xgrid),reshape(y_pred,100,100),subplot=4,clims=(0.0,1.0))
    Plots.scatter!(eachcol(X)...,zcolor=sign.(y),subplot=4,lab="")
    if !fullgp
        scatter!(eachcol(Z_init)...,color=:blue,subplot=4,lab="")
    end
    plot!(1:length(σ_forward),σ_forward,lab="forward",subplot=5)
    plot!(1:length(σ_reverse),ls=:dash,σ_reverse,lab="reverse",title="σ",subplot=5,yaxis=:log)
    plot!(1:length(gradσ_forward),gradσ_forward,lab="forward",subplot=6)
    plot!(1:length(gradσ_reverse),ls=:dash,gradσ_reverse,lab="reverse",title="gradσ",subplot=6)
end
p |> display

Zygote.refresh()
g = Zygote.gradient(()->AGP.ELBO_given_theta(model),Flux.params(model.f[1].kernel))
AGP.update_hyperparameters!(model)
AGP.grads[first(Flux.params(model.f[1].kernel))][1]
g.grads[first(Flux.params(model.f[1].kernel))]
AGP.grads
g.grads

Flux.params(model.f[1].kernel)


text:
# This file was generated, do not modify it. # hide
r = report(mtm)
ss = r.parameter_values[:,1]
cbt = r.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(ss, cbt, r.measurements)

xlabel("Sub sample", fontsize=14)
ylabel("Col sample by tree", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

savefig("assets/literate/EX-crabs-xgb-heatmap2.svg") # hide


text:
#md# The word clouds generated by WordCloud.jl are always with a mask, but we can imitate the no-mask-style through the following steps:
#md# * set a lower density
#md# * set the background color as the mask color
#md# * gathering style placement
#md# * generating with teleporting off
using WordCloud
wc = wordcloud(
    processtext(open(pkgdir(WordCloud) * "/res/Donald Trump's Inaugural Address.txt"), maxweight=1, minweight=0),
    density=0.3,
    maskshape=box,
    masksize=(800, 600, 0),
    backgroundcolor=:maskcolor,
    angles=rand((0, (0, 90))),
    run=identity,
)
#md# place words in the center of the background to prevent encountering the mask boundary
placewords!(wc, style=:gathering, reorder=WordCloud.shuffle, level=6, rt=1) # a proper level is important, and so is luck
paint(wc, "nomask-placewords.svg")
#md# ![](nomask-placewords.svg)
#md# prevent teleporting words to the surrounding blank space
generate!(wc, teleporting=false)
paint(wc, "nomask.svg")
wc
#eval# runexample(:nomask)
#md# ![](nomask.svg)

max_stars_repo_name: guo-yong-zhi/WordCloud

max_issues_repo_name: guo-yong-zhi/WordCloud

max_forks_repo_name: guo-yong-zhi/WordCloud


text:
# ---
# title: 1610. Maximum Number of Visible Points
# id: problem1610
# author: Indigo
# date: 2021-02-17
# difficulty: Hard
# categories: Two Pointers, Geometry
# link: <https://leetcode.com/problems/maximum-number-of-visible-points/description/>
# hidden: true
# ---
#
# You are given an array `points`, an integer `angle`, and your `location`,
# where `location = [posx, posy]` and `points[i] = [xi, yi]` both denote
# **integral coordinates** on the X-Y plane.
#
# Initially, you are facing directly east from your position. You **cannot
# move** from your position, but you can **rotate**. In other words, `posx` and
# `posy` cannot be changed. Your field of view in **degrees** is represented by
# `angle`, determining how wide you can see from any given view direction. Let
# `d` be the amount in degrees that you rotate counterclockwise. Then, your
# field of view is the **inclusive** range of angles `[d - angle/2, d +
# angle/2]`.
#
# Your browser does not support the video tag or this video format.
#
# You can **see** some set of points if, for each point, the **angle** formed by
# the point, your position, and the immediate east direction from your position
# is **in your field of view**.
#
# There can be multiple points at one coordinate. There may be points at your
# location, and you can always see these points regardless of your rotation.
# Points do not obstruct your vision to other points.
#
# Return _the maximum number of points you can see_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/30/89a07e9b-00ab-4967-976a-c723b2aa8656.png)
#
#
#
#     Input: points = [[2,1],[2,2],[3,3]], angle = 90, location = [1,1]
#     Output: 3
#     Explanation: The shaded region represents your field of view. All points can be made visible in your field of view, including [3,3] even though [2,2] is in front and in the same line of sight.
#
#
# **Example 2:**
#
#
#
#     Input: points = [[2,1],[2,2],[3,4],[1,1]], angle = 90, location = [1,1]
#     Output: 4
#     Explanation: All points can be made visible in your field of view, including the one at your location.
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/30/5010bfd3-86e6-465f-ac64-e9df941d2e49.png)
#
#
#
#     Input: points = [[1,0],[2,1]], angle = 13, location = [1,1]
#     Output: 1
#     Explanation: You can only see one of the two points, as shown above.
#
#
#
#
# **Constraints:**
#
#   * `1 <= points.length <= 105`
#   * `points[i].length == 2`
#   * `location.length == 2`
#   * `0 <= angle < 360`
#   * `0 <= posx, posy, xi, yi <= 109`
#
#
## @lc code=start
using LeetCode

function max_visible_points(points::Vector{Vector{Int}}, view_angle::Int,
                            location::Vector{Int})
    points .-= Ref(location)
    len = length(points)
    angles = [atand(point...) for point in points if point != [0, 0]]
    res, samen = 0, len - length(angles)
    sort!(angles)
    angles = vcat(angles, angles .+ 360)
    for i in 1:length(angles)÷2
        idx = searchsortedlast(@view(angles[i+1:end]), angles[i] + view_angle)
        res = max(res, idx+1+samen)
    end
    res
end
## @lc code=end


text:
##
##  p o l y n o m . j l  Polynomials
##

function pval(p::Vector, x::Number)
    local n = length(p)
    if n == 0
        return NaN
    elseif n == 1
        return p[1]
    else
        y = p[1]
        for i in 2:n
            y = y * x + p[i]
        end
    end
    return y
end

function horner(p::Vector, x::Number)
    local n = length(p)
    if n == 0
        return NaN, NaN
    elseif n == 1
        return p[1], 0
    else
        y = p[1]; dy = 0
        for i in 2:n
            dy = dy * x + y
            y  =  y * x + p[i]
        end
    end
    return y, dy
end

function pzero{T<:Real}(p::Vector{T}, x0::T)
    local x = x0, tol = 2*eps(x)
    px, dpx = horner(p, x)
    df = -px/dpx
    niter = 1
    while abs(df) >= tol && niter <= 100
        x += df
        px, dpx = horner(p, x)
        df = -px/dpx
        niter += 1
    end
    if niter > 100
        warn("Number of iterations exceeded.")
        x = NaN
    end
    return x + df
end

function pfit()
    error("Polynomial fit not yet implemented.")
end


text:
# ===============================
# Written by AAE
# TU Delft, Winter 2014
# simulkade.com
# ===============================

# =============================== SOLVERS ===================================
function solveLinearPDE(m::MeshStructure, M::SparseMatrixCSC{Float64, Int64}, RHS::Array{Float64,1})
  N=m.dims
  x=M\RHS # until the problem is solved with Julia "\" solver
  phi = CellValue(m, reshape(x, tuple(N.+2...)))
  phi
end

function solveLinearPDE!(m::MeshStructure, M::SparseMatrixCSC{Float64, Int64}, RHS::Array{Float64,1}, phi::CellValue)
  N=m.dims
  x=M\RHS # until the problem is solved with Julia "\" solver
  phi.value = reshape(x, tuple(N.+2...))
end

function solvePDE(m::MeshStructure, M::SparseMatrixCSC{Float64, Int64}, RHS::Array{Float64,1})
  N=m.dims
  x=M\RHS # until the problem is solved with Julia "\" solver
  phi = CellValue(m, reshape(x, tuple(N.+2...)))
  return phi
end

function solveMUMPSLinearPDE(m::MeshStructure, M::SparseMatrixCSC{Float64, Int64}, RHS::Array{Float64,1})
  N = m.dims
  x = mumps_solver.solveMUMPS(M,RHS) # until the problem is solved with Julia "\" solver
  phi = CellValue(m, reshape(x, tuple(N.+2...)))
  return phi
  error("MUMPS needs to be installed and imported (import MUMPS).")
end

function solveExplicitPDE(phi_old::CellValue, dt::Real, RHS::Array{Float64,1},
  BC::BoundaryCondition)
  d = phi_old.domain.dimension
  N = phi_old.domain.dims
  phi_val=reshape(phi_old.value[:]+dt*RHS, tuple(N.+2...))
  if (d==1) || (d==1.5)
  	phi_val= phi_val[2:N[1]+1]
  elseif (d==2) || (d==2.5) || (d==2.8)
  	phi_val= phi_val[2:N[1]+1, 2:N[2]+1]
  elseif (d==3) || (d==3.2)
    phi_val= phi_val[2:N[1]+1, 2:N[2]+1, 2:N[3]+1]
  end
  return createCellVariable(phi_old.domain, phi_val, BC)
end

function solveExplicitPDE(phi_old::CellValue, dt::Real, RHS::Array{Float64,1},
  BC::BoundaryCondition, alfa::CellValue)
  d = phi_old.domain.dimension
  N = phi_old.domain.dims
  phi_val=reshape(phi_old.value[:]+dt*RHS./alfa.value[:], tuple(N.+2...))
  if (d==1) || (d==1.5)
  	phi_val= phi_val[2:N[1]+1]
  elseif (d==2) || (d==2.5) || (d==2.8)
  	phi_val= phi_val[2:N[1]+1, 2:N[2]+1]
  elseif (d==3) || (d==3.2)
    phi_val= phi_val[2:N[1]+1, 2:N[2]+1, 2:N[3]+1]
  end
  return createCellVariable(phi_old.domain, phi_val, BC)
end

# =========================== Visualization =================================
# function visualizeCells(phi::CellValue; vmin::Real =0.0, vmax::Real = 0.0, cmap="viridis", shading = "gouraud")
# d=phi.domain.dimension
# if d==1 || d==1.5
#   x = [phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
#   phi = [0.5*(phi.value[1]+phi.value[2]); phi.value[2:end-1]; 0.5*(phi.value[end-1]+phi.value[end])]
#   plot(x, phi)
# elseif d==2 || d==2.5
#   x = [phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
#   y = [phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
#   phi0 = Base.copy(phi.value)
#   phi0[:,1] = 0.5*(phi0[:,1]+phi0[:,2])
#   phi0[1,:] = 0.5*(phi0[1,:]+phi0[2,:])
#   phi0[:,end] = 0.5*(phi0[:,end]+phi0[:,end-1])
#   phi0[end,:] = 0.5*(phi0[end,:]+phi0[end-1,:])
#   phi0[1,1] = phi0[1,2]
#   phi0[1,end] = phi0[1,end-1]
#   phi0[end,1] = phi0[end,2]
#   phi0[end,end] = phi0[end,end-1]
#   if vmin == 0.0 && vmax == 0.0
#       vmin = minimum(phi0)
#       vmax = maximum(phi0)
#   end
#   pcolormesh(x, y, phi0', vmin = vmin, vmax = vmax, cmap = cmap, shading = shading)
# elseif d==2.8
#   x = [phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
#   y = [phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
#   phi0 = Base.copy(phi.value)
#   phi0[:,1] = 0.5*(phi0[:,1]+phi0[:,2])
#   phi0[1,:] = 0.5*(phi0[1,:]+phi0[2,:])
#   phi0[:,end] = 0.5*(phi0[:,end]+phi0[:,end-1])
#   phi0[end,:] = 0.5*(phi0[end,:]+phi0[end-1,:])
#   phi0[1,1] = phi0[1,2]
#   phi0[1,end] = phi0[1,end-1]
#   phi0[end,1] = phi0[end,2]
#   phi0[end,end] = phi0[end,end-1]
#   subplot(111, polar="true")
#   pcolor(y, x, phi0)
# elseif d==3
#     Nx = phi.domain.dims[1]
#     Ny = phi.domain.dims[2]
#     Nz = phi.domain.dims[3]
#     x=[phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
#     y=zeros(1,Ny+2)
#     y[:]=[phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
#     z=zeros(1,1,Nz+2)
#     z[:]=[phi.domain.facecenters.z[1]; phi.domain.cellcenters.z; phi.domain.facecenters.z[end]]

#     phi0 = Base.copy(phi.value)
#     phi0[:,1,:]=0.5*(phi0[:,1,:]+phi0[:,2,:])
#     phi0[:,end,:]=0.5*(phi0[:,end-1,:]+phi0[:,end,:])
#     phi0[:,:,1]=0.5*(phi0[:,:,1]+phi0[:,:,1])
#     phi0[:,:,end]=0.5*(phi0[:,:,end-1]+phi0[:,:,end])
#     phi0[1,:,:]=0.5*(phi0[1,:,:]+phi0[2,:,:])
#     phi0[end,:,:]=0.5*(phi0[end-1,:,:]+phi0[end,:,:])

#     vmin = minimum(phi0)
#     vmax = maximum(phi0)
#     mynormalize = a -> (a .- vmin)./(vmax-vmin)

#     a=ones(Nx+2,Ny+2,Nz+2)
#     X = x.*a
#     Y = y.*a
#     Z = z.*a

#     fig = figure()
#     ax = fig[:add_subplot](111, projection = "3d")
#     # r = linspace(1.25, 1.25, 50)
#     # p = linspace(0, 2π, 50)
#     # R = repeat(r, 1, 50)
#     # P = repeat(p', 50, 1)
#     # Zc = rand(50, 50) # (P.^2-1).^2
#     # Z = repeat(linspace(0, 2, 50), 1, 50)
#     # X, Y = R.*cos.(P), R.*sin.(P)
#     ax[:plot_surface](X[1,:,:], Y[1,:,:], Z[1,:,:], facecolors=PyPlot.cm[:viridis](mynormalize(phi0[1,:,:])), alpha=0.8)
#     ax[:plot_surface](X[end,:,:], Y[end,:,:], Z[end,:,:], facecolors=PyPlot.cm[:viridis](mynormalize(phi0[end,:,:])), alpha=0.8)
#     ax[:plot_surface](X[:,1,:], Y[:,1,:], Z[:,1,:], facecolors=PyPlot.cm[:viridis](mynormalize(phi0[:,1,:])), alpha=0.8)
#     ax[:plot_surface](X[:,end,:], Y[:,end,:], Z[:,end,:], facecolors=PyPlot.cm[:viridis](mynormalize(phi0[:,end,:])), alpha=0.8)
#     ax[:plot_surface](X[:,:,1], Y[:,:,1], Z[:,:,1], facecolors=PyPlot.cm[:viridis](mynormalize(phi0[:,:,1])), alpha=0.8)
#     ax[:plot_surface](X[:,:,end], Y[:,:,end], Z[:,:,end], facecolors=PyPlot.cm[:viridis](mynormalize(phi0[:,:,end])), alpha=0.8)
#   # s=mayavis.pipeline[:scalar_field](X,Y,Z,phi0)
#   #
#   # mayavis.pipeline[:image_plane_widget](s, plane_orientation="x_axes", slice_index=0, vmin=vmin, vmax=vmax)
#   # mayavis.pipeline[:image_plane_widget](s, plane_orientation="y_axes", slice_index=0, vmin=vmin, vmax=vmax)
#   # mayavis.pipeline[:image_plane_widget](s, plane_orientation="z_axes", slice_index=0, vmin=vmin, vmax=vmax)
#   # mayavis.pipeline[:image_plane_widget](s, plane_orientation="z_axes", slice_index=floor(Integer,Nz/2.0), vmin=vmin, vmax=vmax)
#   # mayavis.outline()

# #   # 6 surfaces
# #   # surfaces 1,2 (x=x[1], x=x[end])
# #   Y=repeat(y,1,Nz)
# #   Z=repeat(z,1,Ny)
# #   mayavis.mesh(x[1]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[2,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.mesh(x[end]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[end-1,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
# #
# #   # surfaces 3,4 (y=y[1], y=y[end]
# #   X = repeat(x,1,Nz)
# #   Z = repeat(z,1,Nx)
# #   mayavis.mesh(X,y[1]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,2,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.mesh(X,y[end]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,end-1,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.axes()
# #
# #   # surfaces 5,6 (z=z[1], z=z[end]
# #   X = repeat(x,1,Ny)
# #   Y = repeat(y,1,Nx)
# #   mayavis.mesh(X,Y',z[1]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,2], vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.mesh(X,Y',z[end]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,end-1], vmin=vmin, vmax=vmax, opacity=0.8)

#   # mayavis.colorbar()
#   # mayavis.axes()
#   # mshot= mayavis.screenshot()
#   # mayavis.show()
#   # return mshot

# elseif d==3.2
#   Nx = phi.domain.dims[1]
#   Ny = phi.domain.dims[2]
#   Nz = phi.domain.dims[3]
#   r=[phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
#   theta = zeros(1,Ny+2)
#   theta[:]=[phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
#   z=zeros(1,1,Nz+2)
#   z[:]=[phi.domain.facecenters.z[1]; phi.domain.cellcenters.z; phi.domain.facecenters.z[end]]
#   a=ones(Nx+2,Ny+2,Nz+2)
#   R=r.*a
#   TH = theta.*a
#   Z = z.*a

#   X=R.*cos(TH)
#   Y=R.*sin(TH)

#   phi0 = Base.copy(phi.value)
#   phi0[:,1,:]=0.5*(phi0[:,1,:]+phi0[:,2,:])
#   phi0[:,end,:]=0.5*(phi0[:,end-1,:]+phi0[:,end,:])
#   phi0[:,:,1]=0.5*(phi0[:,:,1]+phi0[:,:,1])
#   phi0[:,:,end]=0.5*(phi0[:,:,end-1]+phi0[:,:,end])
#   phi0[1,:,:]=0.5*(phi0[1,:,:]+phi0[2,:,:])
#   phi0[end,:,:]=0.5*(phi0[end-1,:,:]+phi0[end,:,:])

#   vmin = minimum(phi0)
#   vmax = maximum(phi0)
#   # 6 surfaces
#   # surfaces 1,2 (x=x[1], x=x[end])
#   mayavis.mesh(squeeze(X[floor(Integer,Nx/2.0),:,:],1),squeeze(Y[floor(Integer,Nx/2.0),:,:],1),squeeze(Z[floor(Integer,Nx/2.0),:,:],1),
#     scalars=squeeze(phi0[floor(Integer,Nx/2.0)+1,:,:],1), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(squeeze(X[Nx,:,:],1),squeeze(Y[Nx,:,:],1),squeeze(Z[Nx,:,:],1),
#     scalars=squeeze(phi0[Nx+2,:,:],1), vmin=vmin, vmax=vmax, opacity=0.8)

#   # surfaces 3,4 (y=y[1], y=y[end]
#   mayavis.mesh(squeeze(X[:,floor(Integer,Ny/2.0),:],2),squeeze(Y[:,floor(Integer,Ny/2.0),:],2),squeeze(Z[:,floor(Integer,Ny/2.0),:],2),
#     scalars=squeeze(phi0[:,floor(Integer,Ny/2.0)+1,:],2), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(squeeze(X[:,Ny,:],2),squeeze(Y[:,Ny,:],2),squeeze(Z[:,Ny,:],2),
#     scalars=squeeze(phi0[:,Ny+2,:],2), vmin=vmin, vmax=vmax, opacity=0.8)

#   # surfaces 5,6 (z=z[1], z=z[end]
#   mayavis.mesh(X[:,:,floor(Integer,Nz/2.0)],Y[:,:,floor(Integer,Nz/2.0)],Z[:,:,floor(Integer,Nz/2.0)],
#     scalars=phi0[:,:,floor(Integer,Nz/2.0)+1], vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(X[:,:,Nz],Y[:,:,Nz],Z[:,:,Nz],
#     scalars=phi0[:,:,Nz+1], vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.colorbar()
#   mayavis.axes()
#   mshot=mayavis.screenshot()
#   mayavis.show()
#   return mshot
# end
# end

# ########################## Visualize Vectors #####################
# function visualizeCellVectors(phi::CellVector)
# d=phi.domain.dimension
# if d==1 || d==1.5
#   println("Vector visualization works only in 2D and 3D")
# elseif d==2 || d==2.5
#   x = phi.domain.cellcenters.x
#   y = phi.domain.cellcenters.y
#   quiver(repeat(x, 1, length(y)), repeat(y', length(x), 1), phi.xvalue, phi.yvalue)
# elseif d==2.8
#   x = phi.domain.cellcenters.x
#   y = phi.domain.cellcenters.y'
#   subplot(111, polar="true")
#   quiver(repeat(y, length(x), 1), repeat(x, 1, length(y)),
#   phi.xvalue.*cos(y)-phi.yvalue.*sin(y), phi.xvalue.*sin(y)+phi.yvalue.*cos(y))
# elseif d==3
#   Nx = phi.domain.dims[1]
#   Ny = phi.domain.dims[2]
#   Nz = phi.domain.dims[3]
#   x=phi.domain.cellcenters.x
#   y=zeros(1,Ny)
#   y[:]=phi.domain.cellcenters.y
#   z=zeros(1,1,Nz)
#   z[:]=phi.domain.cellcenters.z

#   #vmin = minimum(phi.xvalue)
#   #vmax = maximum(phi0)

#   a=ones(Nx,Ny,Nz)
#   X = x.*a
#   Y = y.*a
#   Z = z.*a

#   mayavis.quiver3d(X,Y,Z, phi.xvalue, phi.yvalue, phi.zvalue)
#   mayavis.outline()

# #   # 6 surfaces
# #   # surfaces 1,2 (x=x[1], x=x[end])
# #   Y=repeat(y,1,Nz)
# #   Z=repeat(z,1,Ny)
# #   mayavis.mesh(x[1]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[2,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.mesh(x[end]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[end-1,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
# #
# #   # surfaces 3,4 (y=y[1], y=y[end]
# #   X = repeat(x,1,Nz)
# #   Z = repeat(z,1,Nx)
# #   mayavis.mesh(X,y[1]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,2,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.mesh(X,y[end]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,end-1,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.axes()
# #
# #   # surfaces 5,6 (z=z[1], z=z[end]
# #   X = repeat(x,1,Ny)
# #   Y = repeat(y,1,Nx)
# #   mayavis.mesh(X,Y',z[1]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,2], vmin=vmin, vmax=vmax, opacity=0.8)
# #   mayavis.mesh(X,Y',z[end]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,end-1], vmin=vmin, vmax=vmax, opacity=0.8)

#   mayavis.colorbar()
#   mayavis.axes()
#   mshot= mayavis.screenshot()
#   mayavis.show()
#   return mshot

# elseif d==3.2
#   Nx = phi.domain.dims[1]
#   Ny = phi.domain.dims[2]
#   Nz = phi.domain.dims[3]
#   r=phi.domain.cellcenters.x
#   theta = zeros(1,Ny)
#   theta[:]=phi.domain.cellcenters.y
#   z=zeros(1,1,Nz)
#   z[:]=phi.domain.cellcenters.z
#   a=ones(Nx,Ny,Nz)
#   R=r.*a
#   TH = theta.*a
#   Z = z.*a

#   X=R.*cos(TH)
#   Y=R.*sin(TH)

#   #vmin = minimum(phi0)
#   #vmax = maximum(phi0)
#   # 6 surfaces
#   # surfaces 1,2 (x=x[1], x=x[end])
#   mayavis.quiver3d(X,Y,Z, phi.xvalue.*cos(TH)-phi.yvalue.*sin(TH),
#   phi.xvalue.*sin(TH)+phi.yvalue.*cos(TH), phi.zvalue)
#   mayavis.colorbar()
#   mayavis.axes()
#   mshot=mayavis.screenshot()
#   mayavis.show()
#   return mshot
# end
# end


text:
"""
Module to simplify working with events and filters.
--jleugeri, 2017

Usage example:
==============

julia```
using Plots
pyplot()
e1 = EventTrain([1.0, 2.0, 3.0])
e2 = 2*EventTrain([0.5, 1.5, 2.5, 3.5])
e3 = 1.5 * e1 + e2

F = LaplaceFilter(s->1/(s+1))
X = F∘F∘e3
M = x - F∘x
plot(X, 0,10, color=:blue, event_color=:blue, show_events=false)
plot!(F∘x, 0, 10, color=:red, show_events=false)
plot!(M, 0, 10, color=:black)
'''
"""
module EventFiltering
using InverseLaplace, Plots

export Filter, LinearFilterCombination, LaplaceFilter, GammaFilter, ExponentialFilter, EventTrain, FilteredEventTrain, freeze, fromto, ℒ, delay
    # Base.:∘,
    # Base.:+,
    # Base.:-,
    # Base.:*,
    # Base.length,
    # Base.getindex,

# Helper function to evaluate a list of functions on given arguments
call_all(fs, args...) = [f(args...) for f in fs]

###########
# Filters #
###########

abstract type Filter{T, E} end

### LaplaceFilter
struct LaplaceFilter{T} <: Filter{T, Float64}
    terms::Vector{<:Function}
end
LaplaceFilter(t::Vector{<:Function}) = LaplaceFilter{Float64}(t)
LaplaceFilter(t::Function) = LaplaceFilter{Float64}(Function[t])

# Convolve LaplaceFilters by convolution theorem
Base.:∘(F1::LaplaceFilter, F2::LaplaceFilter) = LaplaceFilter([F1.terms; F2.terms])

# Evaluate LaplaceFilter
function ℒ(F::LaplaceFilter)
    function LLaplace(s::V)::V where V
        prod(call_all(F.terms,s))
    end
end

function freeze(F::LaplaceFilter{T}; method=talbot, accuracy=32) where T
    ilt = ILt(ℒ(F), method, accuracy)
    function frozen(x::T)::T
        ilt(x)
    end
end

((F::LaplaceFilter{T})(t::T)::T) where T = freeze(F)(t)

### GammaFilter
struct GammaFilter{T,A,K} <: Filter{T, Float64}
    α::A
    k::K
end
GammaFilter(α::A, k::K) where {A,K} = GammaFilter{Float64,A,K}(α, k)

# Convolve GammaFilters (convert to LaplaceFilter if necessary)
function Base.:∘(F1::GammaFilter, F2::GammaFilter)
    if F1.α == F2.α
        GammaFilter(F1.α, F1.k+F2.k)
    else
        LaplaceFilter(F1)∘LaplaceFilter(F2)
    end
end

# Evaluate GammaFilter
function ℒ(F::GammaFilter)
    function Lgamma(s::V)::V where V
        (F.α/(s+F.α))^F.k
    end
end

freeze(F::GammaFilter) = F
((F::GammaFilter{T,A,K})(t::T)::promote_type(T,A,K)) where {T,A,K} = (F.α^F.k)/gamma(F.k)*t^(F.k-1)*exp(-F.α*t)

# Convert GammaFilter to LaplaceFilter
LaplaceFilter(F::GammaFilter) = LaplaceFilter([ℒ(F)])

# ExponentialFilter
ExponentialFilter(α=1.0) = GammaFilter(α, 1)

struct LinearFilterCombination{Ws, Fs} <: Filter{Float64, Float64}
    a::Vector{Ws}
    F::Vector{Fs}
end

function freeze(F::LinearFilterCombination)
    frozen_funs = map(freeze, F.F)
    function (frozen(x::T)::T) where T
        F.a⋅call_all(frozen_funs, x)
    end
end

function ℒ(F::LinearFilterCombination)
    ℒs = map(ℒ, F.F)
    function (LLFC(s::V)::V) where V
        F.a⋅call_all(ℒs, s)
    end
end

((F::LinearFilterCombination)(t::T)::T) where T = freeze(F)(t)

# Convert Filter to LinearFilterCombination
LinearFilterCombination(F::LinearFilterCombination) = F
LinearFilterCombination(F::Filter) = LinearFilterCombination([1.0],[F])

# Add LinearFilterCombinations
Base.:+(F1::LinearFilterCombination, F2::LinearFilterCombination) = LinearFilterCombination([F1.a;F2.a],[F1.F;F2.F])
# Scale LinearFilterCombinations
Base.:*(a::Real, F::LinearFilterCombination) = LinearFilterCombination(F.a*a, F.F)

# Convolve LinearFilterCombinations by linearity
Base.:∘(F1::LinearFilterCombination, F2::LinearFilterCombination) = LinearFilterCombination([a*b for a ∈ F1.a for b ∈ F2.a], [f∘g for f ∈ F1.F for g ∈ F2.F])
Base.:∘(F1::Filter, F2::LinearFilterCombination) = LinearFilterCombination(F1) ∘ F2
Base.:∘(F1::LinearFilterCombination, F2::Filter) = F2∘F1

# Fallback methods enabling linearity for all filter types
Base.:+(F1::Filter, F2::Filter) = LinearFilterCombination(F1)+LinearFilterCombination(F2)
# Scale Filters
Base.:*(a::Real, F::Filter) = LinearFilterCombination([a], [F])
Base.:*(F::Filter, a::Real) = LinearFilterCombination([a], [F])
# Sign-Invert Filters
Base.:-(F1::Filter) = (-1) * F1
# Subtract Filters
Base.:-(F1::Filter, F2::Filter) = F1 + (-F2)

#####################################
# EventTrain and FilteredEventTrain #
#####################################

struct EventTrain{T,E}
    times::Vector{T}
    events::E
    function EventTrain(times::Vector{T}, events::E = 1) where {T,E}
        if ~issorted(times)
            ord = sortperm(times)
            times = times[ord]

            if isa(events, Vector)
                events = events[ord]
            end
        end
        new{T,E}(times, events)
    end
end

struct FilteredEventTrain{T,E,F<:Filter}
    event_train::EventTrain{T,E}
    F::F
end

#FilteredEventTrain(e::EventTrain{T,E,F}, F::F) where {T,E,F<:Filter} = FilteredEventTrain{T,E,F}(e, F)

Base.length(s::EventTrain) = length(s.times)
Base.length(s::FilteredEventTrain) = length(s.event_train)
Base.endof(s::EventTrain) = length(s)
Base.endof(s::FilteredEventTrain) = length(s)

"""Merge two event trains (times must be SORTED!)"""
function Base.:+(s::EventTrain{S,E}, t::EventTrain{T, F}) where {S,T,E,F}
    ST = promote_type(S,T)
    new_times = ST[]

    EF = if (E<:Vector) | (F<:Vector)
        # We need a vector to hold the events
        eltype_E = E<:Vector ? eltype(E) : E
        eltype_F = F<:Vector ? eltype(F) : F
        Vector{promote_type(eltype_E,eltype_F)}
    elseif s.events != t.events
        # the events of s and t are different and must thus be stored in a vector
        Vector{promote_type(E,F)}
    else
        # We don't need a vector to hold the events
        promote_type(E,F)
    end

    new_events = (EF<:Vector) ? EF() : EF(s.events)

    (i1,i2) = (1,1)
    (n1,n2) = length.((s, t))
    while (i1 <= n1) & (i2 <= n2)
        if(s.times[i1] < t.times[i2])
            push!(new_times, s.times[i1])

            if EF<:Vector
                push!(new_events, (E<:Vector) ? s.events[i1] : s.events)
            end
            i1 += 1
        else
            push!(new_times, t.times[i2])
            if EF<:Vector
                push!(new_events, (F<:Vector) ? t.events[i2] : t.events)
            end
            i2 += 1
        end
    end

    while i1 <= n1
        push!(new_times, s.times[i1])
        if EF<:Vector
            push!(new_events,(E<:Vector) ? s.events[i1] : s.events)
        end
        i1 += 1
    end

    while i2 <= n2
        push!(new_times, t.times[i2])
        if EF<:Vector
            push!(new_events,(F<:Vector) ? t.events[i2] : t.events)
        end
        i2 += 1
    end

    return EventTrain(new_times, new_events)
end
# Add FilteredEventTrains
function Base.:+(s1::FilteredEventTrain,s2::FilteredEventTrain)
    return if s1.F == s2.F
        FilteredEventTrain(s1.event_train+s2.event_train, s1.F)
    elseif s1.event_train == s2.event_train
        FilteredEventTrain(s1.event_train, s1.F+s2.F)
    elseif s1.event_train == -s2.event_train
        FilteredEventTrain(s1.event_train, s1.F-s2.F)
    else
        x->(s1(x)+s2(x))
    end
end

# Sign-Invert EventTrains
Base.:-(s1::EventTrain) = (-1) * s1
# Sign-Invert FilteredEventTrains
Base.:-(s1::FilteredEventTrain) = (-1) * s1

# Subtract EventTrains
Base.:-(s1::EventTrain, s2::EventTrain) = s1 + (-s2)
# Subtract FilteredEventTrains
Base.:-(s1::FilteredEventTrain, s2::FilteredEventTrain) = s1 + (-s2)

# Scale EventTrain
Base.:*(a, s::EventTrain) = EventTrain(s.times, a.*s.events)
# Scale FilteredEventTrain
Base.:*(a, s::FilteredEventTrain) = FilteredEventTrain(a.*s.event_train, s.F)

# Subsample events of a (Filtered)EventTrain
Base.getindex(s::EventTrain{T,<:Vector}, r) where {T} = EventTrain(s.times[r], s.events[r])
Base.getindex(s::EventTrain, r) = EventTrain(s.times[r], s.events)
Base.getindex(s::FilteredEventTrain, r) = FilteredEventTrain(s.event_train[r], s.F)

# Slice a time interval out of a (Filtered)EventTrain
function fromto(s::EventTrain, from, to)
    if from == :start
        from = minimum(s.times)-1.0
    end
    if to == :end
        to = maximum(s.times)+1.0
    end

    r = searchsortedfirst(s.times, from):searchsortedlast(s.times, to)
    return s[r]
end
fromto(s::FilteredEventTrain, from, to) = FilteredEventTrain(fromto(s.event_train, from, to), s.F)

# Shift a (Filtered)EventTrain in time
function delay(s::EventTrain{T, E}, dt) where {T,E}
    new_times = s.times+dt
    ord = sortperm(new_times)
    EventTrain(new_times[ord], (E<:Vector) ? s.events[ord] : s.events)
end

delay(s::FilteredEventTrain, dt) = FilteredEventTrain(delay(s.event_train, dt), s.F)

# Convolve two EventTrains
Base.:∘(e1::EventTrain{S, E}, e2::EventTrain{T, F}) where {S,T,E,F} = ((E<:Vector) | (F<:Vector)) ? mapreduce(dtw->delay(dtw[2]*e1, dtw[1]), +, zip(e2.times,e2.events)) : mapreduce(dt->delay(e2.events*e1, dt), +, e2.times)

# Convolve two FilteredEventTrains
Base.:∘(e1::FilteredEventTrain, e2::FilteredEventTrain) = FilteredEventTrain(e1.event_train∘e2.event_train, e1.F∘e2.F)
# Convolve one FilteredEventTrain and one EventTrain
Base.:∘(e1::FilteredEventTrain, e2::EventTrain) = FilteredEventTrain(e1.event_train∘e2, e1.F)
Base.:∘(e1::EventTrain, e2::FilteredEventTrain) = e2∘e1

################################################
# Interaction between filters and event trains #
################################################

# Use linearity to chain filters
Base.:∘(F::Filter, e::FilteredEventTrain) = FilteredEventTrain(e.event_train, F∘e.F)

# Apply Filter to event-train
Base.:∘(F::Filter, e::EventTrain) = FilteredEventTrain(e, F)

# Evaluate FilteredEventTrain
function ((s::FilteredEventTrain{T,E,F})(t::T, args...)) where {T,E,F}
    s_frozen = freeze(s, args...)
    return s_frozen(t)
end

function freeze(s::FilteredEventTrain{T,E,F}, max_range=nothing) where {T,E,F}
    F_frozen = freeze(s.F)

    function frozen(t::T)
        sub_train = fromto(s.event_train, max_range==nothing ? :start : t-max_range, t)
        return sum(F_frozen.(t.-sub_train.times).*sub_train.events)
    end
end

####################
# Plotting recipes #
####################

@recipe function f(f::Filter, xmin, xmax)
    x->f(x), xmin, xmax
end

@recipe function f(fv::Vector{Filter}, xmin, xmax)
    [x->f(x) for f in fv], xmin, xmax
end

@recipe function f(e::EventTrain, y=0.0)
    y = (length(y) == length(e.times)) ? y : fill(y, length(e.times))
    ([e.times'; e.times'; e.times'][:], [y'; (y+e.events)'; fill(NaN, length(e.times))'][:])
end

@recipe function f(e::FilteredEventTrain{T,E,F}, xmin=:start, xmax=:end; show_events=false, show_kernels=false, stack_kernels=true, event_color=:black) where {T,E,F}
    if xmax==:end
        xmax = maximum(e.event_train.times)
    end
    if xmin==:start
        xmin = minimum(e.event_train.times)
    end

    xlims := (xmin, xmax)

    if show_events
        @series begin
            label := ""
            color := event_color
            e.event_train
        end
    end

    if show_kernels
        f = freeze(e.F)
        sub = fromto(e, :start, xmax)

        if stack_kernels
            for t ∈ sub.event_train.times
                if t<xmax
                    @series begin
                        linewidth := 1
                        linestyle := :dot
                        color := event_color
                        label := ""
                        (x->convert(Float64, fromto(sub, :start, t)(x))), t, xmax
                    end
                end
            end
        elseif (E<:Vector)
            for (t,w) ∈ zip(sub.event_train.times, sub.event_train.events)
                if t<xmax
                    @series begin
                        linewidth := 1
                        linestyle := :dot
                        color := event_color
                        label := ""
                        (x->convert(Float64, x>=t ? (w*f(x-t)) : 0.0)), t, xmax
                    end
                end
            end
        else
            for t ∈ sub.event_train.times
                if t<xmax
                    @series begin
                        linewidth := 1
                        linestyle := :dot
                        color := event_color
                        label := ""
                        (x->convert(Float64, x>=t ? (sub.event_train.events*f(x-t)) : 0.0)), t, xmax
                    end
                end
            end
        end
    end

    @series begin
        label --> ""
        x->e(x), xmin, xmax
    end
end
end


text:
import Mustache

template = """
#=
This file is auto-generated. Do not edit.
=#
module Input

import PowerSystems
const PSY = PowerSystems
import EMISApproximateEquilibrium
const SR = EMISApproximateEquilibrium

devices_added = SR.ContinuousRegionalBuildContainer()
{{#types}}
devices_added[PSY.PrimeMovers.{{type}}] = {{{config}}}
{{/types}}
devices_deleted  = SR.ConfigDataContainer()

devices_deleted[PSY.ThermalStandard] = [
    "115_STEAM_1",
    "115_STEAM_2",
    "315_STEAM_1",
    "315_STEAM_2",
    "315_STEAM_3",
    "315_STEAM_4",
    "315_STEAM_5",
    "107_CC_1",
    "113_CT_1",
    "113_CT_2",
    "113_CT_3",
    "113_CT_4",
    "207_CT_2",
    "213_CC_3",
    "213_CT_1",
    "213_CT_2",
    "307_CT_2",
    "313_CC_1",
    "315_CT_6",
    "315_CT_7",
    "315_CT_8",
    "318_CC_1",
    "323_CC_1",
    "323_CC_2",
    "218_CC_1",
    "118_CC_1",
    "101_CT_1",
    "101_CT_2",
    "102_CT_1",
    "102_CT_2",
    "201_CT_1",
    "201_CT_2",
    "202_CT_1",
    "202_CT_2",
    "301_CT_1",
    "301_CT_2",
    "302_CT_1",
    "302_CT_2",
    "207_CT_1",
    "307_CT_1",
    "101_STEAM_4",
    "214_SYNC_COND_1",
    "114_SYNC_COND_1",
]

devices_deleted[PSY.RenewableFix] = [
    "308_RTPV_1",
    "313_RTPV_1",
    "313_RTPV_2",
    "313_RTPV_3",
    "313_RTPV_4",
    "313_RTPV_5",
    "313_RTPV_6",
    "313_RTPV_7",
    "313_RTPV_8",
    "313_RTPV_9",
    "313_RTPV_10",
    "313_RTPV_11",
    "313_RTPV_12",
    "320_RTPV_1",
    "320_RTPV_2",
    "320_RTPV_3",
    "313_RTPV_13",
    "320_RTPV_4",
    "320_RTPV_5",
    "118_RTPV_1",
    "118_RTPV_2",
    "118_RTPV_3",
    "118_RTPV_4",
    "118_RTPV_5",
    "118_RTPV_6",
    "320_RTPV_6",
    "118_RTPV_7",
    "118_RTPV_8",
    "118_RTPV_9",
    "118_RTPV_10",
    "213_RTPV_1",
]

devices_deleted[PSY.RenewableDispatch] = ["309_WIND_1", "212_CSP_1", "317_WIND_1", "122_WIND_1",
"313_PV_1", "314_PV_3","313_PV_2", "312_PV_1", "319_PV_1","215_PV_1", "212_CSP_1"]
devices_deleted[PSY.HydroDispatch] = [
    "122_HYDRO_5",
    "122_HYDRO_6",
    "201_HYDRO_4",
    "215_HYDRO_3",
    "222_HYDRO_3",
    "222_HYDRO_4",
    "222_HYDRO_5",
    "222_HYDRO_6",
]
end
"""

function generate_regional_semi_discrete_configuration(directory, future_config_dict)
    for (id, config) in future_config_dict
        item = Dict()
        params = Vector{Dict}()
        for (type, devices) in config
            if !isempty(devices)
                tmptype = Dict("type" => type, "config" => devices)
                push!(params, tmptype)
            end
        end
        item["types"] = params
        filename = joinpath(directory, "sample_configuration_$(id).jl")
        open(filename, "w") do io
            write(io, Mustache.render(template, item))
        end
        println("Wrote $filename")
    end
    return
end


text:
abstract type AbstractRecorder end
abstract type StateRecorder <: AbstractRecorder end

mutable struct WealthRecorder{T<:Real} <: StateRecorder
    record::DenseAxisArray{T,1}

    function WealthRecorder(dates; zero::T=0.0) where {T}
        date_step = dates[end] - dates[end-1]
        new{T}(DenseAxisArray{T}(undef, vcat(dates,dates[end]+date_step)))
    end
end

function record!(recorder::WealthRecorder, market::VolumeMarket, date)
    recorder.record[date] = market_budget(market)
    return nothing
end

get_records(recorder::StateRecorder) = recorder.record

mutable struct ReturnsRecorder{T<:Real} <: StateRecorder
    record::DenseAxisArray{T,1}

    function ReturnsRecorder(dates; zero::T=0.0) where {T}
        new{T}(DenseAxisArray{T}(undef, dates))
    end
end

function record!(recorder::ReturnsRecorder, market::VolumeMarket, date)
    recorder.record[date] = total_profit(market) / market_budget(market)
    return nothing
end

mutable struct DecisionRecorder{T<:Real} <: StateRecorder
    record::DenseAxisArray{T}

    function DecisionRecorder(numA::Int, dates; zero::T=0.0) where {T}
        new{T}(DenseAxisArray{T}(undef, dates, 1:numA))
    end
end

function record!(recorder::DecisionRecorder, market::VolumeMarket, date)
    recorder.record[date, :] = market.volume_bids
    return nothing
end

default_state_recorders(numA::Int, date_range) = Dict(
    :wealth => WealthRecorder(date_range),
    :returns => ReturnsRecorder(date_range),
    :decisions => DecisionRecorder(numA, date_range)
)

function change_bids!(market::VolumeMarket, formulation::PortfolioFormulation, optimizer_factory::Any; record_measures::Bool=false)
    model, decision_variables = market_model(market, optimizer_factory; sense=sense(formulation))
    pointers = portfolio_model!(model, formulation, decision_variables; record_measures=record_measures)
    change_bids!(market, model, decision_variables)

    status = termination_status(model)
    if status !== MOI.OPTIMAL
        return nothing
    else
        return pointers
    end
end

function day_backtest_market!(
    strategy_logic::Function,
    market_history::MarketHistory,
    date;
    state_recorders::Dict{Symbol,StateRecorder}=Dict(),
    optimization_recorder::Union{DenseAxisArray,Nothing}=nothing,
    provide_state::Bool=false
)
    market =  market_template(market_history, date)
    pointers = if provide_state
        strategy_logic(
            market,
            past_prices(market_history, date),
            Dict(:state_recorders => state_recorders, :date => date)
        )
    else
        strategy_logic(
            market,
            past_prices(market_history, date)
        )
    end

    clear_market!(market, current_prices(market_history, date))

    for recorder in values(state_recorders)
        record!(recorder, market, date)
    end

    if !isnothing(pointers)
        for (uid, pointer) in pointers
            optimization_recorder[date,uid] = value(pointer)
        end
    end

    set_market_budget(market, market_budget(market) + total_profit(market))

    return nothing
end

"""
    sequential_backtest_market(strategy_logic::Function, market_history::VolumeMarketHistory) -> Vector{Real}, Vector{Real}

Simple sequential backtest functionality for strategies that return an array of invested money per asset. Ideal for statefull strategies.

API:
```julia
wealth_strategy, returns_strategy = sequential_backtest_market(market_history; date_range=date_range)
    do market, past_returns
    # ... Strategy definition ...
    return investment_decision
end
```

Arguments:
 - `strategy_logic::Function`: Function that represents the investment strategy.
 - `market_history::MarketHistory`: Market history.

Optional Keywork Arguments:
 - `date_range`: Dates to be simulated with corresponding entries in the market history.
"""
function sequential_backtest_market(
    strategy_logic::Function,
    market_history::MarketHistory,
    date_range;
    optimization_recorder::Union{DenseAxisArray,Nothing}=nothing,
    state_recorders::Dict{Symbol,StateRecorder}=default_state_recorders(number_assets(market_history), date_range)
)
    @assert issorted(date_range)

    for date in date_range
        day_backtest_market!(strategy_logic, market_history, date;
            state_recorders=state_recorders, optimization_recorder=optimization_recorder, provide_state=true
        )
    end

    if haskey(state_recorders, :wealth)
        date_step = date_range[end] - date_range[end-1]
        record!(state_recorders[:wealth], market_template(market_history), date_range[end]+date_step)
    end

    return state_recorders, optimization_recorder
end

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
"""
Plots a subgraph that is part of a greater one
global ids in property :eid.
"""
function plot_subg(sg)
    return graphplot(sg
                    ,names = [ get_prop(sg,i,:eid) for i=1:nv(sg)]
          )
end
"""
 Plots a multiagent control network
"""
function plot_ctl_network_multi( model; kwargs...,)

    nsize = 0.13
    lwidth = 0.5

    method = model.ctrl_model == GraphModel(4) ? :sfdp : :circular
    Random.seed!(model.seed)

    ctl_p = graphplot(
        model.ctl_graph
        ,names = [ i for i in 1:nv(model.ctl_graph) ]
                #[ get_control_agent(i,model) for i in 1:nv(model.ctl_graph) ]
        ,method = method#:sfdp#:stress#:shell# #:spectral #:circular
        #TODO check if required, not working atm
        #,func = NetworkLayout.SFDP.layout(adjacency_matrix(model.ctl_graph),2)
        , curvature_scalar = 0.0
        ,size=(300,200)
        ,node_weights = [ i > 9 ? 1 : 5 for i in 1:nv(model.ctl_graph)]
        ,nodeshape = :circle
        # ,nodecolor = [ has_active_controlled_assets(
        #                 getindex(model,model.ctl_graph[i,:aid]),model
        #                ) ? :lightblue : :lightgray for i in 1:nv(model.ctl_graph) ]
        # ,markerstrokecolor = :dimgray
        # ,edgecolor=:dimgray
        ,markerstrokewidth = 1.1
        ,node_size=nsize
        ,edgestyle = :dot
        ,titlefontcolor=:white
        ,curves = false
    )
    #TODO replace buggy annotation not thread-safe
    # annotate!((-0.7,0.72,Plots.text("Control Network", 11, :black, :center)))
    annotate!((0.0,0.72,Plots.text("Control Network", 11, :black, :center)))

    return ctl_p

end
"""
 Plots a single agent control network
"""
function plot_ctl_network_mono(model; kwargs...,)

    ctl_p = Plots.plot(circle_shape(0,0,0.1)
                 , seriestype = [:shape]
                 , lw=0.5
                 , c=:lightblue
                 , linecolor=:black
                 , legend = false
                 #, fillalpha = 0.5
                 , aspect_ratio = 1
                 , showaxis = false
                 ,xlims=[-1,1]
                 ,ylims=[-1,1]
                 ,grid = false
                 ,ticks=false
                 ,titlefontcolor=:white
                )
    annotate!((-1,0.72,Plots.text("Single Controller", 11, :black, :center)))

    return ctl_p

end
"""
 Plots message exchange among control agents
"""
function plot_ctl_throughput(model;kwargs...)
    tpt_v = get_ag_msg(model)
    max_y = maximum(tpt_v)+5
    tpt_p = Plots.plot(title="tpt",titlefontcolor=:white,ylims=[0,max_y])
    # for i=1:nv(model.ctl_graph)
    #     a = getindex(model,model.ctl_graph[i,:aid])
        # tpt_v = get_throughput_up(a,model)

    # log_info("Plotting...")
    # log_info(tpt_v)
    tpt_p = Plots.plot!(tpt_v,xlims=[0,model.N], linealpha=0.5
        # , line=:stem
        ,ylabel = "Quantity of agent messages"
        ,legend = false#:outerright
        )
    # end

    annotate!((3*(model.N/4),max_y+1,Plots.text("Control Msgs", 11, :black, :center)))

    return tpt_p
end
"""
 Plots the underlying controlled network
"""
function plot_asset_networks(model;kwargs...)

    nsize = 0.2
    lwidth = 0.5
    condition_color = cgrad([:red, :yellow, :green],collect(0.00:0.01:1))
    method = model.ntw_model == GraphModel(4) ? :stress : :circular
    Random.seed!(model.seed)

    edge_color_dict = Dict()
    edge_width_dict = Dict()
    edge_style_dict = Dict()
    node_color_dict = Dict()

    for e in edges(model.base_ntw_graph)
        if model.ticks > 0
            if is_active_flow((e.src,e.dst),model)
                edge_color_dict[(e.src,e.dst)] = :green
                edge_width_dict[(e.src,e.dst)] = 3
                edge_style_dict[(e.src,e.dst)] = model.ticks % 3 > 0 ? model.ticks % 3 > 1 ? :dashdot : :solid : :dot
            else
                edge_color_dict[(e.src,e.dst)] = has_edge(model.ntw_graph,e.src,e.dst) || has_edge(model.ntw_graph,e.dst,e.src) ? :dimgray : :white
                edge_width_dict[(e.src,e.dst)] = 1
                edge_style_dict[(e.src,e.dst)] = :solid
            end
        else
            edge_color_dict[(e.src,e.dst)] = :red
            edge_width_dict[(e.src,e.dst)] = 1
            edge_style_dict[(e.src,e.dst)] = :solid
        end

    end

    ruls = [ Int(round(get_state(getindex(model,i)).rul)) for i in 1:nv(model.ntw_graph)]

    # rearrange ruls splitting array in 2 and starting with second half to match plotting algorithms, IDKW.
    ruls = vcat(ruls[nv(model.ntw_graph)÷2+1:nv(model.ntw_graph)],ruls[1:nv(model.ntw_graph)÷2])

    node_colors = [ ruls[i] > 0 ? condition_color[ruls[i]] : :red for i=1:nv(model.ntw_graph) ]
    # log_info(model.ticks," RULs: $(ruls)")
    ntw_p = graphplot(
        model.base_ntw_graph
        ,names = [get_eid(i,model) for i=1:nv(model.ntw_graph)] #[ i for i in ruls ]
        , method = method
        ,size=(300,200)
        ,dpi=400
        ,node_weights = [ get_eid(i,model) > 9 ? 1 : 10 for i in 1:nv(model.ntw_graph)]  #[ i > 9 ? 1 : 10 for i in 1:nv(model.ntw_graph)]
        ,nodeshape = :hexagon
        ,nodecolor = node_colors
        # ,markerstrokecolor = :dimgray
        ,edgecolor= edge_color_dict
        ,edgewidth= edge_width_dict
        ,edgestyle = edge_style_dict
        #,arrow = arrow(:closed, :tail)
        ,markerstrokewidth = 1.1
        ,node_size=nsize
        ,palette = [:lightgray, :red]
        #,titlefontsize=1
        ,titlefontcolor=:white
    )

    # annotate!((-0.7,0.75,Plots.text("Infrastructure Condition & Flows View", 11, :black, :center)))
    annotate!((0.0,0.75,Plots.text("Asset Condition & Flows", 11, :black, :center)))

    return ntw_p
end
"""
 Plots an empty placeholder
"""
function plot_empty()
    return Plots.plot(title="false", titlefontcolor=:white ,showaxis = false, ticks=false,grid=false)
end
"""
 Plots throughput of the underlying controlled network
"""
function plot_throughput(model;kwargs...)
    max_y = 500
    tpt_p = Plots.plot(title="tpt",titlefontcolor=:white,ylims=[0,max_y])
    for i=1:nv(model.ntw_graph)
        sne = getindex(model,get_eid(i,model))
        #tpt_v = get_throughput_up(sne,model)
        tpt_v = get_throughput_trj(sne)
        # log_info(model.ticks,sne.id,"==> tpt_trj: $tpt_v")
        tpt_p = Plots.plot!(tpt_p,tpt_v
        ,xlims=[0,model.N]
        , linealpha=0.5
        # , line=:stem
        ,label = "$i"
        ,ylabel = "MB"
        ,legend = :outerright#false#
        )
    end
    #TODO: This annotation breaks the multithreading as it does not receive the plot object, it seems to take the last one, which might clash among threads.
    # annotate!((3*(model.N/4),max_y+1,Plots.text("Throughput ($(model.interval_tpt) steps)", 11, :black, :center)))
    # annotate!((3*(model.N/4),max_y+5,Plots.text("Throughput", 11, :black, :center)))
    annotate!((2.5*(model.N/4),max_y+15,Plots.text("Throughput", 11, :black, :center)))

    return tpt_p
end
"""
 Plots packet loss of the underlying controlled network
"""
function plot_packet_loss(model;kwargs...)
    max_y = model.pkt_per_tick
    pktl_p = Plots.plot(title="Packet Loss",titlefontcolor=:white,ylims=[0,max_y])
    for i=1:nv(model.ntw_graph)
        sne = getindex(model,get_eid(i,model))
        pktl_v = get_packet_loss_trj(sne)

        pktl_p = Plots.plot!(pktl_p,pktl_v
        ,xlims=[0,model.N]
        , linealpha=0.5
        ,label = "$i"
        ,ylabel = "Q"
        ,legend = :outerright
        )
    end
    #TODO: This annotation breaks the multithreading as it does not receive the plot object, it seems to take the last one, which might clash among threads.
    annotate!((3*(model.N/4),max_y+1,Plots.text("Packet Loss", 11, :black, :center)))

    return pktl_p
end

"""
Main plotting function calls others
"""
function plotabm_networks(model;kwargs...)

    l =  @layout([A{0.01h}; [B C ; D E]])  #(2,2) #@layout [a{1w} [grid(1,2) b{0.2h}] ] #@layout [a{1w};(1,2)]

    title = Plots.plot(title = "Plot title", grid = false, showaxis = false, ticks=false, bottom_margin = -50Plots.px)

    ctl_p = model.ctrl_model != GraphModel(1) ? # Centralised
            plot_ctl_network_multi(model;kwargs...) :
            plot_ctl_network_mono(model;kwargs...)

    ctl_r = model.ctrl_model != GraphModel(1) ? plot_ctl_throughput(model; kwargs) : plot_empty()

    ntw_p = plot_asset_networks(model; kwargs)

    bottom_right_p = plot_throughput(model; kwargs)#plot_packet_loss(model; kwargs) #

    p = Plots.plot(title,ctl_p,ctl_r,ntw_p,bottom_right_p, layout=l
    , size=(600,400),dpi=400)

    #p = ntw_p

    return p
end

"""
  Plots the bare base map using the shp file passed
"""
function plot_base_map(shp_file)
    shp_map_uk_regions = load_map(shp_file)
    #df_shp_map_uk_regions = DataFrame(shp_map_uk_regions)
    points = []

    dpi = 300
    size = (400,500)

    p = plot(shp_map_uk_regions,alpha=0.07,palette=:cork,dpi = dpi,xlims=[0,8e5],size=size,showaxis=false,grid=false,ticks=false,label="")

    return p
end

"""
    Plot network on top of the given geo plot passed
"""
function plot_geo_network(base_geo_plot,df_net)
    dpi = 300
    size = (400,500)
	p = @df df_net scatter!(base_geo_plot,:lon,:lat,markersize=1,alpha=0.5,dpi=dpi,markerstrokewidth=0,label="",color=:red,legend=:topright, markershape=:octagon)
    return p
end

"""
    Plot network on top of the given geo plot passed
"""
function plot_geo_network(base_geo_plot,g::MetaGraph)
    dpi = 300
    size = (400,500)
    n_v = nv(g)
    lons = [ get_prop(g,i,:lon) for i=1:n_v]
    lats = [ get_prop(g,i,:lat) for i=1:n_v]
    p = graphplot!(base_geo_plot,g,x=lons,y=lats,alpha=0.3,dpi=dpi,size=size,nodesize=4.5,aspect_ratio=1,curvature_scalar=0.5,nodecolor="#FDF3C4"
    ,names=collect(1:n_v)
    ,fontsize=3, edgecolor = :green)
    return p
end

"""
It plots one step at a time, receiving all steps in a vector of vector of NetworkAssetState
"""
function plot_geo_network_step(base_geo_plot,g::MetaGraph,sne_steps::Vector{Vector{NetworkAssetState}},step)
    dpi = 300
    size = (400,500)

    condition_color = cgrad([:red, :yellow, :green],collect(0.01:0.01:1))

    lons = [ get_prop(g,i,:lon) for i=1:nv(g)]
    lats = [ get_prop(g,i,:lat) for i=1:nv(g)]
    ruls = [sne_steps[i][step] for i=1:nv(g) ]

    p = graphplot!(base_geo_plot,g,x=lons,y=lats,alpha=0.3,dpi=dpi,size=size,nodesize=2,aspect_ratio=1,curvature_scalar=0.5,nodecolor=[ condition_color[i] for i in ruls ])
    return p
end

"""
It plots one step at a time, receiving all steps in a DataFrame with data coming from NetworkAssetState
"""
function plot_geo_network_step(base_geo_plot,g::MetaGraph,sne_steps::DataFrame,step)
    dpi = 300
    size = (400,500)

    condition_color = cgrad([:red, :yellow, :green],collect(0.01:0.01:1))

    lons = [ get_prop(g,i,:lon) for i=1:nv(g)]
    lats = [ get_prop(g,i,:lat) for i=1:nv(g)]
    ruls = [ r.rul for r in eachrow(sne_steps[sne_steps[!,:t] .== step,:])]./100

    p = graphplot!(base_geo_plot,g,x=lons,y=lats,alpha=0.3,dpi=dpi,size=size,nodesize=2,aspect_ratio=1,curvature_scalar=0.5,nodecolor=[ condition_color[i] for i in ruls ])
    return p
end

function plot_geo_network_step(base_geo_plot,g::MetaGraph,v_props::Dict{Symbol,Symbol},e_props::Dict{Symbol,Symbol},t::Int64)
    dpi = 300
    size = (400,500)

    v_size = 4.5
    v_shape = :hexagon
    n_v = nv(g)
    lons = [ get_prop(g,i,:lon) for i=1:n_v]
    lats = [ get_prop(g,i,:lat) for i=1:n_v]
    e_color,e_width,e_style = get_edge_plot_props_step(g,t,e_props)
    v_color = get_vertex_plot_props_step(g,t,v_props)
    g_plot = graphplot!(
                        base_geo_plot
                        ,g
                        ,x=lons,y=lats
                        ,alpha=0.3
                        ,names = collect(1:n_v)
                        ,node_weights = [ i >9 ? 1 : 10 for i=1:n_v ]
                        ,node_size = v_size
                        ,nodecolor = v_color
                        ,nodeshape = v_shape
                        ,edgecolor = e_color
                        ,edgewidth = e_width
                        ,edgestyle = e_style
                        ,left_margin = -400Plots.px
                        ,right_margin = -250Plots.px
                        ,fontsize=3
                        ,dpi=dpi,size=size, aspect_ratio=1,curvature_scalar=0.5
                        )
    return g_plot
end

"""
Plot network according to arguments
# """
# function plot_network(g::MetaGraph,)
# end

########TODO REVIEW#############

function prepare_graph!(g::MetaGraph,snes_ts::Vector{Vector{NetworkAssetState}},model_ts::Vector{ModelState},props::Dict{Symbol,Symbol})

    rul_prop = props[:rul]
    up_prop = props[:up]
    tpt_prop = props[:tpt]
    drop_prop = props[:drop]

        #vertices
        for nid=1:nv(g)
            sne_ts = snes_ts[nid]
            set_prop!(g,nid,rul_prop,[ s.rul for s in sne_ts ][1:end-1])#remove last to have equal ticks to edges
            set_prop!(g,nid,up_prop,[ s.up for s in sne_ts ][1:end-1])
            set_prop!(g,nid,drop_prop,[ s.drop_pkt for s in sne_ts ][1:end-1])
        end

        #edges
        for e in edges(g)
            if e.src < e.dst
                edge_ts = [ haskey(s.links_load,(e.src,e.dst)) ? s.links_load[(e.src,e.dst)] : haskey(s.links_load,(e.dst,e.src)) ? s.links_load[(e.dst,e.src)] : 0  for s in model_ts ]
                set_prop!(g,e.src,e.dst,tpt_prop, edge_ts)
                set_prop!(g,e.dst,e.src,tpt_prop, edge_ts)
            end

        end
    return g
end

"""
labels is a dictionary where keys are: :tpt for througput and :up for node alive.
"""
function get_edge_plot_props_step(g::MetaGraph,t::Int64,props::Dict{Symbol,Symbol})
    colors = Dict(1=>:blue,14=>:green,7=>:purple,8=>:orange)
    e_color = Dict()
    e_width = Dict()
    e_style = Dict()

    for e in edges(g)
        if get_prop(g,e,props[:tpt])[t] > 0
            e_color[(e.src,e.dst)] = "#2F9D96"#"#1C5E5A" #:green
            e_width[(e.src,e.dst)] = 3
            e_style[(e.src,e.dst)] = t % 3 > 0 ? t % 3 > 1 ? :dashdot : :solid : :dot
        else
            e_width[(e.src,e.dst)] = 1
            e_style[(e.src,e.dst)] = :solid
            if get_prop(g,e.src,props[:up])[t] && get_prop(g,e.dst,props[:up])[t]
                e_color[(e.src,e.dst)] = :gray
            else
                e_color[(e.src,e.dst)] = :white
            end
        end
    end

    return e_color, e_width, e_style
end

function get_vertex_plot_props_step(g::MetaGraph,t::Int64,props::Dict{Symbol,Symbol})
    condition_color = cgrad(["#FF5964", "#FDF3C4", "#669bbc"],collect(0.00:0.01:1))

    ups = [ get_prop(g,v,props[:up])[t] for v=1:nv(g) ]
    ruls = [ get_prop(g,v,props[:rul])[t] for v=1:nv(g) ]
    ruls = vcat(ruls[nv(g)÷2+1:nv(g)],ruls[1:nv(g)÷2])

    v_color = [ ruls[i] > 0 ? condition_color[Int(round(ruls[i]))+1] : :white for i=1:length(ruls)]

    return v_color
end

function plot_graph_step(g::MetaGraph,v_props::Dict{Symbol,Symbol},e_props::Dict{Symbol,Symbol},t::Int64)
    Random.seed!(seed)
    method = :stress
    v_size = 0.4
    v_shape = :hexagon

    e_color,e_width,e_style = get_edge_plot_props_step(g,t,e_props)
    v_color = get_vertex_plot_props_step(g,t,v_props)
    g_plot = graphplot(
                        g
                        ,method = method
                        ,names = collect(1:nv(g))
                        ,node_weights = [ i >9 ? 1 : 10 for i=1:nv(g) ]
                        ,node_size = v_size
                        ,nodecolor = v_color
                        ,nodeshape = v_shape
                        ,edgecolor = e_color
                        ,edgewidth = e_width
                        ,edgestyle = e_style
                        ,bottom_margin = -30Plots.px
                        )
    return g_plot
end

"""
It plots throughput and reference for one step of an animation of t timesteps
"""
function plot_tpt_step(snes_ts::Vector{Vector{NetworkAssetState}},snes_ref_ts::Vector{Vector{NetworkAssetState}},model_ts::Array{ModelState, 1},t::Int64)
    max_x = 300
    colors = Dict(1=>:blue,14=>:green,7=>:purple,8=>:orange)

    #get end node of each active flow (service)
    end_snes_ts = [[1,7,8,14] for tk=1:t]#[ [f[2] for f in model_ts[t].active_flows if f[3] == f_E || f[3] == f_SE] for t=1:t ]
    sort!(end_snes_ts)
    p = plot(xlims=[0,max_x],ylims=[0,250]
    ,xlabel="Time",ylabel="Throughput. (MB)"
    #,xticks=false
    ,legend=false#:outertop#:outerright
    )
    for end_sne in end_snes_ts[t]
        tpts_ref = hcat([ get_throughput_trj(snes_ref_ts[sne],max_x) for sne=1:length(snes_ref_ts) if sne in end_sne ]...)
        p = plot!(tpts_ref, color = colors[end_sne], alpha=0.2
        #,linestyle=:dot
        )
        tpts = hcat([ get_throughput_trj(snes_ts[sne],t) for sne=1:length(snes_ts) if sne == end_sne ]...)

        p = plot!(p,tpts,color=colors[end_sne])
    end
    if t == 80
        print(get_throughput_trj(snes_ts[1],t))
    end
    return p
end

"""
It plots throughput for one step of an animation of t timesteps
"""
function plot_tpt_step(snes_ts::Vector{Vector{NetworkAssetState}},end_services::Vector{Tuple{Int64,Int64}},t::Int64)
    max_x = 300
    colors = Dict(1=>:blue,14=>:green,7=>:purple,8=>:orange)

    #get end node of each active flow (service)
    end_snes_ts = [end_services for tk=1:t]#[ [f[2] for f in model_ts[t].active_flows if f[3] == f_E || f[3] == f_SE] for t=1:t ]
    sort!(end_snes_ts)
    p = plot(xlims=[0,max_x],ylims=[0,250]
    ,xlabel="Time",ylabel="Throughput (MB)"
    #,xticks=false
    ,legend=:outertop#:outerright
    )
    for end_sne in end_snes_ts[t]
        println("doing $end_sne ...")

        # tpts = hcat([ get_throughput_trj(snes_ts[sne],t) for sne=1:length(snes_ts) if sne == last(end_sne) ]...)

        tpts = get_throughput_trj(snes_ts[last(end_sne)],t)

        # if last(end_sne) in [8,14]
            p = plot!(p,tpts,color=colors[last(end_sne)])
        # end
    end
    if t == 80
        print(get_throughput_trj(snes_ts[1],t))
    end
    return p
end

function plot_maintenance_cost_step(snes_ts::Vector{Vector{NetworkAssetState}},model_ts::Array{ModelState, 1},t::Int64)
    costs = zeros(1:t)
    if t > 1

        ruls = [  snes_ts[sne][tk].on_maintenance ? snes_ts[sne][tk].rul : -1.0  for sne=1:length(snes_ts), tk=1:t ]

        is_starts = hcat(Bool.(zeros(1:length(snes_ts))),[  !snes_ts[sne][tk-1].on_maintenance && snes_ts[sne][tk].on_maintenance ? true : false for sne=1:length(snes_ts),tk=2:t ])

        is_actives = transpose(hcat([ is_in.(sne,
                                            [ unique(vcat([ vcat(f[1],f[2]) for f in model_ts[tk].active_flows ]...)) for tk=1:t ]) for sne in collect(1:length(snes_ts)) ]...))

        costs = cumsum(sum(eachrow(
                    maintenance_cost.(
                        ruls,
                        is_starts,
                        is_actives,
                        5, 4, 10, 3
                    )
                )))

        if t == 140

        println("=========RULS=============")
        print(ruls)
        print("\n")
        println("========IS_STARTS==============")
        print(is_starts)
        print("\n")
        println("========IS_ACTIVEs==============")
        print(is_actives)
        print("\n")
        println("==========COSTS============")
        print(costs)
        print("\n")
        end

    end

    return plot(
                costs
                ,legend = false
                ,ylims=[0,1000]
                ,xlims=[0,300]
            )

    # sum_dp = sum(eachrow([  snes_ts[sne][tk].drop_pkt  for sne=1:length(snes_ts), tk=1:t ]))

end

function plot_maintenance_cost_step(snes_ts1::Vector{Vector{NetworkAssetState}},snes_ts2::Vector{Vector{NetworkAssetState}},snes_ts3::Vector{Vector{NetworkAssetState}},model_ts::Vector{ModelState},t::Int64)
    costs = zeros(3,t)

    snes_ts = [snes_ts1,snes_ts2,snes_ts3]

    if t > 1

        for sr=1:length(snes_ts)
            ruls = [  snes_ts[sr][sne][tk].on_maintenance ? snes_ts[sr][sne][tk].rul : -1.0  for sne=1:length(snes_ts[sr]), tk=1:t ]

            is_starts = hcat(Bool.(zeros(1:length(snes_ts[sr]))),[  !snes_ts[sr][sne][tk-1].on_maintenance && snes_ts[sr][sne][tk].on_maintenance ? true : false for sne=1:length(snes_ts[sr]),tk=2:t ])

            is_actives = transpose(hcat([ is_in.(sne,
                                                [ unique(vcat([ vcat(f[1],f[2]) for f in model_ts[tk].active_flows ]...)) for tk=1:t ]) for sne in collect(1:length(snes_ts[sr])) ]...))

            costs[:,sr] = cumsum(sum(eachrow(
                        maintenance_cost.(
                            ruls,
                            is_starts,
                            is_actives,
                            5, 4, 10, 3
                        )
                    )))
        end

    end

    return plot(
                costs
                ,legend = false
                ,ylims=[0,1000]
                ,xlims=[0,300]
            )

    # sum_dp = sum(eachrow([  snes_ts[sne][tk].drop_pkt  for sne=1:length(snes_ts), tk=1:t ]))

end

"""
Function that prepare serialised time series of control agents for plotting
"""
function prepare_service_paths(services,ca_ts)
    t = length(ca_ts)
    paths = [ collect(values(ca_ts[tk].active_paths)) for tk=1:t ]
    path_lengths = [ [ (first(path),last(path),length(path)) for path in paths[tk] ] for tk=1:t ]

    return hcat([ vcat([ [ p[3] for p in path_lengths[tk] if p[1] == first(s) && p[2] == last(s) ] for tk=1:t ]...) for s in services ]...)
end

"""
It loads serialised data for all the runs in the passed directory
"""
function load_run_data(services,data_dir::String)

    run_data = []
    ca_ts = Dict()
    snes_ts = Dict()
    model_ts = Dict()
    root = data_dir
    sdirs = readdir(root)

    for dir in sdirs
        println("processing : $dir ...")
        files = readdir(joinpath(root,dir))
        for file in files
            mdsp = split(file,"_steps_model.bin")
            nesp = split(file,"_steps_nelements.bin")
            agsp =  split(file,"_steps_ctl_agents.bin")
            mdata = length(mdsp) > 1 ? deserialize(joinpath(root, dir,file)) : nothing

            if !isnothing(mdata)
                #key: last of directory + unique run label
                model_ts[first(mdsp)*"_"*last(split(dir,"_"))] = mdata
            end

            snedata = length(nesp) > 1 ? deserialize(joinpath(root, dir,file)) : nothing

            if !isnothing(snedata)
                snes_ts[first(nesp)*"_"*last(split(dir,"_"))] = snedata
            end

            ctlagdata = length(agsp) > 1 ? deserialize(joinpath(root, dir, file)) : nothing

            if !isnothing(ctlagdata)
                ca_ts[first(agsp)*"_"*last(split(dir,"_"))] = ctlagdata
            end
        end
    end

    for k in keys(model_ts)
        push!(run_data,
            (label = k,
            snes_ts=snes_ts[k]
            ,ca_ts=ca_ts[k]
            ,model_ts=model_ts[k]
            ,ctl_model = split(k,"_")[5]
            ,fail_prop= last(split(k,"_"))
            ,seed = length(split(k,"_")) > 11 ? split(k,"_")[10] : split(k,"_")[7]
            )
        )
    end

    return run_data
end

"""
It plots cost for three time series at a time
"""

function plot_maintenance_cost_step(snes_ts1::Vector{Vector{NetworkAssetState}},snes_ts2::Vector{Vector{NetworkAssetState}},snes_ts3::Vector{Vector{NetworkAssetState}},model_ts1::Array{ModelState, 1},model_ts2::Array{ModelState, 1},model_ts3::Array{ModelState, 1},t::Int64)
    costs = []#zeros(3,t)

    snes_ts = [snes_ts1,snes_ts2,snes_ts3]
    model_ts = [model_ts1,model_ts2,model_ts3]

    if t > 1

        for sr=1:length(snes_ts)
            ruls = [  snes_ts[sr][sne][tk].on_maintenance ? snes_ts[sr][sne][tk].rul : -1.0  for sne=1:length(snes_ts[sr]), tk=1:t ]

            is_starts = hcat(Bool.(zeros(1:length(snes_ts[sr]))),[  !snes_ts[sr][sne][tk-1].on_maintenance && snes_ts[sr][sne][tk].on_maintenance ? true : false for sne=1:length(snes_ts[sr]),tk=2:t ])

            is_actives = transpose(hcat([ is_in.(sne,
                                                [ unique(vcat([ vcat(f[1],f[2]) for f in model_ts[sr][tk].active_flows ]...)) for tk=1:t ]) for sne in collect(1:length(snes_ts[sr])) ]...))

            push!(costs,cumsum(sum(eachrow(
                        maintenance_cost.(
                            ruls,
                            is_starts,
                            is_actives,
                            5, 4, 10, 1 #coeficients for dt, l, p & r costs (see function docs)
                        )
                    ))))
        end
        #print(costs[1])
    end

    p = plot(
        title = "Total Costs"
        ,title_location = :right
        ,titlefont = font(8)
        ,legend = :topleft
        #,legend = false
        ,ylims=[0,1500]
        ,xlims=[0,300]
        ,legendfontsize=6
        ,ylabel="£"
        ,xlabel="Time"
        ,guidefontsize=6
        )

    labels = ["Corrective", "Preventive", "Optimal"]
    colors = [:orange,:purple,:green]
     for c=1:length(costs)
         p = plot!(p,costs[c],label=labels[c],color=colors[c])
     end
    return p

    # sum_dp = sum(eachrow([  snes_ts[sne][tk].drop_pkt  for sne=1:length(snes_ts), tk=1:t ]))

end


text:
const μ = 3.986e14 * 60^2
const r = 42164.0e3
const r2 = r^2
const mc = 500.0
const n2 = μ / r^3
const n = sqrt(n2)

const two_n = 2*n
const μ_r2 = μ/r2

const K₁ = [-28.8287 0.1005 -1449.9754 0.0046;
            -0.087 -33.2562 0.00462 -1451.5013]
const K₂ = [-288.0288 0.1312 -9614.9898 0.0;
            -0.1312 -288.0 0.0 -9614.9883]

const K₁mc = K₁/mc
const K₂mc = K₂/mc

function mymul!(v, A, x)
    for ind = 1:length(v)
        v[ind] = zero(x[1])
        for jind = 1:length(x)
            v[ind] += A[ind, jind] * x[jind]
        end
    end
    return nothing
end

# dynamics in the 'approaching' mode
@taylorize function spacecraft_approaching!(t, x, dx)
    x_1 = r + x[1]
    x_12 = x_1^2
    x_22 = x[2]^2
    rc = sqrt(x_12 + x_22)
    rc3 = rc^3
    μ_rc3 = μ / rc3

    uxy = Array{typeof(x[1])}(undef, size(K₁mc,1))
    mymul!(uxy, K₁mc, x[1:4])

    # x' = vx
    dx[1] = x[3]
    # y' = vy
    dx[2] = x[4]
    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    dx[3] = (n2*x[1] + two_n*x[4]) + ((μ_r2 - μ_rc3*x_1) + uxy[1])
    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    dx[4] = (n2*x[2] - two_n*x[3]) - (μ_rc3*x[2] - uxy[2])
    # t' = 1
    dx[5] = one(t)

    return dx
end

# dynamics in the 'rendezvous attempt' mode
@taylorize function spacecraft_rendezvous_attempt!(t, x, dx)
    x_1 = r + x[1]
    x_12 = x_1^2
    x_22 = x[2]^2
    rc = sqrt(x_12 + x_22)
    rc3 = rc^3
    μ_rc3 = μ / rc3

    uxy = Array{typeof(x[1])}(undef, size(K₂mc,1))
    mymul!(uxy, K₂mc, x[1:4])

    # x' = vx
    dx[1] = x[3]
    # y' = vy
    dx[2] = x[4]
    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    dx[3] = (n2*x[1] + two_n*x[4]) + ((μ_r2 - μ_rc3*x_1) + uxy[1])
    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    dx[4] = (n2*x[2] - two_n*x[3]) - (μ_rc3*x[2] - uxy[2])
    # t' = 1
    dx[5] = one(t)

    return dx
end

# dynamics in the 'aborting' mode
@taylorize function spacecraft_aborting!(t, x, dx)
    x_1 = r + x[1]
    x_12 = x_1^2
    x_22 = x[2]^2
    rc = sqrt(x_12 + x_22)
    rc3 = rc^3
    μ_rc3 = μ / rc3

    # x' = vx
    dx[1] = x[3]
    # y' = vy
    dx[2] = x[4]
    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x)
    dx[3] = (n2*x[1] + two_n*x[4]) + (μ_r2 - μ_rc3*x_1)
    # vy' = n²y - 2n*vx - μ/(rc^3)y
    dx[4] = (n^2*x[2] - 2*n*x[3]) - μ/(rc^3)*x[2]
    dx[4] = (n2*x[2] - two_n*x[3]) - μ_rc3*x[2]

    # t' = 1
    dx[5] = one(t)

    return dx
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks