text:
"""
    function get_flightdata(filename)

Method for collecting flight data from relevant .csv file
"""
function get_flightdata(filename)
    flights=readtable(filename)
    into_mombassa=convert(Vector,flights[:,1])
    into_nairobi=convert(Vector,flights[:,2])
    return into_mombassa,into_nairobi
end
"""
    function get_prevdata(filename)

Method for collecting global prevalence from relevant .csv file
"""
function get_prevdata(filename)
    prev_table=readtable(filename)
    global_prev=convert(Vector,prev_table[:,1])
    return global_prev
end

function transportstructure_params!(P::CoVParameters,ρ::Vector{Float64},transport_matrix)
    d1,d2 = size(transport_matrix)
   #Put in the correct location matrix
   for i = 1:d1,j = 1:d2
       if i != j
           P.T[i,j] = P.ρ[j]*transport_matrix[i,j]
       else
           P.T[i,j] = 1-P.ρ[j]
       end
   end
   return nothing
end

function transportstructure_params!(P::CoVParameters_AS,ρ::Vector{Float64},transport_matrix)
    d1,d2 = size(transport_matrix)
    P.ρ = ρ
   #Put in the correct location matrix
   for i = 1:d1,j = 1:d2
       if i != j
           P.T[i,j] = P.ρ[j]*transport_matrix[i,j]
       else
           P.T[i,j] = 1-P.ρ[j]
       end
   end
   return nothing
end

function transportstructure_params!(P::CoVParameters,ρ::Float64,transport_matrix)
    P.ρ = ρ
    #Put in the correct location matrix
    for i = 1:n,j = 1:n
        if i != j
            P.T[i,j] = P.ρ*transport_matrix[i,j]
        else
            P.T[i,j] = 1-P.ρ
        end
    end
    return nothing
end

function model_ingredients_from_data(datatablename,mixingmatrixname,travelmatrixname,flight_filename,prev_filename)
    KenyaTbl = readtable(datatablename)
    n_data, = size(KenyaTbl)
    @load mixingmatrixname T_opt ρ_county
    @load travelmatrixname P_opt
    if n_data != n
        println("ERROR: difference between number of areas in data, and definition.")
    end
    #Population state array --- read number of areas from T_opt from
    d1,d2 = size(T_opt)
    suspop_kenya = zeros(Int64,d1,n_s,2) #Array by area, state and urban vs rural
    for i = 1:d1
        suspop_kenya[i,1,1] = KenyaTbl.Urban[i]
        if !ismissing(KenyaTbl.Rural[i])
            suspop_kenya[i,1,2] = KenyaTbl.Rural[i]
        end
    end
    N_urb = [sum(suspop_kenya[i,:,1]) for i = 1:d1]
    N_rural = [sum(suspop_kenya[i,:,2]) for i = 1:d1]
    N̂ = T_opt*N_urb + N_rural

    global_prev = get_prevdata(prev_filename)
    into_mom, into_nai = get_flightdata(flight_filename)
    change_matrix(dc,1,1,1,1)
    # Now get flight numbers and global prevalence
    #Parameter definition
    P = CoVParameters(T = T_opt,ρ = ρ_county,
                        into_mom = into_mom, into_nai = into_nai,
                        global_prev = global_prev,
                        N̂=N̂,dc = sparse(dc))
    #Matrix for in-place tau-leap calculations
    return suspop_kenya,P,P_opt

end

function model_ingredients_from_data(agestructuredata_filename,flight_filename,prev_filename)
    @load agestructuredata_filename N_region_age M_Kenya movements_per_person P_dest ρ T σ rel_detection_rates M_Kenya_ho hosp_rate_by_age ICU_rate_by_age_cond_hosp
    n_wa,n_a = size(N_region_age)
    #Population state array
    suspop_kenya = zeros(Int64,n_wa,n_a,n_s) #Array by area, age group and disease state
    for i = 1:n_wa,j=1:n_a
        suspop_kenya[i,j,1] = N_region_age[i,j]
    end

    #Estimate effective population size in each area after mobility
    # N_mobile = N_region_age[:,mobile_age_indices]
    # N_immobile = N_region_age[:,immobile_age_indices]
    N̂ = T*N_region_age
    N̂[immobile_age_indices] .= N_region_age[immobile_age_indices]

    # Now get flight numbers and global prevalence
    global_prev = get_prevdata(prev_filename)
    into_mom, into_nai = get_flightdata(flight_filename)

    #Define the change matrix
    dc = sparse(zeros(Int64,n_wa*n_a*n_s,n_wa*n_a*n_ta))
    change_matrix(dc)

    #Parameter definition
    P = CoVParameters_AS(T = T,ρ = ρ,χ = σ,rel_detection_rate = rel_detection_rates[:,2], #This assumes that the detection rates match ϵ = 0.1
                        into_mom = into_mom, into_nai = into_nai,
                        global_prev = global_prev,
                        M = M_Kenya,
                        M_ho = M_Kenya_ho,
                        hₐ = hosp_rate_by_age,
                        ICUₐ = ICU_rate_by_age_cond_hosp,
                        dc=dc,
                        N̂=N̂)
    #Matrix for in-place tau-leap calculations
    return suspop_kenya,P,P_dest

end

function create_KenyaCoV_prob(u0,tspan,P::CoVParameters)
    u0_vec = u0[:]
    prob_tl = DiscreteProblem(u0_vec,tspan,P)
    return JumpProblem(prob_tl,Direct(),reg_jumps_forKenyaCoV)
end

function solve_KenyaCoV_prob(u0,tspan,P::CoVParameters,dt=1.)
    jump_prob_tl = create_KenyaCoV_prob(u0,tspan,P)
    return solve(jump_prob_tl,SimpleTauLeaping();dt = dt)
end

function create_KenyaCoV_non_neg_prob(u0,tspan,P::CoVParameters)
    u0_vec = u0[:]
    return prob = DiscreteProblem(nonneg_tauleap,u0_vec,tspan,P)
end

function create_KenyaCoV_non_neg_prob(u0,tspan,P::CoVParameters_AS)
    return prob = DiscreteProblem(nonneg_tauleap,u0,tspan,P)
end

function create_KenyaCoV_ode_prob(u0,tspan,P::CoVParameters)
    _u0 = convert.(Float64,u0)
    return prob = ODEProblem(ode_model,_u0,tspan,P)
end

function create_KenyaCoV_ode_prob(u0,tspan,P::CoVParameters_AS)
    _u0 = convert.(Float64,u0)
    return prob = ODEProblem(ode_model,_u0,tspan,P)
end

function calculate_R₀_scale(P::CoVParameters_AS)
    sus_matrix = repeat(P.χ,1,KenyaCoV.n_a)
    trans_matrix = repeat(P.δ*P.rel_detection_rate' + P.ϵ*(1 .- P.δ*P.rel_detection_rate'),KenyaCoV.n_a,1)
    eigs, = eigen(sus_matrix.*M_China.*trans_matrix)
    return Real(eigs[end])
end
function calculate_R₀(P::CoVParameters_AS)
    sus_matrix = repeat(P.χ,1,KenyaCoV.n_a)
    trans_matrix = repeat(P.δ*P.rel_detection_rate' + P.ϵ*(1 .- P.δ*P.rel_detection_rate'),KenyaCoV.n_a,1)
    K = sus_matrix.*P.M.*trans_matrix*P.β/P.γ
    v₀ = zeros(n_a)
    v₀[6] = 1
    eigs,evects = eigen(K)
    return Real(eigs[end]), LinearAlgebra.normalize!(abs.(evects[:,end]),1),(K^5)*v₀
end

function calculate_R₀_homeonly(P::CoVParameters_AS)
    sus_matrix = repeat(P.χ,1,KenyaCoV.n_a)
    trans_matrix = repeat(P.δ*P.rel_detection_rate' + P.ϵ*(1 .- P.δ*P.rel_detection_rate'),KenyaCoV.n_a,1)
    K = sus_matrix.*P.M.*trans_matrix*P.β/P.γ
    v₀ = zeros(n_a)
    v₀[6] = 1
    eigs,evects = eigen(K)
    return Real(eigs[end]), LinearAlgebra.normalize!(abs.(evects[:,end]),1),(K^2)*v₀
end


text:
@doc raw"""
    PolyaGammaPoisson(y::Real, c::Real, λ::Real)

A bivariate distribution, used as hierachical prior as:
```math
    p(\omega, n) = \operatorname{PG}(\omega|y + n, c)\operatorname{Po}(n|\lambda).
```

Random samples as well as statistics from the distribution will returned as `NamedTuple` : `(;ω, n)`.

This structured distributions is needed for example for the [`PoissonLikelihood`](https://juliagaussianprocesses.github.io/GPLikelihoods.jl/dev/#GPLikelihoods.PoissonLikelihood).
"""
struct PolyaGammaPoisson{Ty,Tc,Tλ} <: AbstractNTDist
    y::Ty # Intermediate first parameter for PG(y + n, c)
    c::Tc # Second parameter for PG
    λ::Tλ # Poisson Parameter
end

Distributions.Poisson(d::PolyaGammaPoisson) = Poisson(d.λ)

Distributions.length(::PolyaGammaPoisson) = 2

function ntrand(rng::AbstractRNG, d::PolyaGammaPoisson)
    n = rand(rng, Poisson(d))
    ω = rand(rng, PolyaGamma.(n + d.y, d.c))
    return (; ω, n)
end

function MeasureBase.logdensity_def(d::PolyaGammaPoisson, x::NamedTuple)
    logpdf_n = logpdf(Poisson(d), x.n)
    logpdf_ω = logpdf(PolyaGamma(d.y + x.n, d.c), x.ω)
    return logpdf_ω + logpdf_n
end

function tvmean(ds::AbstractVector{<:PolyaGammaPoisson})
    n = getfield.(ds, :λ)
    ω = map(ds, n) do d, n
        mean(PolyaGamma(d.y + n, d.c))
    end
    return TupleVector((; ω, n))
end

function ntmean(d::PolyaGammaPoisson)
    return (; ω=mean(PolyaGamma(d.y + d.λ, d.c)), n=d.λ)
end

function Distributions.kldivergence(q::PolyaGammaPoisson, p::PolyaGammaPoisson)
    (p.c == 0 && p.y == q.y) || error("No solution for this prior. qΩ = $q, pΩ = $p")
    return kldivergence(PolyaGamma(q.y + q.λ, q.c), PolyaGamma(q.y + q.λ, 0)) +
           kldivergence(Poisson(q), Poisson(p))
end

ext: jl


text:
function update_weights(data, e, closest_same, closest_other, weights, weights_mult, m, k, max_f_vals,
						min_f_vals, dm_vals_same, dm_vals_other, features_msk_same, features_msk_other)

	for t = 1:size(data, 2)

		# Penalty term
		penalty = ((abs.(e[t] .- closest_same[:, t])/(max_f_vals[t] .- min_f_vals[t] .+ eps(Float64))) -
				   dm_vals_same[:,t])[features_msk_same[:, t]]
		penalty = sum(penalty)

		# Reward term
		reward = ((weights_mult .* (abs.(e[t] .- closest_other[:, t])/(max_f_vals[t] .- min_f_vals[t] .+ eps(Float64)))) -
				  dm_vals_other[:,t])[features_msk_other[:, t]]
		reward = sum(reward)

		# Weights update
		weights[t] = weights[t] - penalty/(m*k) + reward/(m*k)
	end

	# Return updated weights.
	return vec(weights)
end

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
isdefined(Base, :__precompile__) && __precompile__()

module CUDAdrv

using Compat
import Compat.String

include("util/logging.jl")

include("pointer.jl")

include("errors.jl")
include("base.jl")
include("version.jl")
include("devices.jl")
include("context.jl")
include("module.jl")
include("memory.jl")
include("stream.jl")
include("execution.jl")
include("events.jl")
include("profile.jl")

include("array.jl")

function __init__()
    __init_logging__()
    __init_library__()
end

end


text:
#=

amss_figures.jl

@author: Shunsuke Hori

=#

include("CES.jl")
include("BGP.jl")
using CES
using BGP

#initialize mugrid for value function iteration
muvec = linspace(-0.7,0.01,200)

#Time Varying Example

CES.M_time_example.transfers = true #Government can use transfers
PP_seq_time = LS.Planners_Allocation_Sequential(CES.M_time_example) #solve sequential problem

PP_im_time = Planners_Allocation_Bellman(CES.M_time_example,muvec)

sHist_h = [1,2,3,4,6,6,6]
sHist_l = [1,2,3,5,6,6,6]

sim_seq_h = LS.simulate(PP_seq_time,1.,1,7,sHist_h)
sim_im_h = simulate(PP_im_time,1.,1,7,sHist_h)
sim_seq_l = LS.simulate(PP_seq_time,1.,1,7,sHist_l)
sim_im_l = simulate(PP_im_time,1.,1,7,sHist_l)

using Plots
pyplot()
p=plot(size = (700,500),layout =(3,2),
        xaxis=(0:6), grid=false, titlefont=Plots.font("sans-serif", 10))
title!(p[1],"Consumption")
plot!(p[1],0:6,sim_seq_l[1],marker=:circle,color=:black,legend=false)
plot!(p[1],0:6,sim_im_l[1],marker=:circle,color=:red,legend=false)
plot!(p[1],0:6,sim_seq_h[1],marker=:utriangle,color=:black,legend=false)
plot!(p[1],0:6,sim_im_h[1],marker=:utriangle,color=:red,legend=false)
title!(p[2],"Labor")
plot!(p[2],0:6,sim_seq_l[2],marker=:circle,color=:black,lab="Complete Markets")
plot!(p[2],0:6,sim_im_l[2],marker=:circle,color=:red,lab="Incomplete Markets")
plot!(p[2],0:6,sim_seq_h[2],marker=:utriangle,color=:black,lab="")
plot!(p[2],0:6,sim_im_h[2],marker=:utriangle,color=:red,lab="")
title!(p[3],"Government Debt")
plot!(p[3],0:6,sim_seq_l[3],marker=:circle,color=:black,legend=false)
plot!(p[3],0:6,sim_im_l[3],marker=:circle,color=:red,legend=false)
plot!(p[3],0:6,sim_seq_h[3],marker=:utriangle,color=:black,legend=false)
plot!(p[3],0:6,sim_im_h[3],marker=:utriangle,color=:red,legend=false)
title!(p[4],"Tax Rate")
plot!(p[4],0:6,sim_seq_l[4],marker=:circle,color=:black,legend=false)
plot!(p[4],0:6,sim_im_l[5],marker=:circle,color=:red,legend=false)
plot!(p[4],0:6,sim_seq_h[4],marker=:utriangle,color=:black,legend=false)
plot!(p[4],0:6,sim_im_h[5],marker=:utriangle,color=:red,legend=false)
title!(p[5],"Government Spending", ylims=(0.05,0.25))
plot!(p[5],0:6,CES.M_time_example.G[sHist_l],marker=:circle,color=:black,legend=false)
plot!(p[5],0:6,CES.M_time_example.G[sHist_l],marker=:circle,color=:red,legend=false)
plot!(p[5],0:6,CES.M_time_example.G[sHist_h],marker=:utriangle,color=:black,legend=false)
plot!(p[5],0:6,CES.M_time_example.G[sHist_h],marker=:utriangle,color=:red,legend=false)
title!(p[6],"Output")
plot!(p[6],0:6,CES.M_time_example.Theta[sHist_l].*sim_seq_l[2],marker=:circle,color=:black,legend=false)
plot!(p[6],0:6,CES.M_time_example.Theta[sHist_l].*sim_im_l[2],marker=:circle,color=:red,legend=false)
plot!(p[6],0:6,CES.M_time_example.Theta[sHist_h].*sim_seq_h[2],marker=:utriangle,color=:black,legend=false)
plot!(p[6],0:6,CES.M_time_example.Theta[sHist_h].*sim_im_h[2],marker=:utriangle,color=:red,legend=false)
savefig("TaxSequence_time_varying_AMSS.png")
p1=p

#=
BGP Example
=#
BGP.M1.transfers = false   #Government can use transfers
PP_seq = LS.Planners_Allocation_Sequential(BGP.M1) #solve sequential problem
PP_bel = LS.Planners_Allocation_Bellman(BGP.M1,muvec) #solve recursive problem
PP_im = Planners_Allocation_Bellman(BGP.M1,muvec)

T = 20
#sHist = utilities.simulate_markov(M1.Pi,0,T)
sHist = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1]

#simulate
sim_seq = LS.simulate(PP_seq,0.5,1,T,sHist)
#sim_bel = PP_bel.simulate(0.5,0,T,sHist)
sim_im = simulate(PP_im,0.5,1,T,sHist)

#plot policies
p=plot(size = (700,500), layout = grid(3,2),
        xaxis=(0:T), grid=false, titlefont=Plots.font("sans-serif", 10))
title!(p[1], "Consumption")
plot!(p[1], sim_seq[1], marker=:circle, color=:black, lab="Complete Markets")
#plt[:plot](sim_bel[0],"-xk")
plot!(p[1], sim_im[1], marker=:utriangle, color=:black, lab="Incomplete Markets")
title!(p[2], "Labor")
plot!(p[2], sim_seq[1], marker=:circle, color=:black, legend=false)
#plt.plot(sim_bel[1],"-xk")
plot!(p[2], sim_im[1], marker=:utriangle, color=:black, legend=false)
title!(p[3], "Government Debt")
plot!(p[3], sim_seq[2], marker=:circle, color=:black, legend=false)
#plt.plot(sim_bel[2],"-xk")
plot!(p[3], sim_im[2], marker=:utriangle, color=:black, legend=false)
title!(p[4], "Tax Rate")
plot!(p[4], sim_seq[4], marker=:circle, color=:black, legend=false)
#plt.plot(sim_bel[3],"-xk")
plot!(p[4], sim_im[5], marker=:utriangle, color=:black, legend=false)
title!(p[5], "Government Spending", ylims=(0.05,0.25))
plot!(p[5], BGP.M1.G[sHist], marker=:circle, color=:black, legend=false)
#plt.plot(M1.G[sHist],"-^k")
title!(p[6], "Output")
plot!(p[6], BGP.M1.Theta[sHist].*sim_seq[2], marker=:circle,
        color=:black, legend=false)
#plt.plot(M1.Theta[sHist]*sim_bel[1],"-xk")
plot!(p[6], BGP.M1.Theta[sHist].*sim_im[2], marker=:utriangle,
        color=:black, legend=false)
savefig("TaxSequence_AMSS.png")
p2=p

#Now long simulations
T_long = 200
sim_seq_long = LS.simulate(PP_seq, 0.5,1,T_long)
sHist_long = sim_seq_long[end-2]
sim_im_long = simulate(PP_im, 0.5,1,T_long,sHist_long)

p=plot(size = (700,500), layout = (3,2), xaxis=(0:50:T_long), grid=false,
        titlefont=Plots.font("sans-serif", 10))
title!(p[1],"Consumption")
plot!(p[1], sim_seq_long[1], color=:black, linestyle=:solid, lab="Complete Markets")
plot!(p[1], sim_im_long[1], color=:black, linestyle=:dot, lab="Incomplete Markets")
title!(p[2],"Labor")
plot!(p[2], sim_seq_long[1], color=:black, linestyle=:solid, legend=false)
plot!(p[2], sim_im_long[1], color=:black, linestyle=:dot, legend=false)
title!(p[3],"Government Debt")
plot!(p[3], sim_seq_long[2], color=:black, linestyle=:solid, legend=false)
plot!(p[3], sim_im_long[2], color=:black, linestyle=:dot, legend=false)
title!(p[4],"Tax Rate")
plot!(p[4], sim_seq_long[3], color=:black, linestyle=:solid, legend=false)
plot!(p[4], sim_im_long[4], color=:black, linestyle=:dot, legend=false)
title!(p[5],"Government Spending",ylims=(0.05,0.25))
plot!(p[5], BGP.M1.G[sHist_long], color=:black, linestyle=:solid, legend=false)
plot!(p[5], BGP.M1.G[sHist_long], color=:black, linestyle=:dot, legend=false)
title!(p[6],"Output")
plot!(p[6],BGP.M1.Theta[sHist_long].*sim_seq_long[2],
    color=:black, linestyle=:solid, legend=false)
plot!(p[6],BGP.M1.Theta[sHist_long].*sim_im_long[2],
    color=:black, linestyle=:dot, legend=false)
# plt[:tight_layout]()
savefig("Long_SimulationAMSS.png")
p3=p

#=
Show Convergence example
=#
M_convergence = CES.M1
muvec = linspace(-0.15,0.0,100) #change
PP_C = Planners_Allocation_Bellman(M_convergence,muvec)
xgrid = PP_C.xgrid
xf = PP_C.policies[end-1] #get x policies
p= plot()
for s in 1:2
    plot!(p,xgrid,xf[1,s].(xgrid)-xgrid)
end
p4=p

sim_seq_convergence = simulate(PP_C,0.5,1,2000)
sHist_long = sim_seq_convergence[end]

p=plot(size = (700,500),layout = grid(3,2), xaxis=(0:200:2000),
        grid=false, titlefont=Plots.font("sans-serif", 10))
title!("Consumption")
plot!(sim_seq_convergence[1],color=:black,
    lab=["Complete Markets" "Incomplete Markets"])
title!(p[2],"Labor")
plot!(p[2],sim_seq_convergence[2],color=:black,legend=false)
title!(p[3],"Government Debt")
plot!(p[3],sim_seq_convergence[3],color=:black,legend=false)
title!(p[4],"Tax Rate")
plot!(p[4],sim_seq_convergence[4],color=:black,legend=false)
title!(p[5],"Government Spending",ylims=(0.05,0.25))
plot!(p[5],M_convergence.G[sHist_long],color=:black,legend=false)
title!(p[6],"Output")
plot!(p[6],M_convergence.Theta[sHist_long].*sim_seq_convergence[2],color=:black,legend=false)
savefig("Convergence_SimulationAMSS.png")
p5=p


text:
module Printer
    export printer_error

    function printer_error(s)
        totals = length(s)
        errors = count( c -> c < 'a' || 'm' < c ,s )
        "$errors/$totals"
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
include("nanop.jl")
include("CESMReader.jl")

using Statistics
using NCDatasets

using ArgParse
using JSON
using Formatting

using .CESMReader

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--ilev-file"
            help = "File of pressure levels. Variable assume to be `ilev`."
            arg_type = String
            default = "ilev"

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--varname"
            help = "Variable name of streamfunction."
            arg_type = String
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

function findHadleyCellsStrength(ψ, lat, lev; hemisphere)

    ψ_flat = reshape(ψ, reduce(*, size(ψ)))

    if hemisphere == :NH
        ψ_m, ψ_m_idx = findmax(ψ_flat)
    elseif hemisphere == :SH
        ψ_m, ψ_m_idx = findmin(ψ_flat)
    else
        throw(ErrorException("Unknown keyword: " * string(hemisphere)))
    end

    lat_m_idx = mod(ψ_m_idx - 1, length(lat)) + 1
    lev_m_idx = floor(Int64, (ψ_m_idx - 1) / length(lat)) + 1

    #println(AMOC_max_idx, " => ( ",  lat_max_idx, ", ", z_max_idx, " )")

    lat_m = lat[lat_m_idx]
    lev_m = lev[lev_m_idx]

    if hemisphere == :NH
        if ! ( 0 < lat_m <= 30 && 250 <= lev_m <= 750  )
            println(format("lat_m: {:f}, lev_m: {:f}", lat_m, lev_m))
            throw(ErrorException("ψ maximum not within desired range."))
        end
    elseif hemisphere == :SH
        if ! ( -30 <= lat_m < 0 && 250 <= lev_m <= 750  )
            println(format("lat_m: {:f}, lev_m: {:f}", lat_m, lev_m))
            throw(ErrorException("ψ maximum not within desired range."))
        end
    end

    return ψ_m, lat_m, lev_m
end

Dataset(parsed["domain-file"], "r") do ds
    global mask = replace(ds["mask"], missing=>NaN)
    global lat  = replace(ds["yc"][1, :], missing=>NaN)
    global Ny = length(lat)
end

Dataset(parsed["ilev-file"], "r") do ds
    global ilev = replace(ds["ilev"], missing=>NaN)
    global Nz = length(ilev)
end

if parsed["data-file-timestamp-form"] == "YEAR"

    filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
    form = :YEAR
elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
    filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
    form = :YEAR_MONTH
end

global beg_t = (parsed["beg-year"] - 1) * 12 + 1
global end_t = (parsed["end-year"] - 1) * 12 + 12
global Nt = end_t - beg_t + 1

fh = FileHandler(filename_format=filename_format, form=form)
ψ = reshape( getData(fh, parsed["varname"], (parsed["beg-year"], parsed["end-year"]), (:, :)), Ny, Nz, Nt)

#Dataset(parsed["data-file"], "r") do ds
#    global ψ = reshape(replace(ds[parsed["varname"]][:], missing=>NaN), Ny, Nz, :)
#    global Nt = size(ψ)[3]
#end

years = Int64(Nt/12)

hadleycell_north_max = zeros(Float64, years)
hadleycell_north_max_lat = zeros(Float64, years)
hadleycell_north_max_lev = zeros(Float64, years)

hadleycell_south_max = zeros(Float64, years)
hadleycell_south_max_lat = zeros(Float64, years)
hadleycell_south_max_lev = zeros(Float64, years)

for y=1:years
    _ψ = mean( view(ψ, :, :, ((y-1)*12+1):(y*12)), dims=3 )[:, :, 1]
    hadleycell_north_max[y], hadleycell_north_max_lat[y], hadleycell_north_max_lev[y] = findHadleyCellsStrength(_ψ, lat, ilev; hemisphere=:NH)
    hadleycell_south_max[y], hadleycell_south_max_lat[y], hadleycell_south_max_lev[y] = findHadleyCellsStrength(_ψ, lat, ilev; hemisphere=:SH)
end

println(hadleycell_south_max)

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in  [

        ("hadleycell_north_max",     hadleycell_north_max,     ("time",), Dict()),
        ("hadleycell_north_max_lat", hadleycell_north_max_lat, ("time",), Dict()),
        ("hadleycell_north_max_lev", hadleycell_north_max_lev, ("time",), Dict()),

        ("hadleycell_south_max",     hadleycell_south_max,     ("time",), Dict()),
        ("hadleycell_south_max_lat", hadleycell_south_max_lat, ("time",), Dict()),
        ("hadleycell_south_max_lev", hadleycell_south_max_lev, ("time",), Dict()),

    ]
        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        println("Writing variable:  ", varname)
        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
let
    lrec = ListRecord(1.0, Float64, Bool, Int)
    append!(lrec.frames, [RecordFrame(1,2), RecordFrame(3,4)])
    append!(lrec.states, [RecordState(1.0,1),RecordState(2.0,2),RecordState(3.0,1),RecordState(4.0,3)])
    lrec.defs[1] = true
    lrec.defs[2] = true
    lrec.defs[3] = false

    sparsemat, id_lookup = get_sparse_lookup(lrec)
    @test sparsemat[1,1] == 1.0
    @test sparsemat[2,1] == 3.0
    @test sparsemat[1,2] == 2.0
    @test sparsemat[2,2] == 0.0
    @test sparsemat[1,3] == 0.0
    @test sparsemat[2,3] == 4.0
    @test id_lookup == Dict(2=>2,3=>3,1=>1)

    qrec = convert(QueueRecord{Entity{Float64, Bool, Int}}, lrec)
    @test qrec[-1][1] == Entity(1.0,true,1)
    @test qrec[-1][2] == Entity(2.0,true,2)
    @test qrec[ 0][1] == Entity(3.0,true,1)
    @test qrec[ 0][2] == Entity(4.0,false,3)
    @test nframes(qrec) == 2
end

let
    qrec = QueueRecord(Entity{Float64, Bool, Int}, 2, 0.1)
    @test nframes(qrec) == 0

    update!(qrec, Frame([Entity(1.0,true,1), Entity(2.0,true,2)]))
    update!(qrec, Frame([Entity(3.0,true,1), Entity(4.0,false,3)]))

    @test nframes(qrec) == 2

    lrec = convert(ListRecord{Float64, Bool, Int}, qrec)
    @test nframes(lrec) == 2
    @test lrec.frames == [RecordFrame(1,2), RecordFrame(3,4)]
    @test lrec.states == [RecordState(1.0,1),RecordState(2.0,2),RecordState(3.0,1),RecordState(4.0,3)]
    @test lrec.defs[1] == true
    @test lrec.defs[2] == true
    @test lrec.defs[3] == false
    @test length(lrec.defs) == 3
end


text:
# using Pkg; Pkg.activate("../..")

# Implement routing with a FVM
using Test, Plots, Printf, Statistics

# physics scales
const ρ̂i  = 910
const ρ̂w  = 1000
const ĝ   = 9.81
const ϕ̂   = 1e7
const x̂   = 1e6
const M̂   = 1e-6
const D̂   = 1
# physics dependent
const t̂   = D̂/M̂
const Ĥ   = ϕ̂/ρ̂w/ĝ
const Ψ̂   = ϕ̂/x̂
const û   = x̂/t̂
const k̂   = û/Ψ̂
const rr  = ρ̂i/ρ̂w
const rdh = D̂/Ĥ    # fixed ice thikness routing
const s2d = 3600*24

@views function flow_routing1D()
# nondim
Lx     = 10e3/x̂
k      = 0.1
ttot   = 1000e10*s2d/t̂ # ttot>>1 && nt=1 steady state
M      = M̂*s2d
tim_p  = 0.0
# numerics
nx     = 100
nt     = 1
nout   = 1000
ε      = 1e-10     # aboslute tolerance
damp   = 0.9
dτ_sc  = 12.0
cn     = 0.0       # crank-Nicolson if cn=0.5 (transient)
dx     = Lx/nx
xc     = dx/2:dx:Lx-dx/2
dt     = ttot/nt
# initial
Zb     = (Lx.-xc)./Lx./100.0 .- 30.0./Ĥ.*exp.(-(xc.-Lx/2.0).^2/(Lx/8.0)^2)
H      = 100.0./Ĥ .- xc.*3.0
# H      = (Lx-xc)/Lx/100 + 30/H_s*exp(-(xc-Lx/2).^2/(Lx/8)^2) + (xc + 50)/H_s;
D      = 1.0.*ones(nx)
# D      = max(Zb)-Zb + 10;
∂D     = zeros(nx)
D_o    = zeros(nx)
∂D_o   = zeros(nx)
ϕ      = zeros(nx)
U      = zeros(nx-1)
qD     = zeros(nx-1)
RD     = zeros(nx)
∂Ddτ   = zeros(nx)
errD   = zeros(nx)
err1=[]; err2=[]
# action
for it = 1:nt
    D_o  .= D
    ∂D_o .= ∂D
    iter = 1; err = 2*ε
    while err > ε
        errD .= D
        # flow routing physics
        # if rout==1
        ϕ           .= Zb .+ rdh.*D .+ rr.*H              # fixed ice thikness routing
        # elseif rout==2
        # Zs          = Zb + H;                           # fixed ice surface routing
        # phi         = Zs + (1 - rr)*Zb + rdz*(1-rr)*D   # fixed ice surface routing
        # end
        U           .= -k*diff(ϕ)./dx
        qD          .= U.*(U.>0.0).*D[1:end-1] .+ U.*(U.<0.0).*D[2:end] # not working without upwind: qD = U.*(D(1:end-1)-D(2:end)).*0.5;
        ∂D[2:end-1] .= .-diff(qD)./dx .+ 1.0/M
        RD          .= .-(D.-D_o)./dt .+ ((1-cn).*∂D + cn*∂D_o)         # transient included if dt<inf
        # solver
        ∂Ddτ        .= damp*∂Ddτ .+ RD                                  # damping the rate of change
        dτ           = 1.0/(1.0/dt + 1.0/(dx/maximum(U)/dτ_sc))         # iterative timestep
        D           .= D .+ dτ.*∂Ddτ                                    # update D
        D[[1,end]]  .= 0.0 # BC
        # Check errs
        if mod(iter,nout)==0
            errD .= D.-errD
            ermb  = 1.0./M*dx*(nx-2) - sum(abs.(qD[[1, end]]))
            push!(err1, maximum(abs.(errD))); push!(err2, abs(ermb))
            err   = err1[end]
            @printf("iter=%d  errD=%1.3e, errMB=%1.3e \n", iter, err1[end], err2[end])
        end
        iter+=1
    end
    tim_p = tim_p + dt
    # ploting
    default(size=(700,800))
    p1 = plot(xc*x̂/1e3, Zb*Ĥ, label="Zb", linewidth=2)
        plot!(xc*x̂/1e3, D*D̂+(H+Zb)*Ĥ, label="D+Zb+H", linewidth=2, title=string("Time = ",tim_p*t̂/s2d," days"), framestyle=:box)
    p2 = plot(xc*x̂/1e3, ϕ, ylabel="ϕ", linewidth=2, framestyle=:box, legend=false)
    p3 = plot(xc[2:end-1]*x̂/1e3, D[2:end-1]*D̂, ylabel="D [m]", yscale=:log10, linewidth=2, framestyle=:box, legend=false)
    p4 = plot(xc*x̂/1e3, D*D̂ + Zb*Ĥ, label="D+Zb [m]", linewidth=2)
        plot!(xc*x̂/1e3, Zb*Ĥ, label="Zb [m]", linewidth=2, framestyle=:box)
    p5 = plot(xc[2:end]*x̂/1e3, U*û, ylabel="U", linewidth=2, framestyle=:box, legend=false)
    p6 = plot(xc[2:end]*x̂/1e3, qD, xlabel="x [km]", ylabel="q", linewidth=2, framestyle=:box, legend=false)
    display(plot(p1, p2, p3, p4, p5, p6, layout=(6,1)))
    # figure(2),clf
    # semilogy(err1,'Linewidth',LW),hold on
    # semilogy(err2,'Linewidth',LW),hold off,legend('errD','errMB')
end

return
end

@time flow_routing1D()


text:
using Clustering
using Distances: SqEuclidean, pairwise
using Test

@testset "seeding" begin

Random.seed!(34568)

@test RandSeedAlg <: SeedingAlgorithm
@test KmppAlg <: SeedingAlgorithm
@test KmCentralityAlg <: SeedingAlgorithm

alldistinct(x::Vector{Int}) = (length(Set(x)) == length(x))

function min_interdist(X::AbstractMatrix)
    dists = pairwise(SqEuclidean(), X, dims=2)
    n = size(X, 2)
    r = Inf
    for i = 1:n, j = 1:n
        if i != j && dists[i,j] < r
            r = dists[i,j]
        end
    end
    return r
end

d = 3
n = 100
k = 5
X = rand(d, n)
C = pairwise(SqEuclidean(), X, dims=2)

Xt = copy(transpose(X))
Ct = copy(transpose(C))

md0 = min_interdist(X)

@testset "RandSeed" begin
    Random.seed!(34568)
    iseeds = initseeds(RandSeedAlg(), X, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds(RandSeedAlg(), Xt', k)
    @test iseeds == iseeds_t

    Random.seed!(34568)
    iseeds = initseeds_by_costs(RandSeedAlg(), C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds_by_costs(RandSeedAlg(), Ct', k)
    @test iseeds == iseeds_t

    R = copyseeds(X, iseeds)
    @test isa(R, Matrix{Float64})
    @test R == X[:, iseeds]
    R_t = copyseeds(Xt', iseeds)
    @test R == R_t
end

@testset "Kmpp" begin
    Random.seed!(34568)
    iseeds = initseeds(KmppAlg(), X, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds(KmppAlg(), Xt', k)
    @test iseeds == iseeds_t

    Random.seed!(34568)
    iseeds = initseeds_by_costs(KmppAlg(), C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds_by_costs(KmppAlg(), Ct', k)
    @test iseeds == iseeds_t

    @test min_interdist(X[:, iseeds]) > 20 * md0
    @test min_interdist((Xt')[:, iseeds]) > 20 * md0

    Random.seed!(34568)
    iseeds = kmpp(X, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = kmpp(Xt', k)
    @test iseeds_t == iseeds

    Random.seed!(34568)
    iseeds = kmpp_by_costs(C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = kmpp_by_costs(Ct', k)
    @test iseeds_t == iseeds
end

@testset "Kmcentrality" begin
    Random.seed!(34568)
    iseeds = initseeds(KmCentralityAlg(), X, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds(KmCentralityAlg(), Xt', k)
    @test iseeds == iseeds_t

    Random.seed!(34568)
    iseeds = initseeds_by_costs(KmCentralityAlg(), C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds_by_costs(KmCentralityAlg(), Ct', k)
    @test iseeds == iseeds_t

    @test min_interdist(X[:, iseeds]) > 2 * md0
    @test min_interdist((Xt')[:, iseeds]) > 2 * md0
end

end


text:
using Plots
using BenchmarkTools

function vfsolve(vnew, kgrid, tolerance, imax)
    α = 0.37
    δ = 0.1
    β = 0.95
    σ = 1.5

    v = vnew .+ 2*tolerance
    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1

    c = kgrid.^α + (1-δ)*kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid.^α + (1-δ)*kgrid .- kgrid' .< 0] .= -Inf

    while maximum(abs.(v - vnew)) > tolerance && i <= imax
        v = vnew;

        (vnew, cartesianindex) = findmax(u .+ β*v', dims = 2);
        i += 1;
    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex)
end

kupper = 5
klower = 0.05
n = 1000
kgrid = collect(range(klower, stop = kupper, length = n))
(v, kprime, kprimeindex) = vfsolve(zeros(n), kgrid, 0.001, 1000);

plot(kgrid, kprime, label = "k'")

function policy(x)
    e=getindex(findall(kgrid .== x),1)
    return kprime[e]
end

T = 100
kpath = zeros(T)
k = kgrid[92]
for i in 1:T
    k = policy(k)
    kpath[i] = k
end
time = collect(range(1, T, length = T))
plot(time, kpath, label = "kpath")

# Check k bounds, stepsize, tolerance, imax
findall(kprime .== maximum(kgrid))
any(kprime .== minimum(kgrid))

function vfsolve1(vnew, kgrid, imax)
    α = 1.5
    δ = 0.1
    β = 0.95
    σ = 1.5

    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1

    c = kgrid.^α + (1-δ)*kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid.^α + (1-δ)*kgrid .- kgrid' .< 0] .= -Inf

    while i <= imax

        (vnew, cartesianindex) = findmax(u .+ β*vnew', dims = 2);
        i += 1;

    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex)
end

(v1, kprime1, kprimeindex1) = vfsolve1(zeros(n), kgrid, 1);
(v5, kprime5, kprimeindex5) = vfsolve1(zeros(n), kgrid, 5);
(v10, kprime10, kprimeindex10) = vfsolve1(zeros(n), kgrid, 10);
scatter(kgrid, v1, label = "v1")
scatter(kgrid, v5, label = "v5")
scatter(kgrid, v10, label = "v10")

function vfsolve2(vnew, kgrid, tolerance, imax)
    α = 1.5
    δ = 0.1
    β = 0.95
    σ = 0.5

    v = vnew .+ 2*tolerance
    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1

    c = kgrid.^α + (1-δ)*kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid.^α + (1-δ)*kgrid .- kgrid' .< 0] .= -Inf

    while maximum(abs.(v - vnew)) > tolerance && i <= imax
        v = vnew;

        (vnew, cartesianindex) = findmax(u .+ β*v', dims = 2);
        i += 1;
    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex)
end

(v, kprime, kprimeindex) = vfsolve2(zeros(n), kgrid, 0.001, 1000);
scatter(kgrid, v, label = "v")
scatter(kgrid, kprime, label = "k'")

max_stars_repo_path: Dynamic Programming/DM.jl

max_issues_repo_path: Dynamic Programming/DM.jl

max_forks_repo_path: Dynamic Programming/DM.jl


text:
module Topsis

import GreyDecision.GreyNumbers.GreyNumber
import GreyDecision.Utility: normalize, unitize, weightize, apply_columns, reverseminmax, euclidean

function topsis(decisionMat::Array{GreyNumber{T}, 2}, weights::Array{Float64,1}, fns::Array{Function,1}) where {T <: Real}
    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    normalizedMat = normalize(decisionMat)

    weightednormalizedMat = weightize(normalizedMat, w)

    col_max = apply_columns(fns, weightednormalizedMat)
    col_min = apply_columns(reverseminmax(fns), weightednormalizedMat)

    distances_plus  = zeros(GreyNumber{Float64}, nalternatives)
    distances_minus = zeros(GreyNumber{Float64}, nalternatives)

    scores = zeros(GreyNumber{Float64}, nalternatives)

    for i in 1:nalternatives
        ithrow = weightednormalizedMat[i,:]
		distances_plus[i]  = euclidean(col_max, ithrow)
		distances_minus[i] = euclidean(col_min, ithrow)
		scores[i] = distances_minus[i] / (distances_minus[i] + distances_plus[i])
    end

    #best_index = sortperm(q) |> last
    return scores
end

end # End of module Topsis


text:
push!(LOAD_PATH,"../src/")

# using MiCRM
using Documenter
using MiCRM

DocMeta.setdocmeta!(MiCRM, :DocTestSetup, :(using MiCRM); recursive=true)

makedocs(;
    modules=[MiCRM],
    authors="Tom <t.clegg17@imperial.ac.uk> and contributors",
    repo="https://github.com/cleggtom/MiCRM.jl/blob/{commit}{path}#{line}",
    sitename="MiCRM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://cleggtom.github.io/MiCRM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => ["Basic Usage" => "pages/guide.md",
                    "Generating Communities" => "pages/community_generation.md",
                    "Building an ODESystem" => "pages/building_ODESystems.md",
                    "Coalescence" => "pages/coalescence.md"],
        "Library" => ["Public" => "pages/library/public.md",
                      "Internal" => "pages/library/internal.md"]
    ],
)

deploydocs(;
    repo="github.com/CleggTom/MiCRM.jl.git"
)


text:
function add_time_series!(
    component::T,
    time_series::TimeSeriesMetadata;
    skip_if_present = false,
) where {T <: InfrastructureSystemsComponent}
    component_name = get_name(component)
    container = get_time_series_container(component)
    if isnothing(container)
        throw(ArgumentError("type $T does not support storing time series"))
    end

    add_time_series!(container, time_series, skip_if_present = skip_if_present)
    @debug "Added $time_series to $(typeof(component)) $(component_name) " *
           "num_time_series=$(length(get_time_series_container(component).data))." _group =
        LOG_GROUP_TIME_SERIES
end

"""
Removes the metadata for a time_series.
If this returns true then the caller must also remove the actual time series data.
"""
function remove_time_series_metadata!(
    component::InfrastructureSystemsComponent,
    ::Type{T},
    name::AbstractString,
) where {T <: TimeSeriesMetadata}
    container = get_time_series_container(component)
    remove_time_series!(container, T, name)
    @debug "Removed time_series from $(get_name(component)):  $name." _group =
        LOG_GROUP_TIME_SERIES
    if T <: DeterministicMetadata &&
       has_time_series_internal(container, SingleTimeSeriesMetadata, name)
        return false
    elseif T <: SingleTimeSeriesMetadata &&
           has_time_series_internal(container, DeterministicMetadata, name)
        return false
    end

    return true
end

function clear_time_series!(component::InfrastructureSystemsComponent)
    container = get_time_series_container(component)
    if !isnothing(container)
        clear_time_series!(container)
        @debug "Cleared time_series in $(get_name(component))." _group =
            LOG_GROUP_TIME_SERIES
    end
end

function _get_columns(start_time, count, ts_metadata::ForecastMetadata)
    offset = start_time - get_initial_timestamp(ts_metadata)
    interval = time_period_conversion(get_interval(ts_metadata))
    window_count = get_count(ts_metadata)
    if window_count > 1
        index = Int(offset / interval) + 1
    else
        index = 1
    end
    if count === nothing
        count = window_count - index + 1
    end

    if index + count - 1 > get_count(ts_metadata)
        throw(
            ArgumentError(
                "The requested start_time $start_time and count $count are invalid",
            ),
        )
    end
    return UnitRange(index, index + count - 1)
end

_get_columns(start_time, count, ts_metadata::StaticTimeSeriesMetadata) = UnitRange(1, 1)

function _get_rows(start_time, len, ts_metadata::StaticTimeSeriesMetadata)
    index =
        Int(
            (start_time - get_initial_timestamp(ts_metadata)) / get_resolution(ts_metadata),
        ) + 1
    if len === nothing
        len = length(ts_metadata) - index + 1
    end
    if index + len - 1 > length(ts_metadata)
        throw(
            ArgumentError(
                "The requested index=$index len=$len exceeds the range $(length(ts_metadata))",
            ),
        )
    end

    return UnitRange(index, index + len - 1)
end

function _get_rows(start_time, len, ts_metadata::ForecastMetadata)
    if len === nothing
        len = get_horizon(ts_metadata)
    end

    return UnitRange(1, len)
end

function _check_start_time(start_time, ts_metadata::TimeSeriesMetadata)
    if start_time === nothing
        return get_initial_timestamp(ts_metadata)
    end

    time_diff = start_time - get_initial_timestamp(ts_metadata)
    if time_diff < Dates.Second(0)
        throw(
            ArgumentError(
                "start_time=$start_time is earlier than $(get_initial_timestamp(ts_metadata))",
            ),
        )
    end

    if typeof(ts_metadata) <: ForecastMetadata
        window_count = get_count(ts_metadata)
        interval = get_interval(ts_metadata)
        if window_count > 1 &&
           Dates.Millisecond(time_diff) % Dates.Millisecond(interval) != Dates.Second(0)
            throw(
                ArgumentError(
                    "start_time=$start_time is not on a multiple of interval=$interval",
                ),
            )
        end
    end

    return start_time
end

"""
Return a time series corresponding to the given parameters.

# Arguments
- `::Type{T}`: Concrete subtype of TimeSeriesData to return
- `component::InfrastructureSystemsComponent`: Component containing the time series
- `name::AbstractString`: name of time series
- `start_time::Union{Nothing, Dates.DateTime} = nothing`: If nothing, use the
  `initial_timestamp` of the time series. If T is a subtype of Forecast then `start_time`
  must be the first timstamp of a window.
- `len::Union{Nothing, Int} = nothing`: Length in the time dimension. If nothing, use the
  entire length.
- `count::Union{Nothing, Int} = nothing`: Only applicable to subtypes of Forecast. Number
  of forecast windows starting at `start_time` to return. Defaults to all available.
"""
function get_time_series(
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::AbstractString;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
    count::Union{Nothing, Int} = nothing,
) where {T <: TimeSeriesData}
    if !has_time_series(component)
        throw(ArgumentError("no forecasts are stored in $component"))
    end

    metadata_type = time_series_data_to_metadata(T)
    ts_metadata = get_time_series(metadata_type, component, name)
    start_time = _check_start_time(start_time, ts_metadata)
    rows = _get_rows(start_time, len, ts_metadata)
    columns = _get_columns(start_time, count, ts_metadata)
    storage = _get_time_series_storage(component)
    return deserialize_time_series(T, storage, ts_metadata, rows, columns)
end

function get_time_series(
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::AbstractString,
) where {T <: TimeSeriesMetadata}
    return get_time_series(T, get_time_series_container(component), name)
end

"""
Return a TimeSeries.TimeArray from storage for the given time series parameters.

If the data are scaling factors then the stored scaling_factor_multiplier will be called on
the component and applied to the data unless ignore_scaling_factors is true.
"""
function get_time_series_array(
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::AbstractString;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
    ignore_scaling_factors = false,
) where {T <: TimeSeriesData}
    ts = get_time_series(T, component, name; start_time = start_time, len = len, count = 1)
    if start_time === nothing
        start_time = get_initial_timestamp(ts)
    end

    return get_time_series_array(
        component,
        ts,
        start_time;
        len = len,
        ignore_scaling_factors = ignore_scaling_factors,
    )
end

"""
Return a TimeSeries.TimeArray for one forecast window from a cached Forecast instance.

If the data are scaling factors then the stored scaling_factor_multiplier will be called on
the component and applied to the data unless ignore_scaling_factors is true.

See also [`ForecastCache`](@ref).
"""
function get_time_series_array(
    component::InfrastructureSystemsComponent,
    forecast::Forecast,
    start_time::Dates.DateTime;
    len = nothing,
    ignore_scaling_factors = false,
)
    return _make_time_array(component, forecast, start_time, len, ignore_scaling_factors)
end

"""
Return a TimeSeries.TimeArray from a cached StaticTimeSeries instance.

If the data are scaling factors then the stored scaling_factor_multiplier will be called on
the component and applied to the data unless ignore_scaling_factors is true.

See also [`StaticTimeSeriesCache`](@ref).
"""
function get_time_series_array(
    component::InfrastructureSystemsComponent,
    time_series::StaticTimeSeries,
    start_time::Union{Nothing, Dates.DateTime} = nothing;
    len::Union{Nothing, Int} = nothing,
    ignore_scaling_factors = false,
)
    if start_time === nothing
        start_time = get_initial_timestamp(time_series)
    end

    if len === nothing
        len = length(time_series)
    end

    return _make_time_array(component, time_series, start_time, len, ignore_scaling_factors)
end

"""
Return a vector of timestamps from storage for the given time series parameters.
"""
function get_time_series_timestamps(
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::AbstractString;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
) where {T <: TimeSeriesData}
    return TimeSeries.timestamp(
        get_time_series_array(T, component, name; start_time = start_time, len = len),
    )
end

"""
Return a vector of timestamps from a cached Forecast instance.
"""
function get_time_series_timestamps(
    component::InfrastructureSystemsComponent,
    forecast::Forecast,
    start_time::Union{Nothing, Dates.DateTime} = nothing;
    len::Union{Nothing, Int} = nothing,
)
    return TimeSeries.timestamp(
        get_time_series_array(component, forecast, start_time; len = len),
    )
end

"""
Return a vector of timestamps from a cached StaticTimeSeries instance.
"""
function get_time_series_timestamps(
    component::InfrastructureSystemsComponent,
    time_series::StaticTimeSeries,
    start_time::Union{Nothing, Dates.DateTime} = nothing;
    len::Union{Nothing, Int} = nothing,
)
    return TimeSeries.timestamp(
        get_time_series_array(component, time_series, start_time; len = len),
    )
end

"""
Return an Array of values from storage for the requested time series parameters.

If the data size is small and this will be called many times, consider using the version
that accepts a cached TimeSeriesData instance.
"""
function get_time_series_values(
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::AbstractString;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
    ignore_scaling_factors = false,
) where {T <: TimeSeriesData}
    return TimeSeries.values(
        get_time_series_array(
            T,
            component,
            name;
            start_time = start_time,
            len = len,
            ignore_scaling_factors = ignore_scaling_factors,
        ),
    )
end

"""
Return an Array of values for one forecast window from a cached Forecast instance.
"""
function get_time_series_values(
    component::InfrastructureSystemsComponent,
    forecast::Forecast,
    start_time::Dates.DateTime;
    len::Union{Nothing, Int} = nothing,
    ignore_scaling_factors = false,
)
    return TimeSeries.values(
        get_time_series_array(
            component,
            forecast,
            start_time;
            len = len,
            ignore_scaling_factors = ignore_scaling_factors,
        ),
    )
end

"""
Return an Array of values from a cached StaticTimeSeries instance for the requested time
series parameters.
"""
function get_time_series_values(
    component::InfrastructureSystemsComponent,
    time_series::StaticTimeSeries,
    start_time::Union{Nothing, Dates.DateTime} = nothing;
    len::Union{Nothing, Int} = nothing,
    ignore_scaling_factors = false,
)
    return TimeSeries.values(
        get_time_series_array(
            component,
            time_series,
            start_time;
            len = len,
            ignore_scaling_factors = ignore_scaling_factors,
        ),
    )
end

function _make_time_array(component, time_series, start_time, len, ignore_scaling_factors)
    ta = make_time_array(time_series, start_time; len = len)
    if ignore_scaling_factors
        return ta
    end

    multiplier = get_scaling_factor_multiplier(time_series)
    if multiplier === nothing
        return ta
    end

    return ta .* multiplier(component)
end

"""
Return true if the component has time series data.
"""
function has_time_series(component::InfrastructureSystemsComponent)
    container = get_time_series_container(component)
    return !isnothing(container) && !isempty(container)
end

function has_time_series(
    component::InfrastructureSystemsComponent,
    type::Type{<:TimeSeriesMetadata},
    name::AbstractString,
)
    container = get_time_series_container(component)
    container === nothing && return false
    return has_time_series_internal(container, type, name)
end

"""
Efficiently add all time_series in one component to another by copying the underlying
references.

# Arguments
- `dst::InfrastructureSystemsComponent`: Destination component
- `src::InfrastructureSystemsComponent`: Source component
- `name_mapping::Dict = nothing`: Optionally map src names to different dst names.
  If provided and src has a time_series with a name not present in name_mapping, that
  time_series will not copied. If name_mapping is nothing then all time_series will be
  copied with src's names.
- `scaling_factor_multiplier_mapping::Dict = nothing`: Optionally map src multipliers to
  different dst multipliers.  If provided and src has a time_series with a multiplier not
  present in scaling_factor_multiplier_mapping, that time_series will not copied. If
  scaling_factor_multiplier_mapping is nothing then all time_series will be copied with
  src's multipliers.
"""
function copy_time_series!(
    dst::InfrastructureSystemsComponent,
    src::InfrastructureSystemsComponent;
    name_mapping::Union{Nothing, Dict{String, String}} = nothing,
    scaling_factor_multiplier_mapping::Union{Nothing, Dict{String, String}} = nothing,
)
    storage = _get_time_series_storage(dst)
    if isnothing(storage)
        throw(
            ArgumentError(
                "Component does not have time series storage. " *
                "It may not be attached to the system.",
            ),
        )
    end

    # There may be time series that share time series arrays as a result of
    # transform_single_time_series! being called.
    # Don't add these references to the storage more than once.
    refs = Set{Tuple{String, Base.UUID}}()

    for ts_metadata in get_time_series_multiple(TimeSeriesMetadata, src)
        name = get_name(ts_metadata)
        new_name = name
        if !isnothing(name_mapping)
            new_name = get(name_mapping, name, nothing)
            if isnothing(new_name)
                @debug "Skip copying ts_metadata" _group = LOG_GROUP_TIME_SERIES name
                continue
            end
            @debug "Copy ts_metadata with" _group = LOG_GROUP_TIME_SERIES new_name
        end
        multiplier = get_scaling_factor_multiplier(ts_metadata)
        new_multiplier = multiplier
        if !isnothing(scaling_factor_multiplier_mapping)
            new_multiplier = get(scaling_factor_multiplier_mapping, multiplier, nothing)
            if isnothing(new_multiplier)
                @debug "Skip copying ts_metadata" _group = LOG_GROUP_TIME_SERIES multiplier
                continue
            end
            @debug "Copy ts_metadata with" _group = LOG_GROUP_TIME_SERIES new_multiplier
        end
        new_time_series = deepcopy(ts_metadata)
        assign_new_uuid!(new_time_series)
        set_name!(new_time_series, new_name)
        set_scaling_factor_multiplier!(new_time_series, new_multiplier)
        add_time_series!(dst, new_time_series)
        ts_uuid = get_time_series_uuid(new_time_series)
        ref = (new_name, ts_uuid)
        if !in(ref, refs)
            add_time_series_reference!(storage, get_uuid(dst), new_name, ts_uuid)
            push!(refs, ref)
        end
    end
end

function get_time_series_keys(component::InfrastructureSystemsComponent)
    return keys(get_time_series_container(component).data)
end

function get_time_series_names(
    ::Type{T},
    component::InfrastructureSystemsComponent,
) where {T <: TimeSeriesData}
    return get_time_series_names(
        time_series_data_to_metadata(T),
        get_time_series_container(component),
    )
end

function get_num_time_series(component::InfrastructureSystemsComponent)
    container = get_time_series_container(component)
    if isnothing(container)
        return (0, 0)
    end

    static_ts_count = 0
    forecast_count = 0
    for key in keys(container.data)
        if key.time_series_type <: StaticTimeSeriesMetadata
            static_ts_count += 1
        elseif key.time_series_type <: ForecastMetadata
            forecast_count += 1
        else
            error("panic")
        end
    end

    return (static_ts_count, forecast_count)
end

function get_time_series(
    component::InfrastructureSystemsComponent,
    time_series::TimeSeriesData,
)
    storage = _get_time_series_storage(component)
    return get_time_series(storage, get_time_series_uuid(time_series))
end

function get_time_series_uuids(component::InfrastructureSystemsComponent)
    container = get_time_series_container(component)

    return [
        (get_time_series_uuid(container.data[key]), key.name) for
        key in get_time_series_keys(component)
    ]
end

"""
This function must be called when a component is removed from a system.
"""
function prepare_for_removal!(component::InfrastructureSystemsComponent)
    # TimeSeriesContainer can only be part of a component when that component is part of a
    # system.
    clear_time_series_storage!(component)
    set_time_series_storage!(component, nothing)
    clear_time_series!(component)
    @debug "cleared all time series data from" _group = LOG_GROUP_SYSTEM get_name(component)
end

"""
Returns an iterator of TimeSeriesData instances attached to the component.

Note that passing a filter function can be much slower than the other filtering parameters
because it reads time series data from media.

Call `collect` on the result to get an array.

# Arguments
- `component::InfrastructureSystemsComponent`: component from which to get time_series
- `filter_func = nothing`: Only return time_series for which this returns true.
- `type = nothing`: Only return time_series with this type.
- `name = nothing`: Only return time_series matching this value.
"""
function get_time_series_multiple(
    component::InfrastructureSystemsComponent,
    filter_func = nothing;
    type = nothing,
    start_time = nothing,
    name = nothing,
)
    container = get_time_series_container(component)
    storage = _get_time_series_storage(component)

    Channel() do channel
        for key in keys(container.data)
            ts_metadata = container.data[key]
            ts_type = time_series_metadata_to_data(ts_metadata)
            if !isnothing(type) && !(ts_type <: type)
                continue
            end
            if !isnothing(name) && key.name != name
                continue
            end
            ts = deserialize_time_series(
                ts_type,
                storage,
                ts_metadata,
                UnitRange(1, length(ts_metadata)),
                UnitRange(1, get_count(ts_metadata)),
            )
            if !isnothing(filter_func) && !filter_func(ts)
                continue
            end
            put!(channel, ts)
        end
    end
end

"""
Returns an iterator of TimeSeriesMetadata instances attached to the component.
"""
function get_time_series_multiple(
    ::Type{TimeSeriesMetadata},
    component::InfrastructureSystemsComponent,
)
    container = get_time_series_container(component)
    Channel() do channel
        for key in keys(container.data)
            put!(channel, container.data[key])
        end
    end
end

"""
Transform all instances of SingleTimeSeries to DeterministicSingleTimeSeries. Do nothing
if the component does not contain any instances.

All required checks must have been completed by the caller.

Return true if a transformation occurs.
"""
function transform_single_time_series_internal!(
    component::InfrastructureSystemsComponent,
    ::Type{T},
    params::TimeSeriesParameters,
) where {T <: DeterministicSingleTimeSeries}
    container = get_time_series_container(component)
    metadata_to_add = []
    for ts_metadata in values(container.data)
        if ts_metadata isa SingleTimeSeriesMetadata
            resolution = get_resolution(ts_metadata)
            new_metadata = DeterministicMetadata(
                name = get_name(ts_metadata),
                resolution = params.resolution,
                initial_timestamp = params.forecast_params.initial_timestamp,
                interval = params.forecast_params.interval,
                count = params.forecast_params.count,
                time_series_uuid = get_time_series_uuid(ts_metadata),
                horizon = params.forecast_params.horizon,
                time_series_type = DeterministicSingleTimeSeries,
                scaling_factor_multiplier = get_scaling_factor_multiplier(ts_metadata),
                internal = get_internal(ts_metadata),
            )
            push!(metadata_to_add, new_metadata)
        end
    end

    isempty(metadata_to_add) && return false

    for new_metadata in metadata_to_add
        add_time_series!(container, new_metadata)
        @debug "Added $new_metadata." _group = LOG_GROUP_TIME_SERIES
    end

    return true
end

function get_single_time_series_transformed_parameters(
    component::InfrastructureSystemsComponent,
    ::Type{T},
    horizon::Int,
    interval::Dates.Period,
) where {T <: Forecast}
    container = get_time_series_container(component)
    for (key, ts_metadata) in container.data
        if ts_metadata isa SingleTimeSeriesMetadata
            return _get_single_time_series_transformed_parameters(
                ts_metadata,
                T,
                horizon,
                interval,
            )
        end
    end

    throw(ArgumentError("component $(get_name(component)) does not have SingleTimeSeries"))
end

function _get_single_time_series_transformed_parameters(
    ts_metadata::SingleTimeSeriesMetadata,
    ::Type{T},
    horizon::Int,
    interval::Dates.Period,
) where {T <: Forecast}
    resolution = get_resolution(ts_metadata)
    len = length(ts_metadata)
    if len < horizon
        throw(
            ConflictingInputsError("existing length=$len is shorter than horizon=$horizon"),
        )
    end

    max_interval = horizon * resolution
    if len == horizon && interval == max_interval
        interval = Dates.Second(0)
        @warn "There is only one forecast window. Setting interval = $interval"
    elseif interval > max_interval
        throw(
            ConflictingInputsError(
                "interval = $interval is bigger than the max of $max_interval",
            ),
        )
    end

    initial_timestamp = get_initial_timestamp(ts_metadata)
    return TimeSeriesParameters(initial_timestamp, resolution, len, horizon, interval)
end

function clear_time_series_storage!(component::InfrastructureSystemsComponent)
    storage = _get_time_series_storage(component)
    if !isnothing(storage)
        # In the case of Deterministic and DeterministicSingleTimeSeries the UUIDs
        # can be shared.
        uuids = Set{Base.UUID}()
        for (uuid, name) in get_time_series_uuids(component)
            if !(uuid in uuids)
                remove_time_series!(storage, uuid, get_uuid(component), name)
                push!(uuids, uuid)
            end
        end
    end
end

function set_time_series_storage!(
    component::InfrastructureSystemsComponent,
    storage::Union{Nothing, TimeSeriesStorage},
)
    container = get_time_series_container(component)
    if !isnothing(container)
        set_time_series_storage!(container, storage)
    end
end

function _get_time_series_storage(component::InfrastructureSystemsComponent)
    container = get_time_series_container(component)
    if isnothing(container)
        return nothing
    end

    return container.time_series_storage
end

function get_time_series_by_key(
    key::TimeSeriesKey,
    component::InfrastructureSystemsComponent;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
    count::Union{Nothing, Int} = nothing,
)
    container = get_time_series_container(component)
    ts_metadata = container.data[key]
    ts_type = time_series_metadata_to_data(ts_metadata)
    return get_time_series(
        ts_type,
        component,
        key.name,
        start_time = start_time,
        len = len,
        count = count,
    )
end


text:
using NetCDF, LinearAlgebra
using Statistics, StatsBase
using PyCall
# xr = pyimport("xarray")
#pygui(:qt)
using PyPlot
using MultivariateStats
# eofs = pyimport("eofs")

path = "/Users/milan/Downloads/cams/ieee"
filelist = filter(x->endswith(x,"no2.nc"),readdir(path))

ncfile = NetCDF.open(joinpath(path,filelist[1]))
no2 = ncfile.vars["no2"][:,:,:,1]
level = ncfile.vars["level"][:]
lat = ncfile.vars["latitude"][:]
lon = ncfile.vars["longitude"][:]

no2m = mean(no2,dims=(1,2))[1,1,:]
no2mi = median(no2,dims=(1,2))[1,1,:]
no2min = minimum(no2,dims=(1,2))[1,1,:]
no2max = maximum(no2,dims=(1,2))[1,1,:]
no2p10 = [percentile(vec(no2[:,:,i]),10) for i in level]
no2p90 = [percentile(vec(no2[:,:,i]),90) for i in level]

## EOF in the vertical
no2f = Matrix(reshape(no2,(:,length(level)))')
M = fit(PCA,log10.(no2f),maxoutdim=20,pratio=0.999)
no2fhat = transform(M,log10.(no2f))
no2fr = 10f0 .^reconstruct(M,no2fhat)

## reconstruction statistics
no2rm = mean(no2fr,dims=2)[:,1]
no2rmi = median(no2fr,dims=2)[:,1]
no2rmin = minimum(no2fr,dims=2)[:,1]
no2rmax = maximum(no2fr,dims=2)[:,1]
no2rp10 = [percentile(no2fr[i,:],10) for i in level]
no2rp90 = [percentile(no2fr[i,:],90) for i in level]

## error
decimal_error = abs.(log10.(reshape(no2,(:,length(level)))' ./ no2fr))

dem = mean(decimal_error,dims=2)[:,1]
demi = median(decimal_error,dims=2)[:,1]
demin = minimum(decimal_error,dims=2)[:,1]
demax = maximum(decimal_error,dims=2)[:,1]
dep10 = [percentile(decimal_error[i,:],10) for i in level]
dep90 = [percentile(decimal_error[i,:],90) for i in level]

## highest loadings
P = projection(M)
Ps = zeros(Int64,size(P))
for i in 1:size(Ps)[2]
    Ps[:,i] = sortperm(abs.(P[:,i]))
end
L = zeros(Int64,size(P)[2])
L[1] = Ps[end,1]
for i in 2:length(L)
    j = 0
    while Ps[end-j,i] in L
        j += 1
    end
    L[i] = Ps[end-j,i]
end

##
fig,(ax,ax2,ax3) = subplots(1,3,figsize=(10,6),sharey=true)
ax.invert_yaxis()

ax.semilogx(no2m,level,"C0",lw=3,label="mean")
ax.plot(no2mi,level,"C0",lw=1.5,label="median")
ax.fill_betweenx(level,no2min,no2max,color="C0",alpha=0.2,label="Min-max range")
ax.fill_betweenx(level,no2p10,no2p90,color="C0",alpha=0.5,label="10-90% range")

ax2.semilogx(no2rm,level,"C0",lw=3,label="mean")
ax2.plot(no2rmi,level,"C0",lw=1.5,label="median")
ax2.fill_betweenx(level,no2rmin,no2rmax,color="C0",alpha=0.2,label="Min-max range")
ax2.fill_betweenx(level,no2rp10,no2rp90,color="C0",alpha=0.5,label="10-90% range")

for i in level
    xvals = no2f[i,(decimal_error[i,:] .>= 0.5f0) .* (decimal_error[i,:] .< 1f0)]
    ax.scatter(xvals,i*ones(length(xvals)),10,"yellow",alpha=0.05)

    xvals = no2f[i,(decimal_error[i,:] .>= 1f0) .* (decimal_error[i,:] .< 2f0)]
    ax.scatter(xvals,i*ones(length(xvals)),10,"C1",alpha=0.1)

    xvals = no2f[i,(decimal_error[i,:] .>= 2f0)]
    ax.scatter(xvals,i*ones(length(xvals)),10,"C3",alpha=0.1)
end

ax.scatter([1,2,3],[0,0,0],10,"yellow",alpha=0.5,label="error > 0.5")
ax.scatter([1,2,3],[0,0,0],10,"C1",alpha=0.5,label="error > 1")
ax.scatter([1,2,3],[0,0,0],10,"C3",alpha=0.5,label="error > 2")

ax3.semilogx(dem,level,"C0",lw=3,label="mean")
ax3.plot(demi,level,"C0",lw=1.5,label="median")
ax3.fill_betweenx(level,demin,demax,color="C0",alpha=0.2,label="Min-max range")
ax3.fill_betweenx(level,dep10,dep90,color="C0",alpha=0.5,label="10-90% range")

cmap = ColorMap("viridis")
for i in 1:length(L)
    ax2.plot([1e-15,1e-14],[L[i],L[i]].-0.5,color=cmap(i/length(L)),lw=2)
    ax2.text(2e-14,L[i],"$i")
end
ax2.text(1.5e-15,8,"modes")

ax.set_title(L"NO$_2$ vertical distribution")
ax2.set_title(L"NO$_2$ reconstructed from 20 PCs")
ax3.set_title("L1 log-error")
ax.set_ylabel("model level")
ax.set_xlabel("mixing ratio kg/kg")
ax2.set_xlabel("mixing ratio kg/kg")
ax.set_ylim(137,1)
ax.legend(loc=2,scatterpoints=3)

ax.set_xlim(1e-15,1e-7)
ax2.set_xlim(1e-15,1e-7)

tight_layout()


text:
using Base.Test
using Dagger
using MemPool
using Base.Test

@testset "cache" begin
    @everywhere gc(true)
    # set available memory to 80KB on each worker
    @everywhere empty!(MemPool.datastore)
    @everywhere empty!(MemPool.lru_order)
    @everywhere MemPool.max_memsize[] = 8*10^4

    thunks1 = map(delayed(_ -> rand(10^3), cache=true), workers())
    sum1 = delayed((x...)->sum([x...]))(map(delayed(sum), thunks1)...)
    thunks2 = map(delayed(-), thunks1)
    sum2 = delayed((x...)->sum([x...]))(map(delayed(sum), thunks2)...)
    s1 = collect(sum1)
    @test -s1 == collect(sum2)
    @test s1 == collect(sum1)
    @test -collect(sum1) == collect(sum2)

    thunks1 = map(delayed(_ -> rand(10^4), cache=true), workers())
    sum1 = delayed((x...)->sum([x...]))(map(delayed(sum), thunks1)...)
    s1 = collect(sum1)
    thunks2 = map(delayed(_ -> rand(10^4), cache=true), workers())
    sum2 = delayed((x...)->sum([x...]))(map(delayed(sum), thunks2)...)
    s2 = collect(sum2) # this should evict thunk1s from memory
    @test s1 != collect(sum1)
    @test s2 != collect(sum2)
end


text:
## Uniform ##

logpdf(d::Uniform, x::TrackedReal) = uniformlogpdf(d.a, d.b, x)
uniformlogpdf(a, b, x) = Tracker.track(uniformlogpdf, a, b, x)
Tracker.@grad function uniformlogpdf(a, b, x)
    xd = Tracker.data(x)
    T = typeof(xd)
    l = logpdf(Uniform(a, b), Tracker.data(x))
    f = isfinite(l)
    da = 1/(b - a)
    n = T(NaN)
    return l, Δ->(f ? da : n, f ? -da : n, f ? zero(T) : n)
end

## Semicircle ##

logpdf(d::Semicircle{<:Real}, x::TrackedReal) = semicirclelogpdf(d.r, x)
logpdf(d::Semicircle{<:TrackedReal}, x::Real) = semicirclelogpdf(d.r, x)
logpdf(d::Semicircle{<:TrackedReal}, x::TrackedReal) = semicirclelogpdf(d.r, x)
semicirclelogpdf(r::TrackedReal, x::Real) = Tracker.track(semicirclelogpdf, r, x)
semicirclelogpdf(r::Real, x::TrackedReal) = Tracker.track(semicirclelogpdf, r, x)
semicirclelogpdf(r::TrackedReal, x::TrackedReal) = Tracker.track(semicirclelogpdf, r, x)
Tracker.@grad function semicirclelogpdf(r, x)
    rd = Tracker.data(r)
    xd = Tracker.data(x)
    xx, rr = promote(xd, float(rd))
    d = Semicircle(rr)
    T = typeof(xx)
    l = logpdf(d, xx)
    f = isfinite(l)
    n = T(NaN)
    return l, function (Δ)
        diffsq = rr^2 - xx^2
        (f ? Δ*(-2/rr + rr/diffsq) : n, f ? Δ*(-xx/diffsq) : n)
    end
end

## Binomial ##

binomlogpdf(n::Int, p::Tracker.TrackedReal, x::Int) = Tracker.track(binomlogpdf, n, p, x)
Tracker.@grad function binomlogpdf(n::Int, p::Tracker.TrackedReal, x::Int)
    return binomlogpdf(n, Tracker.data(p), x),
        Δ->(nothing, Δ * (x / p - (n - x) / (1 - p)), nothing)
end

function binomlogpdf(n::Int, p::ForwardDiff.Dual{T}, x::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val = ForwardDiff.value(p)
    Δ = ForwardDiff.partials(p)
    return FD(binomlogpdf(n, val, x),  Δ * (x / val - (n - x) / (1 - val)))
end

## Negative binomial ##

# Note the definition of NegativeBinomial in Julia is not the same as Wikipedia's.
# Check the docstring of NegativeBinomial, r is the number of successes and
# k is the number of failures
_nbinomlogpdf_grad_1(r, p, k) = k == 0 ? log(p) : sum(1 / (k + r - i) for i in 1:k) + log(p)
_nbinomlogpdf_grad_2(r, p, k) = -k / (1 - p) + r / p

nbinomlogpdf(n::Tracker.TrackedReal, p::Tracker.TrackedReal, x::Int) = Tracker.track(nbinomlogpdf, n, p, x)
nbinomlogpdf(n::Real, p::Tracker.TrackedReal, x::Int) = Tracker.track(nbinomlogpdf, n, p, x)
nbinomlogpdf(n::Tracker.TrackedReal, p::Real, x::Int) = Tracker.track(nbinomlogpdf, n, p, x)
Tracker.@grad function nbinomlogpdf(r::Tracker.TrackedReal, p::Tracker.TrackedReal, k::Int)
    return nbinomlogpdf(Tracker.data(r), Tracker.data(p), k),
        Δ->(Δ * _nbinomlogpdf_grad_1(r, p, k), Δ * _nbinomlogpdf_grad_2(r, p, k), nothing)
end
Tracker.@grad function nbinomlogpdf(r::Real, p::Tracker.TrackedReal, k::Int)
    return nbinomlogpdf(Tracker.data(r), Tracker.data(p), k),
        Δ->(Tracker._zero(r), Δ * _nbinomlogpdf_grad_2(r, p, k), nothing)
end
Tracker.@grad function nbinomlogpdf(r::Tracker.TrackedReal, p::Real, k::Int)
    return nbinomlogpdf(Tracker.data(r), Tracker.data(p), k),
        Δ->(Δ * _nbinomlogpdf_grad_1(r, p, k), Tracker._zero(p), nothing)
end

function nbinomlogpdf(r::ForwardDiff.Dual{T}, p::ForwardDiff.Dual{T}, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_p = ForwardDiff.value(p)
    val_r = ForwardDiff.value(r)

    Δ_r = ForwardDiff.partials(r) * _nbinomlogpdf_grad_1(val_r, val_p, k)
    Δ_p = ForwardDiff.partials(p) * _nbinomlogpdf_grad_2(val_r, val_p, k)
    Δ = Δ_p + Δ_r
    return FD(nbinomlogpdf(val_r, val_p, k),  Δ)
end
function nbinomlogpdf(r::Real, p::ForwardDiff.Dual{T}, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_p = ForwardDiff.value(p)
    Δ_p = ForwardDiff.partials(p) * _nbinomlogpdf_grad_2(r, val_p, k)
    return FD(nbinomlogpdf(r, val_p, k),  Δ_p)
end
function nbinomlogpdf(r::ForwardDiff.Dual{T}, p::Real, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_r = ForwardDiff.value(r)
    Δ_r = ForwardDiff.partials(r) * _nbinomlogpdf_grad_1(val_r, p, k)
    return FD(nbinomlogpdf(val_r, p, k),  Δ_r)
end

## Poisson ##

poislogpdf(v::Tracker.TrackedReal, x::Int) = Tracker.track(poislogpdf, v, x)
Tracker.@grad function poislogpdf(v::Tracker.TrackedReal, x::Int)
      return poislogpdf(Tracker.data(v), x),
          Δ->(Δ * (x/v - 1), nothing)
end

function poislogpdf(v::ForwardDiff.Dual{T}, x::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val = ForwardDiff.value(v)
    Δ = ForwardDiff.partials(v)
    return FD(poislogpdf(val, x), Δ * (x/val - 1))
end

## PoissonBinomial ##

struct TuringPoissonBinomial{T<:Real, TV<:AbstractVector{T}} <: DiscreteUnivariateDistribution
    p::TV
    pmf::TV
end
function TuringPoissonBinomial(p::AbstractArray{<:Real})
    pb = Distributions.poissonbinomial_pdf_fft(p)
    @assert Distributions.isprobvec(pb)
    TuringPoissonBinomial(p, pb)
end
function logpdf(d::TuringPoissonBinomial{T}, k::Int) where T<:Real
    insupport(d, k) ? log(d.pmf[k + 1]) : -T(Inf)
end
quantile(d::TuringPoissonBinomial, x::Float64) = quantile(Categorical(d.pmf), x) - 1
PoissonBinomial(p::Tracker.TrackedArray) = TuringPoissonBinomial(p)
Base.minimum(d::TuringPoissonBinomial) = 0
Base.maximum(d::TuringPoissonBinomial) = length(d.p)

poissonbinomial_pdf_fft(x::Tracker.TrackedArray) = Tracker.track(poissonbinomial_pdf_fft, x)
Tracker.@grad function poissonbinomial_pdf_fft(x::Tracker.TrackedArray)
    x_data = Tracker.data(x)
    T = eltype(x_data)
    fft = poissonbinomial_pdf_fft(x_data)
    return  fft, Δ -> begin
        ((ForwardDiff.jacobian(x -> poissonbinomial_pdf_fft(x), x_data)::Matrix{T})' * Δ,)
    end
end


text:
module TFIMSampler

using BitBasis
using DelimitedFiles
using ProgressMeter
using Random
using StatsBase
using LinearAlgebra
using SparseArrays

set_zero_subnormals(true)

include("new.jl")
include("mcmc.jl")

export sample, DQMC, pairing, sweep!, fast_update!, get_probs, get_prob, get_all_probs

end # module


text:
s = BenchmarkGroup()
suite["generators"] = s

srand(17)
k=3; seed=17
for n in [100,500], G in GLIST
    g = random_regular_graph(n, k, G, seed=seed)
    s["rrg","$g"] = @benchmarkable random_regular_graph($n, $k, $G, seed=$seed)

    g = erdos_renyi(n, k*n, G, seed=seed)
    s["erdos","$g"] = @benchmarkable erdos_renyi($n, $(k*n), $G, seed=$seed)
end

for n in [20,100], G in GLIST
    g = CompleteGraph(n, G)
    s["complete","$g"] = @benchmarkable CompleteGraph($n, $G)

    g = CompleteDiGraph(n, digraphtype(G))
    s["complete","$g"] = @benchmarkable CompleteDiGraph($n, $(digraphtype(G)))
end


text:
using LinearAlgebra

mutable struct svd_kalman
    A # state transition
    G #Observable transition
    Q  #observable variance
    R#variance of the state
    state_mean
    state_var
end

function start!(svd_kalman::svd_kalman,x0,var0)

    R = svd_kalman.R
    Q = svd_kalman.Q

    svd_kalman.state_mean = x0
    svd_var0 = svd(var0)
    svd_kalman.state_var = diagm(sqrt.(svd_var0.S))*svd_var0.Vt

    svd_R = svd(R)
    svd_kalman.R = diagm(sqrt.(svd_R.S))*svd_R.Vt

    svd_Q = svd(Q)
    svd_kalman.Q = diagm(sqrt.(svd_Q.S))*svd_Q.Vt

    return nothing
end

function update_state!(svd_kalman::svd_kalman,y;ll=false)

    R = svd_kalman.R
    old_var = svd_kalman.state_var
    G = svd_kalman.G
    old_state = svd_kalman.state_mean

    pre_gain = [R;old_var*G']
    svd_re = svd(pre_gain)

    K_bar = (old_var'*old_var)*G'*svd_re.V

    err = (y - G*old_state)
    err_bar = svd_re.Vt*err

    if size(old_state,1) == 1
        new_state = old_state + (K_bar*(svd_re.S).^(-2)*err_bar)[1]
    else
        new_state = old_state + K_bar*diagm((svd_re.S).^(-2))*err_bar
    end

    K = K_bar*diagm((svd_re.S).^(-2))*svd_re.Vt

    pre_var = [old_var*(I-K*G)';R*K']
    svd_var = svd(pre_var)

    new_var = diagm(svd_var.S)*svd_var.Vt

    svd_kalman.state_mean = new_state
    svd_kalman.state_var = new_var

    if ll == false
        return nothing
    else
        Re = diagm(svd_re.S)*svd_re.Vt
        return err, Re'*Re
    end
end

function forecast!(svd_kalman::svd_kalman)

    A = svd_kalman.A
    state_mean = svd_kalman.state_mean

    svd_kalman.state_mean = A*state_mean

    vari = svd_kalman.state_var
    Q = svd_kalman.Q

    pre_var = [vari*A';Q]
    svd_var = svd(pre_var)

    svd_kalman.state_var = diagm(svd_var.S)*svd_var.Vt

    return nothing

end

function update_and_forecast!(svd_kalman::svd_kalman,y;ll=false)
    if ll == false
        update_state!(svd_kalman,y)
        forecast!(svd_kalman)
        return nothing
    else
        μ,Σ = update_state!(svd_kalman,y,ll=true)
        forecast!(svd_kalman)
        return μ, Σ
    end
end

#loglikelihood
function loglikelihood(svd_kalman::svd_kalman,data)
    n = size(data,1)
    m = size(data,2)

    llh = zeros(n)

    G = svd_kalman.G
    R = svd_kalman.R
    μ = G*svd_kalman.state_mean
    Σ = G*svd_kalman.state_var*G' + R'*R

    for j in 1:n
        if m == 1
            llh[j] = -1/2*(m*log(2*pi) + logdet(Σ) + μ[1]^2/Σ[1])
        else
            llh[j] = -1/2*(m*log(2*pi) + logdet(Σ) + μ'*inv(Σ)*μ)
        end

        μ,Σ = update_and_forecast!(svd_kalman,data[j,:],ll=true)
    end
    return llh
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using Test
using DataWrangler

@testset "p" begin

    x = [1,2,3,4,5]
    x0 = [0,1,2,3,4]
    @test p(d(x)) == x0
    @test p(d(x),[x[1]]) == x
    @test p(d(hcat(x,x,x))) == hcat(x0,x0,x0)
    @test p(d(hcat(x,x,x)), [ x[1] x[1] x[1] ] ) == hcat(x,x,x)

    x = rand(100);
    @test p(d(x),[ x[1] ]) ≈ x
    @test p(d(hcat(x,x,x)), [ x[1] x[1] x[1] ] ) ≈ hcat(x,x,x)

    x = [repeat(1:10,10) repeat(1:10,10)];
    dx = d(x,3);
    x0 = [[1,1,0] [1,1,0]]
    @test p(dx, x0) == x

    lag = 11
    x = [repeat(1:lag,14) repeat(1:lag,14)
         1 1];
    dx = d(x,1,lag);
    x0 = reshape(repeat(1:lag,inner=2),1,2,lag);
    @test p(dx, x0) == x

    x = repeat(1:2,30);
    dx = d(x,2,2);
    x0 = zeros(2,1,2); # order 2, 1 variable, lag 1
    x0[1,:,:] = collect(1:2);
    @test p(dx,x0) ≈ x

    # calculation of π
    x = 0:0.001:1
    y = sqrt.(1 .- x.^2);
    isapprox(4*p(y)[end]/1000 , π, atol = 0.01)

end


text:
module RingPolymers

include("ring_polymer_array.jl")
export RingPolymerArray,
       get_centroid

include("ring_polymer.jl")
export transform_to_normal_modes!,
       transform_from_normal_modes!,
       nbeads

end # module


text:
"""
    $(SIGNATURES)

An in-memory DistributedFactorGraph based on LightGraphs.jl with parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
mutable struct LightDFG{T <: AbstractParams, V <: AbstractDFGVariable, F <:AbstractDFGFactor} <: AbstractDFG{T}
    g::FactorGraph{Int, V, F}
    description::String
    userId::String
    robotId::String
    sessionId::String
    userData::Dict{Symbol, String}
    robotData::Dict{Symbol, String}
    sessionData::Dict{Symbol, String}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
end

"""
    $(SIGNATURES)

Create an in-memory LightDFG with the following parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
function LightDFG{T,V,F}(g::FactorGraph{Int,V,F}=FactorGraph{Int,V,F}();
                           description::String="LightGraphs.jl implementation",
                           userId::String="DefaultUser",
                           robotId::String="DefaultRobot",
                           sessionId::String="Session_$(string(uuid4())[1:6])",
                           userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           solverParams::T=T()) where {T <: AbstractParams, V <:AbstractDFGVariable, F<:AbstractDFGFactor}
   # Validate the userId, robotId, and sessionId
   !isValidLabel(userId) && error("'$userId' is not a valid User ID")
   !isValidLabel(robotId) && error("'$robotId' is not a valid Robot ID")
   !isValidLabel(sessionId) && error("'$sessionId' is not a valid Session ID")
   return LightDFG{T,V,F}(g, description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], solverParams)
end

# LightDFG{T}(; kwargs...) where T <: AbstractParams = LightDFG{T,DFGVariable,DFGFactor}(;kwargs...)

"""
    $(SIGNATURES)

Create an in-memory LightDFG with the following parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
function LightDFG{T}(g::FactorGraph{Int,DFGVariable,DFGFactor}=FactorGraph{Int,DFGVariable,DFGFactor}();
                     params=nothing,
                     kwargs...) where T <: AbstractParams

    #TODO remove params deprecation error in v0.9
    if params != nothing
        @warn "keyword `params` is deprecated, please use solverParams"
        return LightDFG{T,DFGVariable,DFGFactor}(g; solverParams = params, kwargs...)
    end
    return LightDFG{T,DFGVariable,DFGFactor}(g; kwargs...)
end

function LightDFG(g::FactorGraph{Int,DFGVariable,DFGFactor}=FactorGraph{Int,DFGVariable,DFGFactor}();
                  solverParams::T=NoSolverParams(), params=nothing, kwargs...) where T
    #TODO remove params deprecation error in v0.9
    if params != nothing
        @warn "keyword `params` is deprecated, please use solverParams"
        return LightDFG{typeof(params),DFGVariable,DFGFactor}(g; solverParams=params, kwargs...)
    end
    return LightDFG{T,DFGVariable,DFGFactor}(g; solverParams=solverParams, kwargs...)
end

LightDFG(description::String,
         userId::String,
         robotId::String,
         sessionId::String,
         userData::Dict{Symbol, String},
         robotData::Dict{Symbol, String},
         sessionData::Dict{Symbol, String},
         solverParams::AbstractParams) =
         LightDFG(FactorGraph{Int,DFGVariable,DFGFactor}(), description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], solverParams)

LightDFG{T,V,F}(description::String,
                userId::String,
                robotId::String,
                sessionId::String,
                userData::Dict{Symbol, String},
                robotData::Dict{Symbol, String},
                sessionData::Dict{Symbol, String},
                solverParams::T) where {T <: AbstractParams, V <:AbstractDFGVariable, F<:AbstractDFGFactor} =
                LightDFG(FactorGraph{Int,V,F}(), description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], solverParams)


text:
module SCS

if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    include(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
else
    error("SCS not properly installed. Please run Pkg.build(\"SCS\") and restart julia")
end

function __init__()
    vnum = VersionNumber(SCS_version())
    depsdir = realpath(joinpath(dirname(@__FILE__), "..", "deps"))
    if vnum.major != 2 || vnum.minor != 0
        error("Current SCS version installed is $(SCS_version()), but we require version 2.0.*.")
    end
end

include("types.jl")
include("c_wrapper.jl")
include("MPB_wrapper.jl")
include("MOI_wrapper.jl")

end # module


text:
@testset "Gadfly-based plot Test" begin
  # fist on a tree, then on a network with h=2 and various options"
  # keep ";" below to avoid creating a new file
  tre = readTopology("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
  net = readTopology("(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,(((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);");
  @test_logs plot(tre);
  @test_logs plot(net);
  @test_logs (:warn, "At least one non-missing edge length: plotting any missing length as 1.0") plot(net, useEdgeLength=true);
  @test_logs plot(net, mainTree=true);
  @test_logs plot(net, showTipLabel=false);
  @test_logs plot(net, showNodeNumber=true);
  @test_logs plot(net, showEdgeLength=false, showEdgeNumber=true);
  @test_logs plot(net, showGamma=true);
  @test_logs plot(net, edgeColor=colorant"olive", # uses Colors
          minorHybridEdgeColor=colorant"tan",
          majorHybridEdgeColor=colorant"skyblue");

  # example with 2-cycle (simple loop)
  cui3str2 = "(Xalvarezi,(((Xclemenciae_F2,Xmonticolus):1.458,(((((Xmontezumae,(Xnezahuacoyotl)#H26:0.247::0.804):0.375,((Xbirchmanni_GARC,Xmalinche_CHIC2):0.997,Xcortezi):0.455):0.63,(#H26:0.0::0.196,((Xcontinens,Xpygmaeus):1.932,(Xnigrensis,Xmultilineatus):1.401):0.042):2.439):2.0)#H7:0.787::0.835,(Xmaculatus,(Xandersi,(Xmilleri,((Xxiphidium,#H7:9.563::0.165):1.409,(Xevelynae,(Xvariatus,(Xcouchianus,(Xgordoni,Xmeyeri):0.263):3.532):0.642):0.411):0.295):0.468):0.654):1.022):0.788):1.917)#H27:0.817::0.572,#H27:6.307::0.428);"
  net3  = readTopology(cui3str2);
  @test_logs plot(net3); # k=2 cycle at the root. 3 root edges:
  # one to a leaf, 1 major & 1 minor hybrid edge to the same child.

end


text:
using Turing
import Turing.translate!

ex = quote
  x = 1
  y = rand()
  y ~ Normal(0,1)
end

res = translate!(:(y~Normal(1,1)))

Base.@assert res.head == :macrocall
Base.@assert res.args[1] == Symbol("@~")
Base.@assert res.args[3] == :y
Base.@assert res.args[4] == :(Normal(1, 1))

res2 = translate!(ex)

Base.@assert res2.args[end].head == :macrocall
Base.@assert res2.args[end].args[1] == Symbol("@~")
Base.@assert res2.args[end].args[3] == :y
Base.@assert res2.args[end].args[4] == :(Normal(0, 1))


text:
using SafeTestsets

@safetestset "Unconstrained tests" begin include("unconstrained_tests.jl") end
@safetestset "Electromechanical tests" begin include("electromechanical_tests.jl") end
@safetestset "Constrained tests" begin include("constrained_tests.jl") end