text:
#############
# Adjacency #
#############
abstract type AbstractAdjacency end

# Types
"""
Adjacency type where the adjacency information is represented by an `AbstractMatrix`.
"""
mutable struct MatrixAdjacency{T<:AbstractMatrix} <: AbstractAdjacency
	am::T
	function MatrixAdjacency(am::T) where T<:AbstractMatrix
		@assert issymmetric(am) "Adjacency matrix should be symmeric."
		new{T}(am)
	end
end

"""
Adjacency type where the adjacency information is represented by an `AbstractGraph`.
"""
mutable struct GraphAdjacency{T<:AbstractGraph} <: AbstractAdjacency
	ag::T
end

"""
Adjacency type where the adjacency information is represented by a .
function `f` and `data` on which the function can be applied. The result of
`f(data)` has to be either a matrix or a graph.
"""
mutable struct ComputableAdjacency{T,S} <:AbstractAdjacency
	f::T
	data::S
end

"""
Adjacency type where the adjacency information is represented by a
function `f` which can be used to calculate an adjacency matrix or graph,
given proper data."
"""
mutable struct PartialAdjacency{T<:Function} <: AbstractAdjacency
	f::T
end

"Adjacency type indicating the lack of any adjacency information."
mutable struct EmptyAdjacency <: AbstractAdjacency end

# Constructors
"""
	adjacency([a,data])

Constructs an adjacency object. If `a` is a `AbstractMatrix`, `AbstractGraph` or
`Tuple`, it will return usable adjacencies. If `a` is a `Function`
or `PartialAdjacency`, `data` has to be present. If both arguments are missing,
the function returns an `EmptyAdjacency`.

# Examples
```
julia> using NetworkLearning, LightGraphs

julia> A = [0 1 0; 1 0 0; 0 0 0];

julia> Am = adjacency(A)
Matrix adjacency, 3 obs

julia> Ag = adjacency(Graph(A))
Graph adjacency, 3 obs

julia> Ac = adjacency(x->x,A)
Computable adjacency
```
"""
adjacency(a::T) where T<:AbstractAdjacency = a
adjacency(a::T, data) where T<:PartialAdjacency = adjacency(a.f(data))
adjacency(am::T) where T<:AbstractMatrix = MatrixAdjacency(am)
adjacency(ag::T) where T<:AbstractGraph = GraphAdjacency(ag)
adjacency(f::T, data::S) where {T,S} = ComputableAdjacency(f,data)
adjacency(t::Tuple{T,S}) where {T,S} = ComputableAdjacency(t[1],t[2])
adjacency(f::T) where T<:Function = PartialAdjacency(f)
adjacency(::Nothing) = EmptyAdjacency()
adjacency() = EmptyAdjacency()

# Show methods
Base.show(io::IO, a::MatrixAdjacency) = print(io, "Matrix adjacency, $(size(a.am,1)) obs")
Base.show(io::IO, a::GraphAdjacency) = print(io, "Graph adjacency, $(nv(a.ag)) obs")
Base.show(io::IO, a::ComputableAdjacency) = print(io, "Computable adjacency")
Base.show(io::IO, a::PartialAdjacency) = print(io, "Partial adjacency, not computable")
Base.show(io::IO, a::EmptyAdjacency) = print(io, "Empty adjacency, not computable")
Base.show(io::IO, va::T) where T<:AbstractVector{S} where S<:AbstractAdjacency =
	print(io, "$(length(va))-element Vector{$S} ...")

# Methods that to strip the adjacency infomation (used in training of the NetworkLearner); the methods return a Partial Adjacency
"""
	strip_adjacency(a)

Function that removes adjancency information (i.e. matrix, graph or other data) from adjacency objects
and returns a `PartialAdjacency` that can be used in conjunction with the `adjacency` function
to build a new adjacency object.

# Examples
```
julia> using NetworkLearning, LightGraphs

julia> A = [0 1 0; 1 0 0; 0 0 0];

julia> Am = adjacency(A)
Matrix adjacency, 3 obs

julia> Ac = adjacency(x->x,A)
Computable adjacency

julia> Sm = strip_adjacency(Am)
Partial adjacency, not computable

julia> adjacency(Sm, A) # Sm has to be used with a matrix
Matrix adjacency, 3 obs

julia> Sc = strip_adjacency(Ac)
Partial adjacency, not computable

julia> adjacency(Sc, A) # Sc can be use with any data; calls f(A) where f = x->x
Matrix adjacency, 3 obs
```
"""
strip_adjacency(a::MatrixAdjacency{T}) where T <:AbstractMatrix = adjacency(x->T(x))
strip_adjacency(a::GraphAdjacency{T}) where T<:AbstractGraph = adjacency(x->T(x))
strip_adjacency(a::ComputableAdjacency{T,S}) where {T,S} = adjacency(a.f)
strip_adjacency(a::PartialAdjacency) = a
strip_adjacency(a::EmptyAdjacency) = PartialAdjacency(x->adjacency(x))

# Functions to strip the adjacency infomation (used in training of the NetworkLearner); the methods return a Partial Adjacency
"""
	update_adjacency!(a,f_update)

Function that updates the data of an adjacency object. `a` has to be a `MatrixAdjacency` or
`GraphAdjacency` while `f_update` has to of the form `f_update = x->update_function!(x)`.

# Examples
julia> using NetworkLearning, LightGraphs

julia> A = [0 1 0; 1 0 0; 0 0 0];

julia> Am = adjacency(A)
Matrix adjacency, 3 obs

julia> update_function!(X,x,y) = begin
         X[x,y] += 1
         X[y,x] += 1
         return X
       end
update_function! (generic function with 2 methods)

julia> f_update(x,y) = X->update_function!(X,x,y)
f_update (generic function with 1 method)

julia> for i in 1:3
         update_adjacency!(Am, f_update(1,3)) # call function three times
       end

julia> adjacency_matrix(Am)
3×3 Array{Int64,2}:
 0  1  3
 1  0  0
 3  0  0
"""
update_adjacency!(a::MatrixAdjacency{T}, f_update) where T <:AbstractMatrix = begin
	sam = size(a.am)
	f_update(a.am) # f_update should be a closure calling an in-place modifier i.e.  x->foo!(x, args...)
	@assert size(a.am) == sam "The dimensions of the adjacency matrix cannot change."
	return a
end

update_adjacency!(a::GraphAdjacency{T}, f_update) where T<:AbstractGraph = begin
	n = nv(a.ag)
	f_update(a.ag) # f_update should be a closure calling an in-place modifier i.e.  x->foo!(x, args...)
	@assert nv(a.ag) == n "The number of vertices in the adjacency graph cannot change."
	return a
end

update_adjacency!(a::ComputableAdjacency{T,S}, f_update) where {T,S} = error("Cannot update a computable adjacency.")
update_adjacency!(a::PartialAdjacency, args...) = error("Cannot update a partial adjacency.")
update_adjacency!(a::EmptyAdjacency, args...) = error("Cannot update an empty adjacency.")

# Return an adjacency graph from Adjacency types
"""
	adjacency_graph(a)

Returns an adjacency graph computed from the adjacency information of `a`.
"""
adjacency_graph(a::MatrixAdjacency{T}) where T <:AbstractMatrix = adjacency_graph(a.am)
adjacency_graph(a::GraphAdjacency{T}) where T = adjacency_graph(a.ag)
adjacency_graph(a::ComputableAdjacency{T,S}) where {T,S} = adjacency_graph(a.f, a.data)
adjacency_graph(a::PartialAdjacency) = error("Insufficient information to obtain an adjacency graph.")
adjacency_graph(a::EmptyAdjacency) = error("Insufficient information to obtain an adjacency graph.")
adjacency_graph(am::T) where T <:AbstractMatrix = Graph(am)
adjacency_graph(am::T) where T <:AbstractMatrix{Float64} = SimpleWeightedGraph(am)
adjacency_graph(ag::T) where T <:AbstractGraph = ag
adjacency_graph(f::T, data::S) where {T,S} = adjacency_graph(f(data))

# Return an adjacency matrix from Adjacency types
"""
	adjacency_matrix(a)

Returns an adjacency matrix computed from the adjacency information of `a`.
"""
adjacency_matrix(a::MatrixAdjacency{T}) where T <:AbstractMatrix = adjacency_matrix(a.am)
adjacency_matrix(a::GraphAdjacency{T}) where T = adjacency_matrix(a.ag)
adjacency_matrix(a::ComputableAdjacency{T,S}) where {T,S} = adjacency_matrix(a.f, a.data)
adjacency_matrix(a::PartialAdjacency) = error("Insufficient information to obtain an adjacency matrix.")
adjacency_matrix(a::EmptyAdjacency) = error("Insufficient information to obtain an adjacency matrix.")
adjacency_matrix(am::T) where T <:AbstractMatrix = am
adjacency_matrix(ag::T) where T <:AbstractGraph = sparse(ag)
adjacency_matrix(ag::T) where T <:AbstractSimpleWeightedGraph = weights(ag)
adjacency_matrix(f::T, data::S) where {T,S} = adjacency_matrix(f(data))

"""
	adjacency_obs(A, r, obsdim)

Selects a range `r::UnitRange` of observations from the matrix `A`, along the dimension `obsdim::LearnBase.ObsDimension`.
"""
adjacency_obs(A::T where T<:AbstractMatrix, r::S where S<:UnitRange, ::ObsDim.Constant{2}) = view(A,:,r)
adjacency_obs(A::T where T<:AbstractMatrix, r::S where S<:UnitRange, ::ObsDim.Constant{1}) = view(A,r,:)
adjacency_obs(A::T where T<:SparseMatrixCSC, r::S where S<:UnitRange, ::ObsDim.Constant{2}) = A[:,r]
adjacency_obs(A::T where T<:SparseMatrixCSC, r::S where S<:UnitRange, ::ObsDim.Constant{1}) = A[r,:]


text:
using RobustAndOptimalControl
G1 = ssrand(2,3,4, proper=true)
G2 = ssrand(2,3,4, proper=true)

@test h2norm(G1-G1) < 1e-6
@test hinfnorm2(G1-G1)[1] < 1e-6
@test nugap(G1, G1)[1] < 1e-6
@test hankelnorm(G1-G1)[1] < 1e-6
@test hinfnorm2(baltrunc2(G1, n=4)[1]-G1)[1] < 1e-10

@test count(1:100) do _
    G1 = ssrand(1,1,3, proper=true)
    G2 = ssrand(1,1,3, proper=true)
    try
        νgap(G1\G2, ss(inv(tf(G1))*tf(G2)))[1] < 1e-8
    catch
        false
    end
end >= 95

ext: jl


text:
# crisis.jl
#######################

# Load packages.
using Plots, CSV, Dates, DataFrames

# Download data from the MHLW web site.
download("https://covid19.mhlw.go.jp/public/opendata/confirmed_cases_cumulative_daily.csv","./csv/mhlw_cases.csv");
download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv","./csv/jhu_cases.csv");
download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv","./csv/jhu_deaths.csv");

# Shape up the data.
Ccsv=DataFrame(CSV.File("./csv/mhlw_cases.csv", header=false, delim=','));
Wcsv=DataFrame(CSV.File("./csv/jhu_cases.csv", header=false, delim=','));
Xcsv=DataFrame(CSV.File("./csv/jhu_deaths.csv", header=false, delim=','));
(pa,qa)=size(Ccsv);
l=Int64(floor((pa-1)/48));
(pw,qw)=size(Wcsv);

# Plot the data
# d0: the initial date
# df: the final day
dd0=Date(2021,7,1)
DD=Int64(qw-530);
dd1=dd0+Day(floor((DD-1)/2));
dd2=dd0+Day(DD-1);
ll0=string(dd0);
ll1=string(dd1);
ll2=string(dd2);

L=Ccsv[2:pa,2];
J=parse.(Float64,Ccsv[2:pa,3]);

# Okinawa
POKNW=1.458870;
ROWOKNW=findall(x->x=="Okinawa",L);
COKNW=J[ROWOKNW]/POKNW
NOKNW=zeros(DD);
for j=1:DD
    NOKNW[j]=(COKNW[j+418]-COKNW[j+411])/7
end

# Tokyo
PTKY=14.049146;
ROWTKY=findall(x->x=="Tokyo",L);
CTKY=J[ROWTKY]/PTKY
NTKY=zeros(DD);
for j=1:DD
    NTKY[j]=(CTKY[j+418]-CTKY[j+411])/7
end

# Argentina
PARG=45.670451;
BARG=parse.(Float64,Array(Xcsv[8,5:qw]))/PARG;
NDARG=zeros(DD);
for j=1:DD
    NDARG[j]=(BARG[j+526]-BARG[j+519])/7
end

# NSW
PNSW=8.196;
ANSW=parse.(Float64,Array(Wcsv[11,5:qw]))/PNSW;
NNSW=zeros(DD);
for j=1:DD
    NNSW[j]=(ANSW[526+j]-ANSW[519+j])/7
end

# Victoria
PVIC=6.7;
AVIC=parse.(Float64,Array(Wcsv[16,5:qw]))/PVIC;
NVIC=zeros(DD);
for j=1:DD
    NVIC[j]=(AVIC[526+j]-AVIC[519+j])/7
end

# Brazil
PBRA=214.289417;
BBRA=parse.(Float64,Array(Xcsv[32,5:qw]))/PBRA;
NDBRA=zeros(DD);
for j=1:DD
    NDBRA[j]=(BBRA[j+526]-BBRA[j+519])/7
end

# Brunei Darussalam
PBWN=0.442205;
ABWN=parse.(Float64,Array(Wcsv[33,5:qw]))/PBWN;
NBWN=zeros(DD);
for j=1:DD
    NBWN[j]=(ABWN[526+j]-ABWN[519+j])/7
end

# Colombia
PCOL=51.503463;
BCOL=parse.(Float64,Array(Xcsv[94,5:qw]))/PCOL;
NDCOL=zeros(DD);
for j=1:DD
    NDCOL[j]=(BCOL[j+526]-BCOL[j+519])/7
end

# Indonesia
PIDN=276.833206;
AIDN=parse.(Float64,Array(Wcsv[150,5:qw]))/PIDN;
NIDN=zeros(DD);
for j=1:DD
    NIDN[j]=(AIDN[526+j]-AIDN[519+j])/7
end
BIDN=parse.(Float64,Array(Xcsv[150,5:qw]))/PIDN;
NDIDN=zeros(DD);
for j=1:DD
    NDIDN[j]=(BIDN[j+526]-BIDN[j+519])/7
end

# Israel
PISR=8.789774;
AISR=parse.(Float64,Array(Wcsv[154,5:qw]))/PISR;
NISR=zeros(DD);
for j=1:DD
    NISR[j]=(AISR[526+j]-AISR[519+j])/7
end
BISR=parse.(Float64,Array(Xcsv[154,5:qw]))/PISR;
NDISR=zeros(DD);
for j=1:DD
    NDISR[j]=(BISR[j+526]-BISR[j+519])/7
end

# Japan
PJPN=125.36;
AJPN=parse.(Float64,Array(Wcsv[157,5:qw]))/PJPN;
NJPN=zeros(DD);
for j=1:DD
    NJPN[j]=(AJPN[526+j]-AJPN[519+j])/7
end

# South Korea
PKOR=51.318552;
AKOR=parse.(Float64,Array(Wcsv[162,5:qw]))/PKOR;
NKOR=zeros(DD);
for j=1:DD
    NKOR[j]=(AKOR[526+j]-AKOR[519+j])/7
end

# Malaysia
PMYS=32.66;
AMYS=parse.(Float64,Array(Wcsv[178,5:qw]))/PMYS;
NMYS=zeros(DD);
for j=1:DD
    NMYS[j]=(AMYS[526+j]-AMYS[519+j])/7
end

# Mexico
PMEX=130.482814;
BMEX=parse.(Float64,Array(Xcsv[185,5:qw]))/PMEX;
NDMEX=zeros(DD);
for j=1:DD
    NDMEX[j]=(BMEX[j+526]-BMEX[j+519])/7
end

# Peru
PPER=33.510887;
APER=parse.(Float64,Array(Wcsv[213,5:qw]))/PPER;
NPER=zeros(DD);
for j=1:DD
    NPER[j]=(APER[526+j]-APER[519+j])/7
end
BMYS=parse.(Float64,Array(Xcsv[178,5:qw]))/PMYS;
NDMYS=zeros(DD);
for j=1:DD
    NDMYS[j]=(BMYS[j+526]-BMYS[j+519])/7
end
BPER=parse.(Float64,Array(Xcsv[213,5:qw]))/PPER;
NDPER=zeros(DD);
for j=1:DD
    NDPER[j]=(BPER[j+526]-BPER[j+519])/7
end

# Philippines
PPHI=111.249116;
APHI=parse.(Float64,Array(Wcsv[214,5:qw]))/PPHI;
NPHI=zeros(DD);
for j=1:DD
    NPHI[j]=(APHI[526+j]-APHI[519+j])/7
end
BMYS=parse.(Float64,Array(Xcsv[178,5:qw]))/PMYS;
NDMYS=zeros(DD);
for j=1:DD
    NDMYS[j]=(BMYS[j+526]-BMYS[j+519])/7
end
BPHI=parse.(Float64,Array(Xcsv[214,5:qw]))/PPHI;
NDPHI=zeros(DD);
for j=1:DD
    NDPHI[j]=(BPHI[j+526]-BPHI[j+519])/7
end

# Singapore
PSIN=5.902011;
ASIN=parse.(Float64,Array(Wcsv[232,5:qw]))/PSIN;
NSIN=zeros(DD);
for j=1:DD
    NSIN[j]=(ASIN[526+j]-ASIN[519+j])/7
end

# Sri Lanka
PLKA=21.516097;
ALKA=parse.(Float64,Array(Wcsv[240,5:qw]))/PLKA;
NLKA=zeros(DD);
for j=1:DD
    NLKA[j]=(ALKA[526+j]-ALKA[519+j])/7
end
BLKA=parse.(Float64,Array(Xcsv[240,5:qw]))/PLKA;
NDLKA=zeros(DD);
for j=1:DD
    NDLKA[j]=(BLKA[j+526]-BLKA[j+519])/7
end

# Thailand
PTHA=70.000662;
ATHA=parse.(Float64,Array(Wcsv[250,5:qw]))/PTHA;
NTHA=zeros(DD);
for j=1:DD
    NTHA[j]=(ATHA[526+j]-ATHA[519+j])/7
end
BTHA=parse.(Float64,Array(Xcsv[250,5:qw]))/PTHA;
NDTHA=zeros(DD);
for j=1:DD
    NDTHA[j]=(BTHA[j+526]-BTHA[j+519])/7
end

# United Staes
PUSA=333.225477;
AUSA=parse.(Float64,Array(Wcsv[256,5:qw]))/PUSA;
NUSA=zeros(DD);
for j=1:DD
    NUSA[j]=(AUSA[526+j]-AUSA[519+j])/7
end
BUSA=parse.(Float64,Array(Xcsv[256,5:qw]))/PUSA;
NDUSA=zeros(DD);
for j=1:DD
    NDUSA[j]=(BUSA[j+526]-BUSA[j+519])/7
end

# United Kingdom
PGBR=68.294438;
AGBR=parse.(Float64,Array(Wcsv[271,5:qw]))/PGBR;
NGBR=zeros(DD);
for j=1:DD
    NGBR[j]=(max(AGBR[526+j]-AGBR[519+j],0))/7
end

# Vietnam
PVNM=98.341025;
AVNM=parse.(Float64,Array(Wcsv[276,5:qw]))/PVNM;
NVNM=zeros(DD);
for j=1:DD
    NVNM[j]=(AVNM[526+j]-AVNM[519+j])/7
end

p=plot([NUSA NPHI NOKNW NMYS NTHA NBWN NNSW NLKA NISR NGBR],
    grid=false,
    linewidth=2,
    title="COVID-19 7-day average of daily new cases per 1M \n data sourced by JHU and MOH of Japan",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    xticks = ([1 floor(DD/2) DD;], [ll0 ll1 ll2]),
    xlabel="date",
    yaxis="cases/1M",
    legendfont=font(8),
    label=["United States" "Philippines" "Okinawa" "Malaysia" "Thailand" "Brunei" "New South Wales" "Sri Lanka" "Israel" "United Kingdom"],
    palette = :seaborn_bright,
    legend = :topleft)
savefig("./crisis/cases.png")

q=plot([NJPN NTKY NSIN NKOR NVNM NVIC NIDN],
    grid=false,
    linewidth=2,
    title="COVID-19 7-day average of daily new cases per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    xticks = ([1 floor(DD/2) DD;], [ll0 ll1 ll2]),
    xlabel="date",
    yaxis="cases/1M",
    legendfont=font(8),
    label=["Japan" "Tokyo" "Singapore" "South Korea" "Vietnam" "Victoria" "Indonesia"],
    palette = :seaborn_bright,
    legend = :topleft)
savefig("./crisis/delta.png")

r=plot([NDARG NDBRA NDIDN NDMYS NDMEX NDPER NDLKA NDTHA NDUSA NDISR],
    grid=false,
    linewidth=2,
    title="COVID-19: 7-day average deaths per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(DD/2) DD;], [ll0 ll1 ll2]),
    xlabel="date",
    yaxis="deaths/1M",
    legendfont=font(8),
    label=["Argentina" "Brazil" "Indonesia" "Malaysia" "Mexico" "Peru" "Sri Lanka" "Thailand" "United States" "Israel"],
   palette = :seaborn_bright,
     legend = :topleft)
savefig("./crisis/deaths.png")

# Download data from the MHLW web site.
download("https://raw.githubusercontent.com/govex/COVID-19/master/data_tables/vaccine_data/global_data/time_series_covid19_vaccine_global.csv","./csv/jhu_vaccinated.csv");

Vcsv=DataFrame(CSV.File("./csv/jhu_vaccinated.csv", header=false, delim=','));
(m,n)=size(Vcsv);
ROW=Vcsv[2:m,1];
UID=Vcsv[2:m,7];
VAC=Vcsv[2:m,5];

# Plot the data
# d0: the initial date
# df: the final day
d0=Date(2021,7,1)
d2=Date(Vcsv[m,2]);
diff = d2-d0;
D = 1+Dates.value(diff)
d1=d0+Day(floor((D-1)/2));

l0=string(d0);
l1=string(d1);
l2=string(d2);

# Okinawa
POKNW=1.458870;

# Tokyo
PTKY=14.049146;

# Japan
PJPN=1253600
ROWJPN=findall(x->x=="Japan",ROW);
VACJPN=VAC[ROWJPN];
(DJPN,)=size(VACJPN);
FVACJPN=zeros(DJPN);
for j=1:DJPN
    FVACJPN[j]=parse(Float64,VACJPN[j]);
end
VJPN=FVACJPN[DJPN-D+1:DJPN]/PJPN;

# NSW
#PNSW=8.196;

# Brunei Darussalam
PBWN=0.442205;

# Indonesia
PIDN=2768332.06;
ROWIDN=findall(x->x=="Indonesia",ROW);
PVACIDN=VAC[ROWIDN];
PUIDIDN=UID[ROWIDN];
UIDIDN=findall(x->x=="360",PUIDIDN);
VACIDN=PVACIDN[UIDIDN];
(DIDN,)=size(VACIDN);
FVACIDN=zeros(DIDN);
for j=1:DIDN
    FVACIDN[j]=parse(Float64,VACIDN[j]);
end
VIDN=FVACIDN[DIDN-D+1:DIDN]/PIDN;

# South Korea
#PKOR=51.318552;

# Malaysia
PMYS=326600;
ROWMYS=findall(x->x=="Malaysia",ROW);
PVACMYS=VAC[ROWMYS];
PUIDMYS=UID[ROWMYS];
UIDMYS=findall(x->x=="458",PUIDMYS);
VACMYS=PVACMYS[UIDMYS];
(DMYS,)=size(VACMYS);
FVACMYS=zeros(DMYS);
for j=1:DMYS
    FVACMYS[j]=parse(Float64,VACMYS[j]);
end
VMYS=FVACMYS[DMYS-D+1:DMYS]/PMYS;

# Philippines
PPHI=1112491.16;
ROWPHI=findall(x->x=="Philippines",ROW);
VACPHI=VAC[ROWPHI];
(DPHI,)=size(VACPHI);
FVACPHI=zeros(DPHI);
for j=1:DPHI
    FVACPHI[j]=parse(Float64,VACPHI[j]);
end
VPHI=FVACPHI[DPHI-D+1:DPHI]/PPHI;

# Thailand
PTHA=700006.62;
ROWTHA=findall(x->x=="Thailand",ROW);
VACTHA=VAC[ROWTHA];
(DTHA,)=size(VACTHA);
FVACTHA=zeros(DTHA);
for j=1:DTHA
    FVACTHA[j]=parse(Float64,VACTHA[j]);
end
VTHA=FVACTHA[DTHA-D+1:DTHA]/PTHA;

# NSW
#PNSW=8.196;

# Brunei Darussalam
PBWN=4422.05;
ROWBWN=findall(x->x=="Brunei",ROW);
VACBWN=VAC[ROWBWN];
(DBWN,)=size(VACBWN);
FVACBWN=zeros(DBWN);
for j=1:DBWN
    FVACBWN[j]=parse(Float64,VACBWN[j]);
end
VBWN=FVACBWN[DBWN-D+1:DBWN]/PBWN;

# Singapore
PSIN=59020.11;
ROWSIN=findall(x->x=="Singapore",ROW);
VACSIN=VAC[ROWSIN];
(DSIN,)=size(VACSIN);
FVACSIN=zeros(DSIN);
for j=1:DSIN
    FVACSIN[j]=parse(Float64,VACSIN[j]);
end
VSIN=FVACSIN[DSIN-D+1:DSIN]/PSIN;

# Vietnam
#PVNM=983410.25;
PVNM=59020.11;
ROWVNM=findall(x->x=="Vietnam",ROW);
VACVNM=VAC[ROWVNM];
(DVNM,)=size(VACVNM);
FVACVNM=zeros(DVNM);
for j=1:DVNM
    FVACVNM[j]=parse(Float64,VACVNM[j]);
end
VVNM=FVACVNM[DVNM-D+1:DVNM]/PVNM;

# Sri Lanka
#PLKA=21.516097;
# United Staes
#PUSA=333.225477;
# United Kingdom
#PGBR=68.294438;

s=plot([VJPN VPHI VBWN VIDN VMYS VSIN VTHA VVNM],
    grid=false,
    linewidth=2,
    title="COVID-19 fully vaccinated \n data sourced by JHU",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    xticks = ([1 floor(D/2) D;], [l0 l1 l2]),
    xlabel="date",
    yaxis="%",
    legendfont=font(8),
    label=["Japan" "Philippines" "Brunei Darussalam" "Indonesia" "Malaysia" "Singapore" "Thailand" "Vietnam"],
    palette = :seaborn_bright,
    legend = :topleft)
savefig("./crisis/vaccination_asia.png")

plot(p, q, s, r,
     layout=4,
     size=(1260,840),
     left_margin=Plots.Measures.Length(:mm, 5.0),
     right_margin=Plots.Measures.Length(:mm, 15.0),
     top_margin=Plots.Measures.Length(:mm, 5.0),
     bottom_margin=Plots.Measures.Length(:mm, 5.0))
savefig("./crisis/crisis.png")


text:
#=

             Positive X
      -------------------------------------->
      |
p     |
o     |
s     |
i     |
t     |
i     |
v     |
e     |
      |
y     |
      V

=#
function reset_ball!(arena::Arena, ball::Ball)
    k = -1 * sign(ball.velocity.x)
    ball_vx = 2k
    ball_vy = 0
    while true
        ball_vy = -1 * k * rand(-2:2)
        ball_vy != 0 && break
    end
    ball.velocity = Vector_(ball_vx, ball_vy)
    return arena.dims/2f0
end

human_action(player_b::Player, arena::Arena, dir) = update_paddle_position(player_b, arena, dir)

function ai_action(player_a::Player, arena::Arena, ball::Ball)
    dir = sign(ball.position.y - player_a.length/2f0 - player_a.position.y)

    if sign(ball.position.x - player_a.position.x) == sign(ball.velocity.x) ||
        (player_a.position.y ≤ ball.position.y ≤ player_a.position.y + player_a.length)
        dir = 0
    end

    update_paddle_position(player_a, arena, dir)
end

"""The default game configuration"""
function Env()
    arena_width = 80f0
    arena_height = 80f0
    paddle_length = arena_height/5

    margin_height = floor(arena_height/30f0)
    margin_width = floor(arena_width/30f0)

    ball_radius = floor((margin_width + margin_height)/3f0)

    arena = Arena(Vector_(arena_width, arena_height), Vector_(margin_width, margin_height))

    # The position of the paddle is calculated from its top edge.
    # Therefore, acceptable values of the y coordinate of the paddles are
    # from 0 -------> arena_height - paddle_length
    player_a_pos = Vector_(0, (arena.dims.y-paddle_length)/2)
    player_a = Player(paddle_length, deepcopy(player_a_pos), 4, 0)

    player_b_pos = Vector_(arena.dims.x, (arena.dims.y-paddle_length)/2)
    player_b = Player(paddle_length, deepcopy(player_b_pos), 2, 0)

    # Position of the ball is calculated from its center.
    # Acceptable values of the x and y coordinates for the ball are from
    # ball.radius -----------> arena_width - ball.radius
    # ball.radius -----------> arena_height - ball.radius
    ball_pos = Vector_(arena.dims.x/2f0, ball_radius)
    ball_vel = Vector_(3,  -3)
    ball = Ball(deepcopy(ball_pos), deepcopy(ball_vel), ball_radius)

    return Env(arena, player_a, player_b, ball)
end

function step!(env::Env, player_action)
    #@assert player_action ∈ action_space

    # 1 = Noop -1
    # 2 = Up    0
    # 3 = Down  1

    player_action -= 2

    arena = env.arena
    player_a = env.player_a
    player_b = env.player_b
    ball = env.ball

    new_player_b_position = human_action(player_b, arena, player_action)  # Find gradients here (Done)
    new_player_a_position = ai_action(player_a, arena, ball)
    new_ball_position = ball.position + ball.velocity
    new_ball_position = wall_collision!(new_ball_position, ball, arena)

    new_ball_position = paddle_collision!(new_ball_position, new_player_a_position, new_player_b_position,
                                            player_b.length, ball, arena)  # Find gradients here.

    reward = 0

    if new_ball_position.x ≤ 0
        player_b.score += 1
        reward += 80
        new_ball_position = reset_ball!(arena, ball)  # Drop gradients here...
    end

    if new_ball_position.x ≥ arena.dims.x
        player_a.score += 1
        reward -= 80
        new_ball_position = reset_ball!(arena, ball)  # Drop gradients here...
    end

    done = (player_a.score + player_b.score) ≥ 21
    reward += -1 * abs(new_ball_position.y - new_player_b_position.y - player_b.length/2f0)

    # Update positions...
    ball.position = new_ball_position
    player_b.position = new_player_b_position
    player_a.position = new_player_a_position

    return get_obs(env), reward, done, Dict()
end

function reset!(env::Env)
    reset_ball!(env.arena, env.ball)
    env.player_a.score = 0
    env.player_b.score = 0
    get_obs(env)
end


text:
using Oceananigans.Utils: work_layout
using Oceananigans.Architectures: device
using Oceananigans.TimeSteppers: rk3_substep_field!

import Oceananigans.TimeSteppers: rk3_substep!

function rk3_substep!(model::ShallowWaterModel, Δt, γⁿ, ζⁿ)

    workgroup, worksize = work_layout(model.grid, :xyz)

    barrier = Event(device(model.architecture))

    substep_solution_kernel! = rk3_substep_solution!(device(model.architecture), workgroup, worksize)
    substep_tracer_kernel! = rk3_substep_tracer!(device(model.architecture), workgroup, worksize)

    solution_event = substep_solution_kernel!(model.solution,
                                              Δt, γⁿ, ζⁿ,
                                              model.timestepper.Gⁿ,
                                              model.timestepper.G⁻;
                                              dependencies=barrier)

    events = [solution_event]

    for i in 1:length(model.tracers)
        @inbounds c = model.tracers[i]
        @inbounds Gcⁿ = model.timestepper.Gⁿ[i+3]
        @inbounds Gc⁻ = model.timestepper.G⁻[i+3]

        tracer_event = substep_tracer_kernel!(c, Δt, γⁿ, ζⁿ, Gcⁿ, Gc⁻, dependencies=barrier)

        push!(events, tracer_event)
    end

    wait(device(model.architecture), MultiEvent(Tuple(events)))

    return nothing
end

"""
Time step solution fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_solution!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U.uh[i, j, k] += Δt * (γⁿ * Gⁿ.uh[i, j, k] + ζⁿ * G⁻.uh[i, j, k])
        U.vh[i, j, k] += Δt * (γⁿ * Gⁿ.vh[i, j, k] + ζⁿ * G⁻.vh[i, j, k])
        U.h[i, j, k]  += Δt * (γⁿ * Gⁿ.h[i, j, k]  + ζⁿ * G⁻.h[i, j, k])
    end
end

"""
Time step solution fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_solution!(U, Δt, γ¹, ::Nothing, G¹, G⁰)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U.uh[i, j, k] += Δt * γ¹ * G¹.uh[i, j, k]
        U.vh[i, j, k] += Δt * γ¹ * G¹.vh[i, j, k]
        U.h[i, j, k]  += Δt * γ¹ * G¹.h[i, j, k]
    end
end


text:
function readlibsvm(filename::ASCIIString, shape)
  #=
  Read in an libsvm formatted ascii file.
  =#
  data = zeros(Float32, shape)
  label = Float32[]

  f = open(filename, "r")

  count = 1
  for line in eachline(f)
    line = split(line, " ")
    push!(label, float(line[1]))
    line = line[2:end]
    for item in line
      itm = split(itm, ":")
      data[count, int(item[1])+1] = float(int(itm[2]))
    end
    count += 1
  end
  close(f)

  return (data, label)
end

max_stars_repo_name: potterzot/jstats

max_issues_repo_name: potterzot/jstats

max_forks_repo_name: potterzot/jstats


text:
#using GeneralizedGenerated
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

include("expression.jl")

#include("coordinates.jl")
#using .Coordinatesm

let next_id::Int64 = 1
    function randomize!(params::Array{Float64,1},delta=0.0)
        for index in eachindex(params)
            params[index] = params[index] * rand( (1-delta):0.01:(1+delta))
        end
    end
    mutable struct Organism
        id::Int64
        energy::Int64
        birth_time::Float64
        is_sibling::Bool
        generation::Int64

        trigger_split_expression::Expr
        trigger_split_expression_parameter::Array{Float64,1}
        trigger_split_function

        housekeeping_expression::Expr
        housekeeping_expression_parameter::Array{Float64,1}
        housekeeping_function

        energy_expression::Expr
        energy_expression_parameter::Array{Float64,1}
        energy_function

        ϕ_expression::Expr
        ϕ_expression_parameter::Array{Float64,1}
        #ϕ_function::Function
        ϕ_function
        θ_expression::Expr
        θ_expression_parameter::Array{Float64,1}
        #θ_function::Function
        θ_function
        current_position::Coordinates

        target_position::Coordinates
        external_time::Float64

        function Organism(posCoord::PossibleCoordinates, birth_time::Float64,org_max_energie::Int)
            id=next_id
            next_id+=1
            energy=rand(10:org_max_energie)
            is_sibling=false
            generation=0

            (trigger_split_expression,pcount)=create_expression( :life_time, :p, 0, 1 )
            pcount-=1
            trigger_split_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
            #trigger_split_function=mk_function( :(  (life_time, p) -> nothing ) )
            trigger_split_function=@RuntimeGeneratedFunction( :(  (life_time, p) -> nothing ) )
            checked=false
            check_set=zeros(Float64,100)
            entry=floatmin()
            for x in 1:100
                push!(check_set,entry)
                entry=nextfloat(entry)
            end
            while !checked
                checked=true
                try
                    #trigger_split_function=mk_function( :(  (life_time, p) -> $trigger_split_expression ) )
                    trigger_split_function=@RuntimeGeneratedFunction( :(  (life_time, p) -> $trigger_split_expression ) )
                catch
                    checked=false
                end
                if checked
                    for lt in check_set
                        try
                            trigger_split_function(lt,trigger_split_expression_parameter)
                        catch
                            checked=false
                            break
                        end
                    end
                end
                if !checked
                    (trigger_split_expression,pcount)=create_expression( :life_time, :p, 0, 1 )
                    pcount-=1
                    trigger_split_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
                end
            end

            (housekeeping_expression,pcount)=create_expression( :energy, :p, 0, 1 )
            pcount-=1
            housekeeping_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
            #housekeeping_function=mk_function( :(  (energy, p) -> nothing ) )
            housekeeping_function=@RuntimeGeneratedFunction( :(  (energy, p) -> nothing ) )
            checked=false
            while !checked
                checked=true
                try
                    #housekeeping_function=mk_function( :(  (energy, p) -> $housekeeping_expression ) )
                    housekeeping_function=@RuntimeGeneratedFunction( :(  (energy, p) -> $housekeeping_expression ) )
                catch
                    checked=false
                end
                if checked
                    for e in org_max_energie:-1:1
                        try
                            floor(Int,housekeeping_function(e,housekeeping_expression_parameter))
                        catch
                            checked=false
                            break
                        end
                    end
                end
                if !checked
                    (housekeeping_expression,pcount)=create_expression( :energy, :p, 0, 1 )
                    pcount-=1
                    housekeeping_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
                    #housekeeping_function=mk_function( :(  (energy, p) -> $housekeeping_expression ) )
                end
            end

            (energy_expression,pcount)=create_expression( :external_energy, :p, 0, 1 )
            pcount-=1
            energy_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
            #energy_function=mk_function( :(  (external_energy, p) -> nothing ) )
            energy_function=@RuntimeGeneratedFunction( :(  (external_energy, p) -> nothing ) )
            checked=false
            while !checked
                checked=true
                try
                    #energy_function=mk_function( :(  (external_energy, p) -> $energy_expression ) )
                    energy_function=@RuntimeGeneratedFunction( :(  (external_energy, p) -> $energy_expression ) )
                catch e
                    checked=false
                end
                if checked
                    for external_energy in 0:10:100
                        try
                            floor(Int,energy_function(external_energy,energy_expression_parameter))
                        catch e
                            checked=false
                            break
                        end
                    end
                end
                if !checked
                    (energy_expression,pcount)=create_expression( :external_energy, :p, 0, 1 )
                    pcount-=1
                    energy_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
                    #energy_function=mk_function( :(  (external_energy, p) -> $energy_expression ) )
                end
            end

            complexity=1

            (ϕ_expression,pcount)=create_expression( :ϕ, :p, complexity, 1 )
            pcount-=1
            ϕ_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
            #ϕ_function=mk_function( :(  (ϕ, p) -> nothing ) )
            ϕ_function=@RuntimeGeneratedFunction( :(  (ϕ, p) -> nothing ) )
            #ϕ_function=Base.eval( :( (ϕ,p) -> $ϕ_expression ) )
            checked=false
            while !checked
                checked=true
                try
                    #ϕ_function=mk_function( :(  (ϕ, p) -> $ϕ_expression ) )
                    ϕ_function=@RuntimeGeneratedFunction( :(  (ϕ, p) -> $ϕ_expression ) )
                catch e
                    checked=false
                end
                #println("ϕ_expression : ",ϕ_expression," : ",ϕ_expression_parameter)
                if checked
                    for pos_ϕ in posCoord.ϕs
                        #Base.invokelatest(ϕ_function,pos_ϕ,ϕ_expression_parameter)
                        try
                            #Base.invokelatest(ϕ_function,pos_ϕ,ϕ_expression_parameter)
                            new_ϕ=ϕ_function(pos_ϕ,ϕ_expression_parameter)
                            if isnan(new_ϕ)
                                checked=false
                            end
                        catch e
                            checked=false
                            break
                        end
                    end
                end
                if ! checked
                    (ϕ_expression,pcount)=create_expression( :ϕ, :p, complexity, 1 )
                    pcount-=1
                    ϕ_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
                    #ϕ_function=Base.eval( :( (ϕ,p) -> $ϕ_expression ) )
                    #ϕ_function=mk_function( :(  (ϕ, p) -> $ϕ_expression ) )
                end
            end

            (θ_expression,pcount)=create_expression( :θ, :p, complexity, 1 )
            pcount-=1
            θ_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
            #θ_function=mk_function( :( (θ,p) -> nothing ) )
            θ_function=@RuntimeGeneratedFunction( :( (θ,p) -> nothing ) )
            #θ_function=Base.eval( :( (θ,p) -> $θ_expression ) )
            checked=false
            while !checked
                checked=true
                try
                    #θ_function=mk_function( :( (θ,p) -> $θ_expression ) )
                    θ_function=@RuntimeGeneratedFunction( :( (θ,p) -> $θ_expression ) )
                catch e
                    checked=false
                end
                #println("θ_expression : ",θ_expression," : ",θ_expression_parameter)
                if checked
                    for pos_θ in posCoord.θs
                        try
                            new_θ=θ_function(pos_θ,θ_expression_parameter)
                            #Base.invokelatest(θ_function,pos_θ,θ_expression_parameter)
                            if isnan(new_θ)
                                checked=false
                            end
                        catch e
                            checked=false
                            break
                        end
                    end
                end
                if ! checked
                    (θ_expression,pcount)=create_expression( :θ, :p, complexity, 1 )
                    pcount-=1
                    θ_expression_parameter=rand(-5.0:0.1:5.0)*rand(Float64,pcount)
                    #θ_function=Base.eval( :( (θ,p) -> $θ_expression ) )
                    #θ_function=mk_function( :( (θ,p) -> $θ_expression ) )
                end
            end

            current_position=get_random_coordinates(posCoord)
            new(id,
                energy,birth_time,
                is_sibling,generation,
                trigger_split_expression,trigger_split_expression_parameter,trigger_split_function,
                housekeeping_expression,housekeeping_expression_parameter,housekeeping_function,
                energy_expression,energy_expression_parameter,energy_function,
                ϕ_expression,ϕ_expression_parameter,ϕ_function,
                θ_expression,θ_expression_parameter,θ_function,
                current_position
            )
        end

        function Organism(org, randomize = false, delta = 0.0)
            id=next_id
            next_id+=1
            trigger_split_expression_parameter=deepcopy(org.trigger_split_expression_parameter)
            housekeeping_expression_parameter=deepcopy(org.housekeeping_expression_parameter)
            energy_expression_parameter=deepcopy(org.energy_expression_parameter)
            ϕ_expression_parameter=deepcopy(org.ϕ_expression_parameter)
            θ_expression_parameter=deepcopy(org.θ_expression_parameter)
            if randomize
                randomize!(trigger_split_expression_parameter,delta)
                randomize!(housekeeping_expression_parameter,delta)
                randomize!(energy_expression_parameter,delta)
                randomize!(ϕ_expression_parameter,delta)
                randomize!(θ_expression_parameter,delta)
            end
            new(id,
                org.energy,org.birth_time,
                org.is_sibling,org.generation,
                org.trigger_split_expression,trigger_split_expression_parameter,org.trigger_split_function,
                org.housekeeping_expression,housekeeping_expression_parameter,org.housekeeping_function,
                org.energy_expression,energy_expression_parameter,org.energy_function,
                org.ϕ_expression,ϕ_expression_parameter,org.ϕ_function,
                org.θ_expression,θ_expression_parameter,org.θ_function,
                org.current_position
            )
        end
    end
end


text:
using DiscreteEvents, Random, Distributions, GLMakie, LaTeXStrings, DataFrames

abstract type MState end

struct Idle <: MState end
struct Busy <: MState end

mutable struct Job
    no::Int64
    ts::Float64
    t1::Float64
    t2::Float64
    t3::Float64
end

mutable struct Machine
    state::MState
    job
end

Q = Job[]   # input queue
S = Job[]   # stock
M = Machine(Idle(), 0)
df = DataFrame(time = Float64[], buffer=Int[], machine=Int[], finished=Int[])
count = 1
printing = true

stats() = push!(df, (tau(), length(Q), M.state == Busy() ? 1 : 0, length(S)))

function arrive(μ, σ, c)
    @assert μ ≥ 1 "μ must be ≥ 1"
    ts = rand(Normal(μ, σ))/c
    job = Job(count, ts, tau(), 0, 0)
    global count += 1
    push!(Q, job)
    ta = rand(Erlang())*μ
    event!(𝐶, fun(arrive, μ, σ, c), after, ta)  # we schedule the next arrival
    printing ? println(tau(), ": job $(job.no) has arrived") : nothing # tau() is the current time
    if M.state == Idle()
        load()
    else
        stats()
    end
end

function load()
    M.state = Busy()
    M.job = popfirst!(Q)
    M.job.t2 = tau()
    event!(𝐶, fun(unload), after, M.job.ts)  # we schedule the unload
    printing ? println(tau(), ": job $(M.job.no) has been loaded") : nothing
    stats()
end

function unload()
    M.state = Idle()
    M.job.t3 = tau()
    push!(S, M.job)
    printing ? println(tau(), ": job $(M.job.no) has been finished") : nothing
    stats()
    M.job = 0
    if !isempty(Q)
        load()
    end
end

sample_time!(𝐶, 0.1)  # we determine the sample rate
periodic!(𝐶, fun(stats));  # we register stats() as sampling function

Random.seed!(2019)
#resetClock!(𝐶)
arrive(5, 1/5, 1)  # we schedule the first event
run!(𝐶, 30)        # and run the simulation

# FigureAxisPlot takes figure and axis keywords
fig, ax, p = lines((df.time, df.buffer),
    figure = (resolution = (1600, 600),),
    axis = (xlabel = "time (min)", ylabel = "buffer",title = "Single Server"),
    color = :blue,
    linewidth = 3)
fig

fig, ax, p = lines((df.time, df.machine),
    figure = (resolution = (1600, 600),),
    axis = (xlabel = "time (min)", ylabel = "machine",title = "Single Server"),
    color = :red,
    linewidth = 3)
fig

fig, ax, p = lines((df.time, df.finished),
    figure = (resolution = (1600, 600),),
    axis = (xlabel = "time (min)", ylabel = "stock",title = "Single Server"),
    color = :red,
    linewidth = 3)
fig

max_stars_repo_name: JMSoukup/github-demo

max_issues_repo_name: JMSoukup/github-demo

max_forks_repo_name: JMSoukup/github-demo


text:
function get_data(params::Dict{String, Any})
    n = params["num_nodes"]
    instance = params["instance"]

    file_path = joinpath(dirname(pathof(LaplacianOpt)),"..", "examples/instances/$(n)_nodes/$(n)_$(instance).json")

    data_dict = JSON.parsefile(file_path)

    if data_dict["num_nodes"] != n
        Memento.error(_LOGGER, "Mismatch in number of input nodes")
    end

    # Input data
    if "data_type" in keys(params)
        if params["data_type"] == "old"
            edge_weights = data_dict["old_data"]["edge_weights"]
        elseif params["data_type"] == "new"
            edge_weights = data_dict["new_data"]["edge_weights"]
        end
    else
        # default value
        edge_weights = data_dict["old_data"]["edge_weights"]
    end

    edge_weigths_matrix = LOpt.convert_array_to_matrix(n, edge_weights)

    # Solution type
    if "solution_type" in keys(params)
        solution_type = params["solution_type"]
    else
        solution_type = "exact"
    end

    # Relax Integrality
    if "relax_integrality" in keys(params)
        relax_integrality = params["relax_integrality"]
    else
        # default value
        relax_integrality = false
    end

    # Tolerance to verify zero values
    if "tol_zero" in keys(params)
        tol_zero = params["tol_zero"]

        if params["tol_zero"] >= 1E-4
            Memento.warn(_LOGGER, "Zero tolerance value may be too high for numerical accuracy in solutions")
        end

    else
        # default value
        tol_zero = 1E-6
    end

    # Tolerance to verify PSD-ness of a matrix
    if "tol_psd" in keys(params)
        tol_psd = params["tol_psd"]

        if params["tol_psd"] >= 1E-3
            Memento.warn(_LOGGER, "PSD tolerance value may be too high for PSD feasibility")
        end

    else
        # default value
        tol_psd = 1E-6
    end

    if "eigen_cuts_full" in keys(params)
        eigen_cuts_full = params["eigen_cuts_full"]
        if eigen_cuts_full
            Memento.info(_LOGGER, "Applying full-sized eigen cuts")
        end
    else
        #default value
        Memento.info(_LOGGER, "Applying full-sized eigen cuts")
        eigen_cuts_full = true
    end

    if "soc_linearized_cuts" in keys(params)
        soc_linearized_cuts = params["soc_linearized_cuts"]
        if soc_linearized_cuts
            Memento.info(_LOGGER, "Applying linearized SOC cuts (2x2 minors)")
        end
    else
        #default value
        soc_linearized_cuts = false
    end

    if "topology_flow_cuts" in keys(params)
        topology_flow_cuts = params["topology_flow_cuts"]
    else
        #default value
        Memento.info(_LOGGER, "Applying topology flow cuts")
        topology_flow_cuts = true
    end

    if eigen_cuts_full || topology_flow_cuts || soc_linearized_cuts
        lazy_callback_status = true
    end

    if "lazycuts_logging" in keys(params)
        lazycuts_logging = params["lazycuts_logging"]
    else
        #default value
        lazycuts_logging = false
    end

    data = Dict{String, Any}("num_nodes" => n,
                             "instance" => instance,
                             "edge_weights" => edge_weigths_matrix,
                             "solution_type" => solution_type,
                             "tol_zero" => tol_zero,
                             "tol_psd" => tol_psd,
                             "eigen_cuts_full" => eigen_cuts_full,
                             "soc_linearized_cuts" => soc_linearized_cuts,
                             "lazycuts_logging" => lazycuts_logging,
                             "topology_flow_cuts" => topology_flow_cuts,
                             "lazy_callback_status" => lazy_callback_status,
                             "relax_integrality" => relax_integrality)

    # Optimizer
    if "optimizer" in keys(params)
        data["optimizer"] = params["optimizer"]
    end

    return data
end

function convert_array_to_matrix(n::Int, edge_weights::Any)

    instance_matrix = Matrix{Float64}(zeros(n,n))

    for k = 1:length(edge_weights)
        i, j = edge_weights[k][1]
        w_ij = edge_weights[k][2]

        if w_ij < 0
            Memento.error(_LOGGER, "LaplacianOpt does not support graphs with negative weights")
        end

        instance_matrix[i,j] = w_ij
        instance_matrix[j,i] = w_ij
    end

    return instance_matrix
end


text:
using Simplices, SimplexSplitting, StaticArrays

function potentially_intersecting_simplices(t::Triangulation, image_i::Int)
    inds_potential_simplices = Int[]

    n_simplices = length(t.radii)

    @inbounds for i = 1:n_simplices
        dist_difference = ((t.centroids_im[image_i] - t.centroids[i]).' *
                            (t.centroids_im[image_i] - t.centroids[i]) - (t.radii_im[image_i] + t.radii[i])^2)[1]
        if dist_difference < 0
            push!(inds_potential_simplices, i)
        end
    end
    return inds_potential_simplices
end

function contains_point!(signs, s_arr, s, point)

    # Keep track of signs of determinants. If all have the same sign, the point
    # lies inside the simplex. If not, it lies outside.
    n_pts = size(s, 1)
    # Compute the first sign
    @views s_arr[1, 1:(end-1)] = point
    @views s_arr[2:end, 1:(end-1)] = s[2:end, :]

    signs[1] = sign(det(s_arr))

    for i = 2:(n_pts - 1) # Check remaining signs and stop if sign changes
        @views s_arr[1:(i - 1), 1:(end-1)] = s[1:(i - 1), :]
        @views s_arr[i, 1:(end-1)] = point
        @views s_arr[(i + 1):end, 1:(end-1)] = s[(i + 1):end, :]

        signs[i] = sign(det(s_arr))
        if !(signs[i] == signs[i - 1])
           return false
        end
    end

    @views s_arr[1:(end - 1), 1:(end-1)] = s[1:(end - 1), :]
    @views s_arr[end, 1:(end-1)] = point

    signs[end] = sign(det(s_arr))
    if !(signs[end] == signs[end - 1])
       return false
    end
    return true
end

function mm_dd(t::Triangulation;
    n_randpts::Int = 100,
    dist::Distributions.Distribution = Distributions.Uniform(0, 1),
    sample_randomly::Bool = false)

    n_simplices = size(t.simplex_inds, 1)
    n_vertices = size(t.simplex_inds, 2)
    dim = n_vertices - 1
    simplices = Size(n_simplices, dim + 1, dim)(zeros(n_simplices, dim + 1, dim))
    imsimplices = Size(n_simplices, dim + 1, dim)(zeros(n_simplices, dim + 1, dim))

    for i in 1:n_simplices
        simplices[i, :, :] = t.points[t.simplex_inds[i, :], :]
        imsimplices[i, :, :] = t.impoints[t.simplex_inds[i, :], :]
    end

    convex_coeffs = get_coeffs(dim, n_randpts, sample_randomly).'
    n_coeffs = maximum(size(convex_coeffs))

    pt = Size(dim)(zeros(Float64, dim))
    cache_arr = Size(dim + 1, dim + 1)(zeros(Float64, 4, 4))
    cache_arr[:, dim+1] = ones(dim+1)

    s_arr = Size(dim, dim + 1)(zeros(Float64, dim, dim + 1))
    signs = Size(dim + 1)(zeros(Float64, dim + 1))

    M = zeros(Float64, n_simplices, n_simplices)

    for i in 1:n_simplices
        inds = potentially_intersecting_simplices(t, i)
        image = imsimplices[i, :, :]

        for c in 1:n_coeffs
            @views pt = convex_coeffs[c, :].' * image
            for j in inds
                simplex = @views simplices[j, :, :]
                if contains_point!(signs, cache_arr, simplex, pt)
                    M[j, i] += 1.0
                end
            end
        end
    end
    return M.' ./ n_coeffs
end


text:
using Test

include("./data/turtle_ex1.jl")

@testset "CSerd" begin
  include("CSerd.jl")
end

@testset "Serd" begin
  include("Serd.jl")
end


text:
"""
    $(SIGNATURES)

Calculate ΔG' of a reaction described by `rxn_string` using the settings in
`equilibrator`. Optionally, set the concentrations (abundances) of the species involved
in the reaction through `concentration`. If `skip_unbalanced` is true then
return `nothing` if the reaction is unbalanced. If `balance_warn` is false
then do not emit a warning when the reaction is unbalanced. Note, concentrations
that are smaller than 0.001 mM are capped at 0.001 mM (an eQuilibrator/equilibrator_api)
restriction.

# Example
```
equilibrator = Equilibrator()

rxn_string = "bigg.metabolite:atp + bigg.metabolite:h2o = bigg.metabolite:adp + bigg.metabolite:pi"
concens = Dict("bigg.metabolite:atp"=>1u"mM", "bigg.metabolite:adp"=>100u"μM", "bigg.metabolite:pi"=>0.003u"M")
dg_prime(equilibrator, rxn_string; concentrations=concens)
```
"""
function dg_prime(
    equilibrator::Equilibrator,
    rxn_string::String;
    concentrations = Dict{String,Real}(),
    balance_warn = true,
    skip_unbalanced = false,
)
    rxn = equilibrator.cc.parse_reaction_formula(rxn_string)

    if balance_warn && !rxn.is_balanced()
        @warn "Reaction unbalanced."
        skip_unbalanced && return nothing
    end

    for (k, v) in concentrations
        val = uconvert(u"M", _lower_bound(float(v)))
        rxn.set_abundance(equilibrator.cc.get_compound(k), equilibrator.eq.Q_(string(val)))
    end
    pydg = equilibrator.cc.dg_prime(rxn)
    dg = pydg.value.magnitude ± pydg.error.magnitude
    units = uparse(_parse_units(py"str"(pydg.units)))
    return (dg)units
end

"""
    $(SIGNATURES)

Calculate ΔG' at physiological conditions (all concentrations set to 1 mM, 1mbar) of a reaction
described by `rxn_string` using the settings in `equilibrator`. If `skip_unbalanced` is true
then return `nothing` if the reaction is unbalanced. If `balance_warn` is false then do not
emit a warning when the reaction is unbalanced.

# Example
```
equilibrator = Equilibrator()

rxn_string = "bigg.metabolite:atp + bigg.metabolite:h2o = bigg.metabolite:adp + bigg.metabolite:pi"
physiological_dg_prime(equilibrator, rxn_string)
```
"""
function physiological_dg_prime(
    equilibrator::Equilibrator,
    rxn_string::String;
    balance_warn = true,
    skip_unbalanced = false,
)
    rxn = equilibrator.cc.parse_reaction_formula(rxn_string)

    if balance_warn && !rxn.is_balanced()
        @warn "Reaction unbalanced."
        skip_unbalanced && return nothing
    end

    pydg = equilibrator.cc.physiological_dg_prime(rxn)
    dg = pydg.value.magnitude ± pydg.error.magnitude
    units = uparse(_parse_units(py"str"(pydg.units)))
    return (dg)units
end

"""
    $(SIGNATURES)

Calculate ΔG⁰ at standard conditons of a reaction described by `rxn_string` using
the settings in `equilibrator`. If `skip_unbalanced` is true then return `nothing` if the
reaction is unbalanced. If `balance_warn` is false then do not emit a warning
when the reaction is unbalanced.

Note, standard conditions means that all concentrations are set at 1 M.

# Example
```
equilibrator = Equilibrator()

rxn_string = "bigg.metabolite:atp + bigg.metabolite:h2o = bigg.metabolite:adp + bigg.metabolite:pi"
standard_dg_prime(equilibrator, rxn_string)
```
"""
function standard_dg_prime(
    equilibrator::Equilibrator,
    rxn_string::String;
    balance_warn = true,
    skip_unbalanced = false,
)
    rxn = equilibrator.cc.parse_reaction_formula(rxn_string)

    if balance_warn && !rxn.is_balanced()
        @warn "Reaction unbalanced."
        skip_unbalanced && return nothing
    end

    pydg = equilibrator.cc.standard_dg_prime(rxn)
    dg = pydg.value.magnitude ± pydg.error.magnitude
    units = uparse(_parse_units(py"str"(pydg.units)))
    return (dg)units
end

"""
    $(SIGNATURES)

Calculate the log of the reversibility index of a reaction described by `rxn_string`
using settings in `equilibrator`. If `skip_unbalanced` is true then return `nothing` if the reaction is
unbalanced. If `balance_warn` is false then do not emit a warning when the
reaction is unbalanced.

# Example
```
equilibrator = Equilibrator()

rxn_string = "bigg.metabolite:atp + bigg.metabolite:h2o = bigg.metabolite:adp + bigg.metabolite:pi"
ln_reversibility_index(equilibrator, rxn_string)
```
"""
function ln_reversibility_index(
    equilibrator::Equilibrator,
    rxn_string::String;
    balance_warn = true,
    skip_unbalanced = false,
)
    rxn = equilibrator.cc.parse_reaction_formula(rxn_string)

    if balance_warn && !rxn.is_balanced()
        @warn "Reaction unbalanced."
        skip_unbalanced && return nothing
    end

    rev_obj = equilibrator.cc.ln_reversibility_index(rxn)
    return rev_obj.value.magnitude ± rev_obj.error.magnitude
end


text:
## Some simple tests of the package

using Mustache
mtrender = render

tpl = mt"the value of x is {{x}} and that of y is {{y}}"

## a dict
out = mtrender(tpl, Dict("x"=>1, "y"=>2))
println(out)

## A module
x = 1; y = "two"
mtrender(tpl, Main)

## a CompositeKind
mutable struct ThrowAway
    x
    y
end

mtrender(tpl, ThrowAway("ex","why"))

## a more useful CompositeKind
using Distributions
tpl = mt"Beta distribution with alpha={{alpha}}, beta={{beta}}"
mtrender(tpl, Beta(1, 2))

## conditional text
using Mustache
tpl = "{{#b}}this doesn't show{{/b}}{{#a}}this does show{{/a}}"
mtrender(tpl, Dict("a" => 1))

## We can iterate over data frames. Handy for making tables
using Mustache
using DataFrames

## SHow values in Main in a web page

_names = Array(String, 0)
_summaries = Array(String, 0)
m = Main
for s in sort(map(string, names(m)))
    v = Symbol(s)
    if isdefined(m,v)
        push!(_names, s)
        push!(_summaries, summary(eval(m,v)))
    end
end

using DataFrames
d = DataFrame(names = _names, summs = _summaries)

tpl = "

<title>{{Title}}</title>

name	summary

{{#d}}

{{names}}	{{summs}}

{{/d}}

";

out = mtrender(tpl, {"Title" => "A quick table", "d" => d})
## show in browser (on Mac)
f = tempname()
io = open("$f.html", "w")
print(io, out)
close(io)
run(`open $f.html`)

## A dict using symbols
d = { :a => 1, :b => 2}
tpl = "symbol {{:a}} and {{:b}}"
mtrender(tpl, d)

## array of Dicts
using Mustache

A = [Dict("a" => "eh", "b" => "bee"),
     Dict("a" => "ah", "b" => "buh")]

## Contrast to data frame:
# D = DataFrame(quote
#   a = ["eh", "ah"]
#   b = ["bee", "buh"]
# end);

tpl = mt"{{#A}} pronounce a as {{a}} and b as {{b}}.{{/A}}"

mtrender(tpl, Dict("A" => A))

## use .[ind] to index within a vector:

tpl = mt"{{#:vec}}{{.[1]}}{{/:vec}}" # just first one
mtrender(tpl, vec=["A1","B2","C3"])

tpl = mt"{{#:vec}}{{.}}{{^.[end]}}, {{/.[end]}}{{/:vec}}" # serial commas
mtrender(tpl, vec=["A1","B2","C3"])


text:
######## Internal Simulation Object Structs ########
mutable struct StageInternal
    number::Int
    executions::Int
    execution_count::Int
    end_of_interval_step::Int
    # This line keeps track of the executions of a stage relative to other stages.
    # This might be needed in the future to run multiple stages. For now it is disabled
    #synchronized_executions::Dict{Int, Int} # Number of executions per upper level stage step
    psi_container::PSIContainer
    # Caches are stored in set because order isn't relevant and they should be unique
    caches::Set{CacheKey}
    chronolgy_dict::Dict{Int, <:FeedForwardChronology}
    built::BUILD_STATUS
    stage_path::String
    ext::Dict{String, Any}
    function StageInternal(
        number,
        executions,
        execution_count,
        psi_container;
        ext = Dict{String, Any}(),
    )
        new(
            number,
            executions,
            execution_count,
            0,
            psi_container,
            Set{CacheKey}(),
            Dict{Int, FeedForwardChronology}(),
            false,
            "",
            ext,
        )
    end
end

# TODO: Add DocString
@doc raw"""
    Stage({M<:AbstractOperationsProblem}
        template::OperationsProblemTemplate
        sys::PSY.System
        optimizer::JuMP.MOI.OptimizerWithAttributes
        internal::Union{Nothing, StageInternal}
        )
"""
mutable struct Stage{M <: AbstractOperationsProblem}
    template::OperationsProblemTemplate
    sys::PSY.System
    internal::Union{Nothing, StageInternal}

    function Stage{M}(
        template::OperationsProblemTemplate,
        sys::PSY.System,
        settings::PSISettings,
        jump_model::Union{Nothing, JuMP.AbstractModel} = nothing,
    ) where {M <: AbstractOperationsProblem}
        internal = StageInternal(0, 0, 0, PSIContainer(sys, settings, jump_model))
        new{M}(template, sys, internal)
    end
end

function Stage{M}(
    template::OperationsProblemTemplate,
    sys::PSY.System,
    optimizer::JuMP.MOI.OptimizerWithAttributes,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
) where {M <: AbstractOperationsProblem}
    check_kwargs(kwargs, STAGE_ACCEPTED_KWARGS, "Stage")
    settings = PSISettings(sys; optimizer = optimizer, use_parameters = true, kwargs...)
    return Stage{M}(template, sys, settings, jump_model)
end

"""
    Stage(::Type{M},
    template::OperationsProblemTemplate,
    sys::PSY.System,
    optimizer::JuMP.MOI.OptimizerWithAttributes,
    jump_model::Union{Nothing, JuMP.AbstractModel}=nothing;
    kwargs...) where {M<:AbstractOperationsProblem}
This builds the optimization problem of type M with the specific system and template for the simulation stage
# Arguments
- `::Type{M} where M<:AbstractOperationsProblem`: The abstract operation model type
- `template::OperationsProblemTemplate`: The model reference made up of transmission, devices,
                                          branches, and services.
- `sys::PSY.System`: the system created using Power Systems
- `jump_model::Union{Nothing, JuMP.AbstractModel}`: Enables passing a custom JuMP model. Use with care
# Output
- `Stage::Stage`: The operation model containing the model type, unbuilt JuMP model, Power
Systems system.
# Example
```julia
template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
stage = Stage(MyOpProblemType template, system, optimizer)
```
# Accepted Key Words
- `initial_time::Dates.DateTime`: Initial Time for the model solve
- `PTDF::PTDF`: Passes the PTDF matrix into the optimization model for StandardPTDFModel networks.
- `warm_start::Bool` True will use the current operation point in the system to initialize variable values. False initializes all variables to zero. Default is true
- `slack_variables::Bool` True will add slacks to the system balance constraints
"""
function Stage(
    ::Type{M},
    template::OperationsProblemTemplate,
    sys::PSY.System,
    optimizer::JuMP.MOI.OptimizerWithAttributes,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
) where {M <: AbstractOperationsProblem}
    return Stage{M}(template, sys, optimizer, jump_model; kwargs...)
end

function Stage(
    template::OperationsProblemTemplate,
    sys::PSY.System,
    optimizer::JuMP.MOI.OptimizerWithAttributes,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
)
    return Stage{GenericOpProblem}(template, sys, optimizer, jump_model; kwargs...)
end

stage_built(s::Stage) = s.internal.built == BUILT
stage_empty(s::Stage) = s.internal.built == EMPTY
get_execution_count(s::Stage) = s.internal.execution_count
get_executions(s::Stage) = s.internal.executions
get_sys(s::Stage) = s.sys
get_template(s::Stage) = s.template
get_number(s::Stage) = s.internal.number
get_psi_container(s::Stage) = s.internal.psi_container
get_end_of_interval_step(s::Stage) = s.internal.end_of_interval_step
warm_start_enabled(s::Stage) = get_warm_start(s.internal.psi_container.settings)
get_initial_time(s::Stage{T}) where {T <: AbstractOperationsProblem} =
    get_initial_time(s.internal.psi_container.settings)

function reset!(stage::Stage{M}) where {M <: AbstractOperationsProblem}
    @assert stage_built(stage)
    if stage_built(stage)
        @info("Stage $(stage.internal.number) will be reset by the simulation build call")
    end
    stage.internal.execution_count = 0
    stage.internal.psi_container =
        PSIContainer(stage.sys, stage.internal.psi_container.settings, nothing)
    stage.internal.built = EMPTY
    return
end

function build!(
    stage::Stage{M},
    initial_time::Dates.DateTime,
    horizon::Int,
    stage_interval::Dates.Period,
) where {M <: PowerSimulationsOperationsProblem}
    !stage_empty(stage) && reset!(stage)
    settings = get_settings(get_psi_container(stage))
    # Horizon and initial time are set here because the information is specified in the
    # Simulation Sequence object and not at the stage creation.
    set_horizon!(settings, horizon)
    set_initial_time!(settings, initial_time)
    stage.internal.built = IN_PROGRESS
    psi_container = get_psi_container(stage)
    _build!(psi_container, stage.template, stage.sys)
    @assert get_horizon(psi_container.settings) == length(psi_container.time_steps)
    stage_resolution = PSY.get_forecasts_resolution(stage.sys)
    stage.internal.end_of_interval_step = Int(stage_interval / stage_resolution)
    stage_path = stage.internal.stage_path
    _write_psi_container(
        stage.internal.psi_container,
        joinpath(stage_path, "Stage$(stage.internal.number)_optimization_model.json"),
    )
    if get_system_to_file(settings)
        PSY.to_json(
            stage.sys,
            joinpath(stage_path, "Stage$(stage.internal.number)_sys_data.json"),
        )
    end
    stage.internal.built = BUILT
    return
end

function run_stage(
    stage::Stage{M},
    start_time::Dates.DateTime,
    results_path::String,
) where {M <: PowerSimulationsOperationsProblem}
    @assert stage.internal.psi_container.JuMPmodel.moi_backend.state != MOIU.NO_OPTIMIZER
    timed_log = Dict{Symbol, Any}()
    model = stage.internal.psi_container.JuMPmodel
    _, timed_log[:timed_solve_time], timed_log[:solve_bytes_alloc], timed_log[:sec_in_gc] =
        @timed JuMP.optimize!(model)

    @info "JuMP.optimize! completed" timed_log

    model_status = JuMP.primal_status(stage.internal.psi_container.JuMPmodel)
    if model_status != MOI.FEASIBLE_POINT::MOI.ResultStatusCode
        error("Stage $(stage.internal.number) status is $(model_status)")
    end
    # TODO: Add Fallback when optimization fails
    # if is_milp(stage.internal.psi_container)
    _export_model_result(stage, start_time, results_path)
    _export_optimizer_log(timed_log, stage.internal.psi_container, results_path)
    stage.internal.execution_count += 1
    # Reset execution count
    if stage.internal.execution_count == stage.internal.executions
        stage.internal.execution_count = 0
    end
    return
end

# Here because requires the stage to be defined
# This is a method a user defining a custom cache will have to define. This is the definition
# in PSI for the building the TimeStatusChange
function get_initial_cache(cache::AbstractCache, stage::Stage)
    throw(ArgumentError("Initialization method for cache $(typeof(cache)) not defined"))
end

function get_initial_cache(cache::TimeStatusChange, stage::Stage)
    ini_cond_on = get_initial_conditions(
        stage.internal.psi_container,
        TimeDurationON,
        cache.device_type,
    )

    ini_cond_off = get_initial_conditions(
        stage.internal.psi_container,
        TimeDurationOFF,
        cache.device_type,
    )

    device_axes = Set((
        PSY.get_name(ic.device) for ic in Iterators.Flatten([ini_cond_on, ini_cond_off])
    ),)
    value_array = JuMP.Containers.DenseAxisArray{Dict{Symbol, Float64}}(undef, device_axes)

    for ic in ini_cond_on
        device_name = PSY.get_name(ic.device)
        condition = get_condition(ic)
        status = (condition > 0.0) ? 1.0 : 0.0
        value_array[device_name] = Dict(:count => condition, :status => status)
    end

    for ic in ini_cond_off
        device_name = PSY.get_name(ic.device)
        condition = get_condition(ic)
        status = (condition > 0.0) ? 0.0 : 1.0
        if value_array[device_name][:status] != status
            throw(IS.ConflictingInputsError("Initial Conditions for $(device_name) are not compatible. The values provided are invalid"))
        end
    end

    return value_array
end

function get_initial_cache(cache::StoredEnergy, stage::Stage)
    ini_cond_level =
        get_initial_conditions(stage.internal.psi_container, EnergyLevel, cache.device_type)

    device_axes = Set((PSY.get_name(ic.device) for ic in ini_cond_level),)
    value_array = JuMP.Containers.DenseAxisArray{Float64}(undef, device_axes)
    for ic in ini_cond_level
        device_name = PSY.get_name(ic.device)
        condition = get_condition(ic)
        value_array[device_name] = condition
    end
    return value_array
end

function get_time_stamps(stage::Stage, start_time::Dates.DateTime)
    resolution = PSY.get_forecasts_resolution(stage.sys)
    horizon = stage.internal.psi_container.time_steps[end]
    range_time = collect(start_time:resolution:(start_time + resolution * horizon))
    time_stamp = DataFrames.DataFrame(Range = range_time[:, 1])

    return time_stamp
end

function write_data(stage::Stage, save_path::AbstractString; kwargs...)
    write_data(stage.internal.psi_container, save_path; kwargs...)
    return
end

# These functions are writing directly to the feather file and skipping printing to memory.
function _export_model_result(stage::Stage, start_time::Dates.DateTime, save_path::String)
    duals = Dict()
    if is_milp(stage.internal.psi_container)
        @warn("Stage $(stage.internal.number) is an MILP, duals can't be exported")
    else
        for c in get_constraint_duals(get_psi_container(stage).settings)
            v = get_constraint(get_psi_container(stage), c)
            duals[c] = axis_array_to_dataframe(v)
        end
    end
    write_data(stage, save_path)
    write_data(duals, save_path; duals = true)
    write_data(get_parameters_value(stage.internal.psi_container), save_path; params = true)
    write_data(get_time_stamps(stage, start_time), save_path, "time_stamp")
    files = collect(readdir(save_path))
    compute_file_hash(save_path, files)
    return
end

struct StageSerializationWrapper
    template::OperationsProblemTemplate
    sys::String
    settings::PSISettings
    stage_type::DataType
end


text:
using Dates
using FileIO
using MPI
using Printf
using Random
using StaticArrays
using Test

using CLIMA
using CLIMA.Atmos
using CLIMA.Diagnostics
using CLIMA.GenericCallbacks
using CLIMA.ODESolvers
using CLIMA.Mesh.Filters
using CLIMA.MoistThermodynamics
using CLIMA.ODESolvers
using CLIMA.VariableTemplates
using CLIMA.Writers

using CLIMAParameters
using CLIMAParameters.Planet: grav, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function init_sin_test!(bl, state, aux, (x, y, z), t)
    FT = eltype(state)

    z = FT(z)
    _grav::FT = grav(bl.param_set)
    _MSLP::FT = MSLP(bl.param_set)

    # These constants are those used by Stevens et al. (2005)
    qref = FT(9.0e-3)
    q_pt_sfc = PhasePartition(qref)
    Rm_sfc = FT(gas_constant_air(param_set, q_pt_sfc))
    T_sfc = FT(292.5)
    P_sfc = _MSLP

    # Specify moisture profiles
    q_liq = FT(0)
    q_ice = FT(0)
    zb = FT(600)         # initial cloud bottom
    zi = FT(840)         # initial cloud top
    dz_cloud = zi - zb
    q_liq_peak = FT(0.00045)     # cloud mixing ratio at z_i

    if z > zb && z <= zi
        q_liq = (z - zb) * q_liq_peak / dz_cloud
    end

    if z <= zi
        θ_liq = FT(289.0)
        q_tot = qref
    else
        θ_liq = FT(297.5) + (z - zi)^(FT(1 / 3))
        q_tot = FT(1.5e-3)
    end

    w = FT(10 + 0.5 * sin(2 * π * ((x / 1500) + (y / 1500))))
    u = (5 + 2 * sin(2 * π * ((x / 1500) + (y / 1500))))
    v = FT(5 + 2 * sin(2 * π * ((x / 1500) + (y / 1500))))

    # Pressure
    H = Rm_sfc * T_sfc / _grav
    p = P_sfc * exp(-z / H)

    # Density, Temperature
    ts = LiquidIcePotTempSHumEquil_given_pressure(bl.param_set, θ_liq, p, q_tot)
    #ρ = air_density(ts)
    ρ = one(FT)

    e_kin = FT(1 / 2) * FT((u^2 + v^2 + w^2))
    e_pot = _grav * z
    E = ρ * total_energy(e_kin, e_pot, ts)

    state.ρ = ρ
    state.ρu = SVector(ρ * u, ρ * v, ρ * w)
    state.ρe = E
    state.moisture.ρq_tot = ρ * q_tot
end

function config_sin_test(FT, N, resolution, xmax, ymax, zmax)
    ode_solver =
        CLIMA.ExplicitSolverType(solver_method = LSRK54CarpenterKennedy)
    config = CLIMA.AtmosLESConfiguration(
        "Diagnostics SIN test",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_sin_test!;
        solver_type = ode_solver,
    )

    return config
end

function config_diagnostics(driver_config)
    interval = "100steps"
    dgngrp = setup_atmos_default_diagnostics(
        interval,
        replace(driver_config.name, " " => "_"),
        writer = JLD2Writer(),
    )
    return CLIMA.DiagnosticsConfiguration([dgngrp])
end

function main()
    CLIMA.init()

    # Disable driver diagnostics as we're testing it here
    CLIMA.Settings.diagnostics = "never"

    FT = Float64

    # DG polynomial order
    N = 4

    # Domain resolution and size
    Δh = FT(50)
    Δv = FT(20)
    resolution = (Δh, Δh, Δv)

    xmax = FT(1500)
    ymax = FT(1500)
    zmax = FT(1500)

    t0 = FT(0)
    dt = FT(0.01)
    timeend = dt

    driver_config = config_sin_test(FT, N, resolution, xmax, ymax, zmax)
    solver_config = CLIMA.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
        init_on_cpu = true,
    )
    dgn_config = config_diagnostics(driver_config)

    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    Q = solver_config.Q
    solver = solver_config.solver

    outdir = mktempdir()
    currtime = ODESolvers.gettime(solver)
    starttime = replace(string(now()), ":" => ".")
    Diagnostics.init(mpicomm, dg, Q, starttime, outdir)
    dgn_config.groups[1](currtime, init = true)
    dgn_config.groups[1](currtime)

    CLIMA.invoke!(solver_config)

    # Check results
    mpirank = MPI.Comm_rank(mpicomm)
    if mpirank == 0
        dgngrp = dgn_config.groups[1]
        nm = @sprintf(
            "%s_%s_%s_num%04d.jld2",
            replace(dgngrp.out_prefix, " " => "_"),
            dgngrp.name,
            starttime,
            1,
        )
        ds = load(joinpath(outdir, nm))
        N = size(ds["u"], 1)
        err = 0
        err1 = 0
        for i in 1:N
            u = ds["u"][i]
            cov_w_u = ds["cov_w_u"][i]
            err += (cov_w_u - 0.5)^2
            err1 += (u - 5)^2
        end
        err = sqrt(err / N)
        @test err1 <= 1e-16
        @test err <= 2e-15
    end
end

main()


text:
concept_ordered(x) = issorted(x)

"""
    ordered(x)
Global constraint ensuring that all the values of `x` are ordered.
"""
const ordered = Constraint(
    concept = concept_ordered,
    error = make_error(:ordered),
)


text:
using Distributions: Normal

"""
    MomentumHarmonicWigner(ω, β)

Wigner distribution in a 1D harmonic potential for the momentum
"""
function MomentumHarmonicWigner(ω, β)
    σ = sqrt(Q(ω, β) / β)
    Normal(0, σ)
end

"""
    PositionHarmonicWigner(ω, β)

Wigner distribution in a 1D harmonic potential for the position
"""
function PositionHarmonicWigner(ω, β, m; centre=0)
    σ = sqrt(Q(ω, β) / β / m) / ω
    Normal(centre, σ)
end

"""
    VelocityHarmonicWigner(ω, β, m)

Wigner distribution in a 1D harmonic potential for the velocity
"""
function VelocityHarmonicWigner(ω, β, m)
    σ = sqrt(Q(ω, β) / β) / m
    Normal(0, σ)
end

"Quantum corrector for the Wigner distribution"
function Q(ω, β)
    halfβħω = β*ω/2
    return halfβħω / tanh(halfβħω)
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Shapes (Modia3D/Shapes/_module.jl)
#
# Content:
#   This file contains mutable visual geometries that are only used for animation
#   (immutable solid geometries that have a volume are defined in module Modia3D.Shapes.
#    They can be additionally used for collision handling and other purposes).

"""
    Sphere(; diameter=1.0)

Generate a new solid or visual shape representing a sphere.
The reference frame = Object3D frame is located at the center of the sphere.

# Arguments
- `diameter` defines the diameter of the sphere.
"""
mutable struct Sphere <: Modia3D.AbstractGeometry
    diameter::Float64

    function Sphere(; diameter=1.0)
        @assert(diameter >= 0.0)
        new(diameter)
    end
end

"""
    Ellipsoid(; lengthX=1.0, lengthY=1.0, lengthZ=1.0)

Generate a new solid or visual shape representing an ellipsoid.
The reference frame = Object3D frame is located at the center of the ellipsoid.

# Arguments
- `lengthX` defines the length of the ellipsoid in x-direction.
- `lengthY` defines the length of the ellipsoid in y-direction.
- `lengthZ` defines the length of the ellipsoid in z-direction.
"""
mutable struct Ellipsoid <: Modia3D.AbstractGeometry
    lengthX::Float64
    lengthY::Float64
    lengthZ::Float64

    function Ellipsoid(; lengthX=1.0, lengthY=1.0, lengthZ=1.0)
        @assert(lengthX >= 0.0)
        @assert(lengthY >= 0.0)
        @assert(lengthZ >= 0.0)
        if lengthX == lengthY == lengthZ
            Sphere(diameter=lengthX)
        else
            new(lengthX, lengthY, lengthZ)
        end
    end
end

"""
    Box(; lengthX=1.0, lengthY=1.0, lengthZ=1.0)

Generate a new solid or visual shape representing a box.
The reference frame = Object3D frame is located at the center of the box.

# Arguments
- `lengthX` defines the length of the box in x-direction.
- `lengthY` defines the length of the box in y-direction.
- `lengthZ` defines the length of the box in z-direction.
"""
mutable struct Box <: Modia3D.AbstractGeometry
    lengthX::Float64
    lengthY::Float64
    lengthZ::Float64

    function Box(; lengthX=1.0, lengthY=1.0, lengthZ=1.0)
        @assert(lengthX >= 0.0)
        @assert(lengthY >= 0.0)
        @assert(lengthZ >= 0.0)
        new(lengthX, lengthY, lengthZ)
    end
end

"""
    Cylinder(; axis=3, diameter=1.0, length=1.0, innerDiameter=0.0)

Generate a new solid or visual shape representing a cylinder.
The reference frame = Object3D frame is located at the center of the cylinder.

# Arguments
- `axis` defines the rotation axis of the cylinder: 1=x; 2=y; 3=z.
- `diameter` defines the diameter of the cylinder.
- `length` defines the length of the cylinder.
- `innerDiameter` defines the inner diameter of the cylinder (where `innerDiameter=0` defines a full cylinder).

# Notes
- `innerDiameter` is not supported by collision and animation export.
"""
mutable struct Cylinder <: Modia3D.AbstractGeometry
    axis::Int
    diameter::Float64
    length::Float64
    innerDiameter::Float64

    function Cylinder(; axis=3, diameter=1.0, length=1.0, innerDiameter=0.0)
        @assert(1 <= axis <= 3)
        @assert(diameter >= 0.0)
        @assert(length >= 0.0)
        @assert(0.0 <= innerDiameter < diameter)
        new(axis, diameter, length, innerDiameter)
    end
end

"""
    Cone(; axis=3, diameter=1.0, length=1.0, topDiameter=0.0)

Generate a new solid or visual shape representing a circular cone or a frustum of a circular cone.
The reference frame = Object3D frame is located at the center of the base circle of the cone/frustum.

# Arguments
- `axis` defines the rotation axis of the cone/frustum: 1=x; 2=y; 3=z.
- `diameter` defines the diameter of the base circle of the cone/frustum.
- `length` defines the length of the cone/frustum.
- `topDiameter` defines the diameter of the top circle of the frustum (where `topDiameter=0` defines a right cone).
"""
mutable struct Cone <: Modia3D.AbstractGeometry
    axis::Int
    diameter::Float64
    length::Float64
    topDiameter::Float64

    function Cone(; axis=3, diameter=1.0, length=1.0, topDiameter=0.0)
        @assert(1 <= axis <= 3)
        @assert(diameter >= 0.0)
        @assert(length >= 0.0)
        @assert(0 <= topDiameter < diameter)
        new(axis, diameter, length, topDiameter)
    end
end

"""
    Capsule(; axis=3, diameter=1.0, length=1.0)

Generate a new solid or visual shape representing a capsule assembled by a cylinder and two half spheres.
The reference frame = Object3D frame is located at the center of the capsule.

# Arguments
- `axis` defines the rotation axis of the capsule: 1=x; 2=y; 3=z.
- `diameter` defines the diameter of the capsule.
- `length` defines the length of the middle cylindrical part of the capsule.

# Notes
- SimVis visualizes [half ellipsoids with semi-axis length `length/2` instead of spheres with diameter `diameter`](https://github.com/ModiaSim/PrivateModia3D.jl/issues/54).
- Capsule is [not supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).
"""
mutable struct Capsule <: Modia3D.AbstractGeometry
    axis::Int
    diameter::Float64
    length::Float64

    function Capsule(; axis=3, diameter=1.0, length=1.0)
        @assert(1 <= axis <= 3)
        @assert(diameter >= 0.0)
        @assert(length >= 0.0)
        new(axis, diameter, length)
    end
end

"""
    Beam(; axis=3, length=1.0, width=0.2, thickness=0.1)

Generate a new solid or visual shape representing a beam assembled by a box and two half cylinders.
The reference frame = Object3D frame is located at the center of the beam.

# Arguments
- `axis` defines the longitudinal axis of the beam: 1=x; 2=y; 3=z.
- `length` defines the box length of the beam. The total length of the beam is `length + width`.
- `width` defines the width (= cylinder diameter) of the beam.
- `thickness` defines the thickness of the beam.

The dimension directions depend on `axis` by circularly shift:

| `axis` | `length` | `width` | `thickness` |
|:------:|:--------:|:-------:|:-----------:|
|   1    |   x      |   y     |   z         |
|   2    |   y      |   z     |   x         |
|   3    |   z      |   x     |   y         |

# Notes
- Beam is [not supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).
"""
mutable struct Beam <: Modia3D.AbstractGeometry
    axis::Int
    length::Float64
    width::Float64
    thickness::Float64

    function Beam(; axis=3, length=1.0, width=0.2, thickness=0.1)
        @assert(1 <= axis <= 3)
        @assert(length >= 0.0)
        @assert(width >= 0.0)
        @assert(thickness >= 0.0)
        new(axis, length, width, thickness)
    end
end

"""
    FileMesh(; filename::AbstractString="", scale=Basics.onesMVector(),
    useMaterialColor::Bool=false, smoothNormals::Bool=false, convexPartition::Bool=false)

Generate a new solid or visual shape representing a mesh.
The reference frame = Object3D frame is defined by the mesh data.

# Arguments
- `filename` defines the name of the mesh file.
- `scale` defines the scaling factors to be applied in x-, y- and z-direction.
- `useMaterialColor` defines if the material color of the shape is to be considered.
- `smoothNormals` defines if smoothing of mesh normals is active.
- `convexPartition` defines if partitioning into convex sub meshes is active.

# Notes
- [Solid](@ref) features are supported only for [.obj files](https://en.wikipedia.org/wiki/Wavefront_.obj_file).
- [Visual](@ref) features are supported for (3ds, dxf, obj, stl) - formats.
- FileMesh is [not directly supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).
  You have to convert your .obj - FileMesh into a .json format first and store it in the same folder as your obj file.
"""
mutable struct FileMesh <: Modia3D.AbstractGeometry
    filename::AbstractString
    scaleFactor::MVector{3,Float64}
    useMaterialColor::Bool
    smoothNormals::Bool
    convexPartition::Bool

    # for solids only
    centroid::MVector{3,Float64}
    longestEdge::Float64
    objPoints::Vector{MVector{3,Float64}}
    facesIndizes::Vector{AbstractVector{}}
    volume::Float64
    centroidAlgo::SVector{3,Float64}
    inertia::SMatrix{3,3,Float64,9}
    function FileMesh(; filename::AbstractString="", scale=Basics.onesMVector(),
        useMaterialColor::Bool=false, smoothNormals::Bool=false, convexPartition::Bool=false)
        if !isfile(filename)
        error("FileMesh(\"$filename\",...): file not found.")
        end
        @assert(scale[1] >= 0.0)
        @assert(scale[2] >= 0.0)
        @assert(scale[3] >= 0.0)
        fileMesh = new(filename, scale, useMaterialColor, smoothNormals, convexPartition)
        return fileMesh
    end
end


text:
using Plots
using Reproduce
using JLD2
using FileIO

function plot_single_autostep(results, log_interval)
    previous_GKSwstype = get(ENV, "GKSwstype", "")
    ENV["GKSwstype"] = "100"

    save_file = "plotting/chunlok/generated_plots/single_autostep_step_sizes.svg"

    p = plot()

    step_sizes = results[:autostep_stepsize]

    # println(findall(x->x!=1.0, step_sizes))
    # dasd
    # println(step_sizes[:, :, 1])
    # println(step_sizes[:, :, 1] - step_sizes[:, :, 100])

    # println(minimum(step_sizes))
    # println(maximum(step_sizes))

    function plot_gvf_step_sizes(range)
        gvf_stepsizes = minimum(step_sizes[range, :, :], dims=(1, 2))[1, 1, :]
        plot!(p, 1:log_interval:size(gvf_stepsizes)[1] * log_interval, gvf_stepsizes)
    end

    plot_gvf_step_sizes(1:4)
    plot_gvf_step_sizes(5:8)
    plot_gvf_step_sizes(9:12)
    plot_gvf_step_sizes(13:16)

    # for i in 1:4
    #     # println(size(error[i, :]))
    #     plot!(p, 1:size(error)[2], error[i, :])
    # end

    # # print(error[3,10])
    savefig(save_file)
    println("Single autostep stepsize plot saved to $(save_file)")

    ENV["GKSwstype"] = previous_GKSwstype
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using Pkg
for p in ("Knet","ArgParse","Images","ImageMagick")
    haskey(Pkg.installed(),p) || Pkg.add(p)
end

"""

julia resnet.jl image-file-or-url

This example implements the ResNet-50, ResNet-101 and ResNet-152 models from
'Deep Residual Learning for Image Regocnition', Kaiming He, Xiangyu Zhang,
Shaoqing Ren, Jian Sun, arXiv technical report 1512.03385, 2015.

* Paper url: https://arxiv.org/abs/1512.03385
* Project page: https://github.com/KaimingHe/deep-residual-networks
* MatConvNet weights used here: http://www.vlfeat.org/matconvnet/pretrained

"""
module ResNet
using Knet, ArgParse
include(Knet.dir("data","imagenet.jl"))

const modelurl = "http://www.vlfeat.org/matconvnet/models/imagenet-resnet-101-dag.mat"
const imgurl = "https://github.com/BVLC/caffe/raw/master/examples/images/cat.jpg"

function main(args)
    s = ArgParseSettings()
    s.description = "resnet.jl (c) Ilker Kesen, 2017. Classifying images with Deep Residual Networks."

    @add_arg_table s begin
        ("image"; default=imgurl; help="Image file or URL.")
        ("--model"; default="imagenet-resnet-101-dag"; help="resnet model name")
        ("--top"; default=5; arg_type=Int; help="Display the top N classes")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array and float type to use")
    end

    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    println(s.description)
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    atype = eval(Meta.parse(o[:atype]))
    model = matconvnet(o[:model])
    avgimg = model["meta"]["normalization"]["averageImage"]
    avgimg = convert(Array{Float32}, avgimg)
    description = model["meta"]["classes"]["description"]
    w, ms = get_params(model["params"], atype)

    @info("Reading $(o[:image])")
    img = imgdata(o[:image], avgimg)
    img = convert(atype, img)

    # get model by length of parameters
    modeldict = Dict(
        162 => (resnet50, "resnet50"),
        314 => (resnet101, "resnet101"),
        467 => (resnet152, "resnet152"))
    !haskey(modeldict, length(w)) && error("wrong resnet MAT file")
    resnet, name = modeldict[length(w)]

    @info("Classifying with $name")
    @time y1 = resnet(w,img,ms)
    z1 = vec(Array(y1))
    s1 = sortperm(z1,rev=true)
    p1 = exp.(logp(z1))
    display(hcat(p1[s1[1:o[:top]]], description[s1[1:o[:top]]]))
    println()
end

# mode, 0=>train, 1=>test
function resnet50(w,x,ms; mode=1)
    # layer 1
    conv1  = conv4(w[1],x; padding=3, stride=2) .+ w[2]
    bn1    = batchnorm(w[3:4],conv1,ms; mode=mode)
    pool1  = pool(bn1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[5:34], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[35:73], r2, ms; mode=mode)
    r4 = reslayerx5(w[74:130], r3, ms; mode=mode) # 5
    r5 = reslayerx5(w[131:160], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[161] * mat(pool5) .+ w[162]
end

# mode, 0=>train, 1=>test
function resnet101(w,x,ms; mode=1)
    # layer 1
    conv1 = reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = pool(conv1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[34:72], r2, ms; mode=mode)
    r4 = reslayerx5(w[73:282], r3, ms; mode=mode)
    r5 = reslayerx5(w[283:312], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[313] * mat(pool5) .+ w[314]
end

# mode, 0=>train, 1=>test
function resnet152(w,x,ms; mode=1)
    # layer 1
    conv1 = reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = pool(conv1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[34:108], r2, ms; mode=mode)
    r4 = reslayerx5(w[109:435], r3, ms; mode=mode)
    r5 = reslayerx5(w[436:465], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[466] * mat(pool5) .+ w[467]
end

# Batch Normalization Layer
# works both for convolutional and fully connected layers
# mode, 0=>train, 1=>test
function batchnorm(w, x, ms; mode=1, epsilon=1e-5)
    mu, sigma = nothing, nothing
    if mode == 0
        d = ndims(x) == 4 ? (1,2,4) : (2,)
        s = prod(size(x,d...))
        mu = sum(x,d) / s
        x0 = x .- mu
        x1 = x0 .* x0
        sigma = sqrt(epsilon + (sum(x1, d)) / s)
    elseif mode == 1
        mu = popfirst!(ms)
        sigma = popfirst!(ms)
    end

    # we need value in backpropagation
    push!(ms, value(mu), value(sigma))
    xhat = (x.-mu) ./ sigma
    return w[1] .* xhat .+ w[2]
end

function reslayerx0(w,x,ms; padding=0, stride=1, mode=1)
    b  = conv4(w[1],x; padding=padding, stride=stride)
    bx = batchnorm(w[2:3],b,ms; mode=mode)
end

function reslayerx1(w,x,ms; padding=0, stride=1, mode=1)
    relu.(reslayerx0(w,x,ms; padding=padding, stride=stride, mode=mode))
end

function reslayerx2(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    ba = reslayerx1(w[1:3],x,ms; padding=pads[1], stride=strides[1], mode=mode)
    bb = reslayerx1(w[4:6],ba,ms; padding=pads[2], stride=strides[2], mode=mode)
    bc = reslayerx0(w[7:9],bb,ms; padding=pads[3], stride=strides[3], mode=mode)
end

function reslayerx3(w,x,ms; pads=[0,0,1,0], strides=[2,2,1,1], mode=1) # 12
    a = reslayerx0(w[1:3],x,ms; stride=strides[1], padding=pads[1], mode=mode)
    b = reslayerx2(w[4:12],x,ms; strides=strides[2:4], pads=pads[2:4], mode=mode)
    relu.(a .+ b)
end

function reslayerx4(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    relu.(x .+ reslayerx2(w,x,ms; pads=pads, strides=strides, mode=mode))
end

function reslayerx5(w,x,ms; strides=[2,2,1,1], mode=1)
    x = reslayerx3(w[1:12],x,ms; strides=strides, mode=mode)
    for k = 13:9:length(w)
        x = reslayerx4(w[k:k+8],x,ms; mode=mode)
    end
    return x
end

function get_params(params, atype)
    len = length(params["value"])
    ws, ms = [], []
    for k = 1:len
        name = params["name"][k]
        value = convert(Array{Float32}, params["value"][k])

        if endswith(name, "moments")
            push!(ms, reshape(value[:,1], (1,1,size(value,1),1)))
            push!(ms, reshape(value[:,2], (1,1,size(value,1),1)))
        elseif startswith(name, "bn")
            push!(ws, reshape(value, (1,1,length(value),1)))
        elseif startswith(name, "fc") && endswith(name, "filter")
            push!(ws, transpose(reshape(value,(size(value,3),size(value,4)))))
        elseif startswith(name, "conv") && endswith(name, "bias")
            push!(ws, reshape(value, (1,1,length(value),1)))
        else
            push!(ws, value)
        end
    end
    map(wi->convert(atype, wi), ws),
    map(mi->convert(atype, mi), ms)
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia vgg.jl cat.jpg
# julia> ResNet.main("cat.jpg")
PROGRAM_FILE=="resnet.jl" && main(ARGS)

end # module


text:
"""
Generate coefficients related to PML boundaries, e.g., damping profiles.
This function outputs a, b, and k variables in eq. 25-26 from Komatitsch 2007 (Geophysics).

	TODO
* At the moment, K=1, and is dummy, unlike in the case of EM when it was useful. Removing K and K_inv might save us memory when using GPUs.
"""
function update_pml!(
    pml::NamedVector{T},
    exmgrid, # 1D grid extended
    mgrid, # 1D grid without extension
    flags::Vector{Bool},
    δt::Float64,
    velavg::Float64,
    freqpeak::Float64, # dominant frequency in Hz
) where {T<:Data.Array}

    nx = length(exmgrid)

    # origin, where the PMLs start (offset b/w the domain of interest and PML region depends on order
    xoriginleft = mgrid[1] - 2 * (step(mgrid))
    xoriginright = mgrid[end] + 2 * (step(mgrid))

    NPOWER = 2.e0
    K_MAX_PML = 1.e0 # from Gedney page 8.11
    ALPHA_MAX_PML = pi * freqpeak # from Festa and Vilotte

    # thickness of the PML layer in meters
    thickness_PML = abs(xoriginleft - exmgrid[1])

    "reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf"
    Rcoef = 0.001e0

    #! check that NPOWER is okay
    #if(NPOWER < 1) stop "NPOWER must be greater than 1"

    # compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
    d0 = -(NPOWER + 1) * (velavg) * log(Rcoef) / (2.e0 * thickness_PML)

    k = zeros(nx)
    fill!(k, 1)
    d = zeros(nx)
    alpha = zero(d)
    a = zero(d)
    b = zero(d)

    "damping in the X direction"
    "origin of the PML layer (position of right edge minus thickness, in meters)"

    for ix = 1:nx
        #---------- left edge
        if (flags[1])
            # define damping profile at the grid points
            abscissa_in_PML = xoriginleft - exmgrid[ix]
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML
                d[ix] = d0 * abscissa_normalized .^ NPOWER
                # this taken from Gedney page 8.2
                k[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized .^ NPOWER
                alpha[ix] =
                    ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.01e0 * ALPHA_MAX_PML
            end
        end
        #---------- right edge
        if (flags[2])
            # define damping profile at the grid points
            abscissa_in_PML = exmgrid[ix] - xoriginright
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML
                d[ix] = d0 * abscissa_normalized^NPOWER
                # this taken from Gedney page 8.2
                k[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized^NPOWER
                alpha[ix] =
                    ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.01e0 * ALPHA_MAX_PML
            end
        end
        #
        # just in case, for -5 at the end
        (alpha[ix] < 0.0) ? alpha[ix] = 0.0 : nothing

        # see equation 25 Komatitsch, 2007, get b and a (need k, alpha, and d before this )
        b[ix] = exp(-(d[ix] / k[ix] + alpha[ix]) * δt)

        # this to avoid division by zero outside the PML
        (abs(d[ix]) > 1.e-6) ?
        a[ix] = d[ix] * (b[ix] - 1.e0) / (k[ix] * (d[ix] + k[ix] * alpha[ix])) : nothing
    end
    pkI = zeros(2 * _fd.npml)
    pa = zeros(2 * _fd.npml)
    pb = zeros(2 * _fd.npml)
    for i = 1:_fd.npml
        j = i + _fd.npml
        pkI[i] = inv(k[i])
        pkI[j] = inv(k[end-_fd.npml+i])
        pa[i] = a[i]
        pa[j] = a[end-_fd.npml+i]
        pb[i] = b[i]
        pb[j] = b[end-_fd.npml+i]
    end
    copyto!(pml[:kI], pkI)
    copyto!(pml[:a], pa)
    copyto!(pml[:b], pb)
end

"""
Generate a NamedArray with PML coefficients for all the dimensions that are then stored in the FDTD structs.
"""
function get_pml(attrib_mod, mgrid)
    dfields = Fields(attrib_mod, "d", ndims = length(mgrid)) # derivative fields that need PML memory
    pnames = [:a, :b, :kI]
    np = 2 * _fd.npml
    return NamedArray(
        [
            NamedArray(Data.Array.([zeros(np), zeros(np), ones(np)]), pnames) for
            df in dfields
        ],
        dfields,
    )
end

"""
Just a loop over dims of pml
"""
function update_pml!(
    pml::NamedVector{T},
    exmgrid,
    mgrid,
    pml_edges::Vector{Symbol},
    attrib_mod,
    args...,
) where {T<:NamedVector}
    for df in names(pml)[1]
        dim = string(last(string(df)))
        i = findfirst(x -> string(x) == dim, dim_names(length(exmgrid)))
        exmgridf = get_mgrid(eval(df)(), attrib_mod, exmgrid...)[i]
        mgridf = mgrid[i]
        update_pml!(
            pml[df],
            exmgridf,
            mgridf,
            [
                any(pml_edges .== Symbol(string(dim), "min")),
                any(pml_edges .== Symbol(string(dim), "max")),
            ],
            args...,
        )
    end
end

function update_pml!(pac)
    update_pml!(
        pac.pml,
        pac.exmedium.mgrid,
        pac.medium.mgrid,
        pac.pml_edges,
        pac.attrib_mod,
        pac.fc[:dt],
        Statistics.mean(pac.exmedium.bounds[:vp]),
        pac.fc[:freqpeak],
    )
end

for dimnames in [zip([:1, :2, :3], dim_names(3)), zip([:1, :2], dim_names(2))]
    is = broadcast(x -> Symbol(string("i", x)), getindex.(collect(dimnames), 2))
    ist = Meta.parse(string("(", [string(s, ",") for s in is]..., ")"))
    N = Meta.parse(string(length(is)))

    for (idim, dim) in dimnames
        i = Symbol("i", string(dim))
        ismoff = replace(is, i => :($i + moff))
        isdoff = replace(is, i => :(doff + $i))
        for (fname, fnamenp, imoff) in zip(
            [Symbol("memory", string(dim), "!"), Symbol("memory1", string(dim), "!")],
            [Symbol("memorynp", string(dim), "!"), Symbol("memorynp1", string(dim), "!")],
            [:($i + moff), :($i + moff + 1)],
        )
            @eval @parallel_indices(
                $ist,
                function $fnamenp(memory::Data.Array{$N}, d, a, b, kI, moff, doff)
                    memory[$(ismoff...)] =
                        b[$imoff] * memory[$(ismoff...)] + a[$imoff] * d[$(isdoff...)]
                    d[$(isdoff...)] = d[$(isdoff...)] * kI[$imoff] + memory[$(ismoff...)]
                    return
                end
            )
            @eval function $fname(memory::Data.Array{$N}, d, a, b, kI)
                sm = collect(size(memory))
                setindex!(sm, _fd.npml, $idim)
                # first _fd.npml points
                @parallel map(x -> (:)(1, x), Tuple(sm)) $fnamenp(memory, d, a, b, kI, 0, 0)
                # last _fd.npml points independent of d
                @parallel map(x -> (:)(1, x), Tuple(sm)) $fnamenp(
                    memory,
                    d,
                    a,
                    b,
                    kI,
                    _fd.npml,
                    getindex(size(d), $idim) - _fd.npml,
                )
            end
        end
    end
end