text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsUserProfileShort(;
        avatar_hash=nothing,
        display_name=nothing,
        display_name_normalized=nothing,
        image_72=nothing,
        is_restricted=nothing,
        is_ultra_restricted=nothing,
        name=nothing,
        real_name=nothing,
        real_name_normalized=nothing,
        team=nothing,
    )

    - avatar_hash::String
    - display_name::String
    - display_name_normalized::String
    - image_72::String
    - is_restricted::Bool
    - is_ultra_restricted::Bool
    - name::String
    - real_name::String
    - real_name_normalized::String
    - team::DefsWorkspaceId
"""
mutable struct ObjsUserProfileShort <: SwaggerModel
    avatar_hash::Any # spec type: Union{ Nothing, String } # spec name: avatar_hash
    display_name::Any # spec type: Union{ Nothing, String } # spec name: display_name
    display_name_normalized::Any # spec type: Union{ Nothing, String } # spec name: display_name_normalized
    image_72::Any # spec type: Union{ Nothing, String } # spec name: image_72
    is_restricted::Any # spec type: Union{ Nothing, Bool } # spec name: is_restricted
    is_ultra_restricted::Any # spec type: Union{ Nothing, Bool } # spec name: is_ultra_restricted
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    real_name::Any # spec type: Union{ Nothing, String } # spec name: real_name
    real_name_normalized::Any # spec type: Union{ Nothing, String } # spec name: real_name_normalized
    team::Any # spec type: Union{ Nothing, DefsWorkspaceId } # spec name: team

    function ObjsUserProfileShort(;avatar_hash=nothing, display_name=nothing, display_name_normalized=nothing, image_72=nothing, is_restricted=nothing, is_ultra_restricted=nothing, name=nothing, real_name=nothing, real_name_normalized=nothing, team=nothing)
        o = new()
        validate_property(ObjsUserProfileShort, Symbol("avatar_hash"), avatar_hash)
        setfield!(o, Symbol("avatar_hash"), avatar_hash)
        validate_property(ObjsUserProfileShort, Symbol("display_name"), display_name)
        setfield!(o, Symbol("display_name"), display_name)
        validate_property(ObjsUserProfileShort, Symbol("display_name_normalized"), display_name_normalized)
        setfield!(o, Symbol("display_name_normalized"), display_name_normalized)
        validate_property(ObjsUserProfileShort, Symbol("image_72"), image_72)
        setfield!(o, Symbol("image_72"), image_72)
        validate_property(ObjsUserProfileShort, Symbol("is_restricted"), is_restricted)
        setfield!(o, Symbol("is_restricted"), is_restricted)
        validate_property(ObjsUserProfileShort, Symbol("is_ultra_restricted"), is_ultra_restricted)
        setfield!(o, Symbol("is_ultra_restricted"), is_ultra_restricted)
        validate_property(ObjsUserProfileShort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ObjsUserProfileShort, Symbol("real_name"), real_name)
        setfield!(o, Symbol("real_name"), real_name)
        validate_property(ObjsUserProfileShort, Symbol("real_name_normalized"), real_name_normalized)
        setfield!(o, Symbol("real_name_normalized"), real_name_normalized)
        validate_property(ObjsUserProfileShort, Symbol("team"), team)
        setfield!(o, Symbol("team"), team)
        o
    end
end # type ObjsUserProfileShort

const _property_map_ObjsUserProfileShort = Dict{Symbol,Symbol}(Symbol("avatar_hash")=>Symbol("avatar_hash"), Symbol("display_name")=>Symbol("display_name"), Symbol("display_name_normalized")=>Symbol("display_name_normalized"), Symbol("image_72")=>Symbol("image_72"), Symbol("is_restricted")=>Symbol("is_restricted"), Symbol("is_ultra_restricted")=>Symbol("is_ultra_restricted"), Symbol("name")=>Symbol("name"), Symbol("real_name")=>Symbol("real_name"), Symbol("real_name_normalized")=>Symbol("real_name_normalized"), Symbol("team")=>Symbol("team"))
const _property_types_ObjsUserProfileShort = Dict{Symbol,String}(Symbol("avatar_hash")=>"String", Symbol("display_name")=>"String", Symbol("display_name_normalized")=>"String", Symbol("image_72")=>"String", Symbol("is_restricted")=>"Bool", Symbol("is_ultra_restricted")=>"Bool", Symbol("name")=>"String", Symbol("real_name")=>"String", Symbol("real_name_normalized")=>"String", Symbol("team")=>"DefsWorkspaceId")
Base.propertynames(::Type{ ObjsUserProfileShort }) = collect(keys(_property_map_ObjsUserProfileShort))
Swagger.property_type(::Type{ ObjsUserProfileShort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsUserProfileShort[name]))}
Swagger.field_name(::Type{ ObjsUserProfileShort }, property_name::Symbol) =  _property_map_ObjsUserProfileShort[property_name]

function check_required(o::ObjsUserProfileShort)
    (getproperty(o, Symbol("avatar_hash")) === nothing) && (return false)
    (getproperty(o, Symbol("display_name")) === nothing) && (return false)
    (getproperty(o, Symbol("image_72")) === nothing) && (return false)
    (getproperty(o, Symbol("is_restricted")) === nothing) && (return false)
    (getproperty(o, Symbol("is_ultra_restricted")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("real_name")) === nothing) && (return false)
    (getproperty(o, Symbol("team")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsUserProfileShort }, name::Symbol, val)
end


text:
## compare my chull implementation with that of the general qhull library

function testqhull(tar, non)
    xo, tc = sortscores(tar, non)
    θ, tc, nc, = binscores(xo, tc)
    pfa = [1, 1 - cumsum(nc)//length(non)]
    pmiss = [0, cumsum(tc)//length(tar)]
    @time ch, llr = chllr(tc, nc, xo, laplace=false, fast=false)
    @time myqh = rochull(pfa, pmiss)
    pfa, pmiss, find(ch), myqh
end

function test()
    ntar = rand(1:10000)
    nnon = rand(1:10000)
    tar = 1 + randn(ntar)
    non = -1 + randn(nnon)
    pfa, pmiss, a, b = testqhull(tar, non)
    if a != b
        println("Difference found:")
        if length(a) > length(b)
            println("qhull more than rochull")
            println(setdiff(a,b))
        else
            println("rochull more than qhull")
            println(setdiff(b,a))
        end
    end
    return pfa, pmiss, a, b
end

function prof()
    N = 1000_000
    tar = 2 + 2randn(N)
    non = -2 + 2randn(N)
    Profile.clear()
    @profile roc(tar, non)
    Profile.print(format=:flat)
end

ext: jl


text:
using MagmaThermoKinematics
using MagmaThermoKinematics.Diffusion3D
using MagmaThermoKinematics.Units
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
using Plots
using WriteVTK

# Initialize
@init_parallel_stencil(Threads, Float64, 3);    # initialize parallel stencil in 2D

#------------------------------------------------------------------------------------------
@views function MainCode_3D();
# Model parameters
W,L,H                   =   30km,30km,30km;         # Width, Length, Height
ρ                       =   2800;                   # Density
cp                      =   1050;                   # Heat capacity
k_rock, k_magma         =   1.5, 1.2;               # Thermal conductivity of host rock & magma
La                      =   350e3;                  # Latent heat J/kg/K
GeoT                    =   20.0/km;                # Geothermal gradient [K/km]
x_in,y_in,z_in          =   20km,20km,-15km;        # Center of dike [x,y,z coordinates]
W_in, H_in              =   5km,  0.5km;            # Width and thickness of dike
T_in                    =   900;                    # Intrusion temperature
InjectionInterval       =   0.1kyr;                 # Inject a new dike every X kyrs
maxTime                 =   15kyr;                  # Maximum simulation time
H_ran, W_ran            =   H*0.4, W*0.3;           # Size of domain in which we randomly place dikes and range of angles
DikeType                =   "ElasticDike"           # Type to be injected ("ElasticDike","SquareDike")

Nx, Ny, Nz              =   250, 250, 250;                          # Resolution
dx,dy,dz                =   W/(Nx-1),L/(Nx-1),H/(Nz-1);             # Grid size
κ                       =   k_rock./(ρ*cp);                         # Thermal diffusivity
dt                      =   min(dx^2,dy^2,dz^2)./κ/10;              # Stable timestep (required for explicit FD)
nt::Int64               =   floor(maxTime/dt);                      # Number of required timesteps
nTr_dike                =   1000;                                   # Number of tracers inserted per dike

# Array initializations
T                       =   @zeros(Nx,Ny,Nz);
K                       =   @ones(Nx,Ny, Nz)*k_rock;
Rho                     =   @ones(Nx,Ny, Nz)*ρ;
Cp                      =   @ones(Nx,Ny, Nz)*cp;

# Work array initialization
Tnew, qx,qy,qz          =   @zeros(Nx,Ny,Nz),   @zeros(Nx-1,Ny,Nz), @zeros(Nx,Ny-1,Nz), @zeros(Nx,Ny,Nz-1)  # Thermal solver
Kx, Ky, Kz              =   @zeros(Nx-1,Ny,Nz), @zeros(Nx,Ny-1,Nz), @zeros(Nx,Ny,Nz-1)                      # Thermal conductivities
Phi_o, Phi, dPhi_dt     =   @zeros(Nx,Ny,Nz),   @zeros(Nx,  Ny,Nz), @zeros(Nx,Ny,  Nz)                      # Solid fraction

# Set up model geometry & initial T structure
x,y,z                   =   (0:Nx-1)*dx, (0:Ny-1)*dy, (-(Nz-1):0)*dz;           # 1D coordinate arrays
crd                     =   collect(Iterators.product(x,y,z))                   # Generate coordinates from 1D coordinate vectors
X,Y,Z                   =   (x->x[1]).(crd),(x->x[2]).(crd),(x->x[3]).(crd);    # Transfer coords to 3D arrays
Grid                    =   (x,y,z);                                            # Grid
dike                    =   Dike(W=W_in,H=H_in,Type=DikeType,T=T_in);           # "Reference" dike with given thickness,radius and T
Tracers                 =   StructArray{Tracer}(undef, 1)                       # Initialize tracers
T                       .=   -Z.*GeoT;                                          # Initial (linear) temperature profile
Phi, dPhi_dt            =   SolidFraction(T, Phi_o, dt);                        # Compute solid fraction

# Preparation of VTK/Paraview output
if isdir("viz3D_out")==false mkdir("viz3D_out") end; loadpath = "./viz3D_out/"; pvd = paraview_collection("Example3D");

time, dike_inj, InjectVol, Time_vec,Melt_Time = 0.0, 0.0, 0.0,zeros(nt,1),zeros(nt,1);
for it = 1:nt   # Time loop

    if floor(time/InjectionInterval)> dike_inj                                                  # Add new dike every X years
        dike_inj        =   floor(time/InjectionInterval)                                       # Keeps track on whether we injected already
        cen             =   [W/2.;L/2.;-H/2.] + rand(-0.5:1e-3:0.5, 3).*[W_ran;W_ran;H_ran];    # Randomly vary center of dike
        if cen[end]<-12km;  Angle_rand = [rand(80.0:0.1:100.0); rand(0:360)]                    # Dikes at depth
        else                Angle_rand = [rand(-10.0:0.1:10.0); rand(0:360)] end                # Sills at shallower depth
        dike            =   Dike(dike,Center=cen[:],Angle=Angle_rand);                          # Specify dike with random location/angle but fixed size
        Tracers, T, Vol =   InjectDike(Tracers, T, Grid, dike, nTr_dike);                       # Add dike, move hostrocks
        InjectVol       +=  Vol                                                                 # Keep track of injected volume
        println("Added new dike; total injected magma volume = $(round(InjectVol/km³,digits=2)) km³; rate Q=$(round(InjectVol/(time),digits=2)) m³/s")
    end
    Phi, dPhi_dt        =   SolidFraction(T, Phi_o, dt);                                        # Compute solid fraction
    K                   .=  Phi.*k_rock .+ (1 .- Phi).*k_magma;                                 # Thermal conductivity

    # Perform a diffusion step
    @parallel diffusion3D_step_varK!(Tnew, T, qx, qy, qz, K, Kx, Ky, Kz, Rho, Cp, dt, dx, dy, dz,  La, dPhi_dt);
    @parallel (1:size(T,2), 1:size(T,3)) bc3D_x!(Tnew);                                         # Set lateral boundary conditions (flux-free)
    @parallel (1:size(T,1), 1:size(T,3)) bc3D_y!(Tnew);                                         # Set lateral boundary conditions (flux-free)
    Tnew[:,:,1] .= GeoT*H; Tnew[:,:,end] .= 0.0;                                                # Bottom & top temperature (constant)

    Tracers         =   UpdateTracers(Tracers, Grid, Tnew, Phi);                                # Update info on tracers
    T, Tnew         =   Tnew, T;                                                                # Update temperature
    time            =   time + dt;                                                              # Keep track of evolved time
    Melt_Time[it]   =   sum( 1.0 .- Phi)/(Nx*Ny*Nz)                                             # Melt fraction in crust
    Time_vec[it]    =   time/kyr;                                                               # Vector with time
    println(" Timestep $it = $(round(time/kyr,digits=2)) kyrs")

    if mod(it,5)==0  # Visualisation
        Phi_melt        =   1.0 .- Phi;
        vtkfile = vtk_grid("./viz3D_out/ex3D_$(Int32(it+1e4))", Vector(x/km), Vector(y/km), Vector(z/km)) # 3-D vtk file
        vtkfile["Temperature"] = T; vtkfile["MeltFraction"] = Phi_melt;                         # Store fields in file
        outfiles = vtk_save(vtkfile); pvd[time/kyr] = vtkfile                                   # Save file & update pvd file
    end
end
vtk_save(pvd)   # save Example3D.pvd file, which you can open with, for example, paraview
return Time_vec, Melt_Time
end             # end of main function
Time_vec, Melt_Time = MainCode_3D(); # start the main code
plot(Time_vec, Melt_Time, xlabel="Time [kyrs]", ylabel="Fraction of crust that is molten", label=:none); png("Time_vs_Melt_Example3D") #Create plot


text:
using MaximinOPF
using PowerModels
using JuMP
using SCS, Ipopt
using Test

include("testcases.jl")
PowerModels.silence()

function run_instance(testcase, PMOption, maximin = true, is_conic = true)
    pm_data = PowerModels.parse_file( testcase["file"] )
    pm_form = PMOption
    pm_data["name"]=testcase["name"]
    pm_data["attacker_budget"] = testcase["attack_budget"] ###Adding another key and entry
    pm_data["inactive_branches"] = testcase["inactive_indices"] ###Adding another key and entry
    pm_data["protected_branches"] = testcase["protected_indices"] ###Adding another key and entry
    pm_data["pm_form"] = pm_form

    #Create JUMP Model
    if maximin
        model = MaximinOPF.MaximinOPFModel(pm_data, pm_form; enforce_int=false) # "SCS does not have mixed-integer support"
    else
        pm = MaximinOPF.MinimaxOPFModel(pm_data, pm_form)
        model = pm.model
    end

    if is_conic
        set_optimizer(model, optimizer_with_attributes(SCS.Optimizer, "verbose" => 0))
    else
        set_optimizer(model, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    end
    result = @elapsed JuMP.optimize!(model)
    status = JuMP.termination_status(model)
    println(" [$(testcase["name"])]: $(result) seconds, status $(status).")
    @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL, MOI.LOCALLY_SOLVED]

    # Note: We do not test writing files, because the functions are defined in MOI.

    # t_sdpa = @elapsed JuMP.write_to_file(model, "$(testcase["name"]).dat-s", format = MOI.FileFormats.FORMAT_SDPA)
    # println("  Wrote $(testcase["name"]).dat-s file: $(t_sdpa) seconds")
    # rm("$(testcase["name"]).dat-s")

    # t_cbf = @elapsed JuMP.write_to_file(model, "$(testcase["name"]).cbf", format = MOI.FileFormats.FORMAT_CBF)
    # println("  Wrote $(testcase["name"]).cbf file: $(t_cbf) seconds")
    # rm("$(testcase["name"]).cbf")
end

@testset "Maximin Model Tests" begin
    println("Maximin Models:")
    @testset "SOC Model Tests" begin
        println("Type: SOCWRConicPowerModel")
        for case in testcases
            run_instance(case, SOCWRConicPowerModel)
        end
    end
    @testset "SDP Model Tests" begin
        println("Type: SDPWRMPowerModel")
        for case in testcases
            run_instance(case, SDPWRMPowerModel)
        end
    end
    @testset "Sparse SDP Model Tests" begin
        println("Type: SparseSDPWRMPowerModel")
        for case in testcases
            run_instance(case, SparseSDPWRMPowerModel)
        end
    end
end

@testset "Minimax Model Tests" begin
    println("Minimax Models:")
    @testset "Nonconvex AC Models" begin
        @testset "ACP Model Tests" begin
            println("Type: ACPPowerModel")
            for case in testcases
                run_instance(case, ACPPowerModel, false, false)
            end
        end
        @testset "ACR Model Tests" begin
            println("Type: ACRPowerModel")
            for case in testcases
                run_instance(case, ACRPowerModel, false, false)
            end
        end
        @testset "ACT Model Tests" begin
            println("Type: ACTPowerModel")
            for case in testcases
                run_instance(case, ACTPowerModel, false, false)
            end
        end
    end
    @testset "Linear Approximation Models" begin
        @testset "DCP Model Tests" begin
            println("Type: DCPPowerModel")
            for case in testcases
                run_instance(case, DCPPowerModel, false, false)
            end
        end
        @testset "DCMP Model Tests" begin
            println("Type: DCMPPowerModel")
            for case in testcases
                run_instance(case, DCMPPowerModel, false, false)
            end
        end
        @testset "NFA Model Tests" begin
            println("Type: NFAPowerModel")
            for case in testcases
                run_instance(case, NFAPowerModel, false, false)
            end
        end
    end
    @testset "Quadratic Approximation Models" begin
        @testset "DCPLL Model Tests" begin
            println("Type: DCPLLPowerModel")
            for case in testcases
                run_instance(case, DCPLLPowerModel, false, false)
            end
        end
        @testset "LPACC Model Tests" begin
            println("Type: LPACCPowerModel")
            for case in testcases
                run_instance(case, LPACCPowerModel, false, false)
            end
        end
    end
    @testset "Quadratic Relaxation Models" begin
        @testset "SOCWR Model Tests" begin
            println("Type: SOCWRPowerModel")
            for case in testcases
                run_instance(case, SOCWRPowerModel, false, false)
            end
        end
        @testset "SOCBF Model Tests" begin
            println("Type: SOCBFPowerModel")
            for case in testcases
                run_instance(case, SOCBFPowerModel, false, false)
            end
        end
        @testset "QCRM Model Tests" begin
            println("Type: QCRMPowerModel")
            for case in testcases
                run_instance(case, QCRMPowerModel, false, false)
            end
        end
        @testset "QCLS Model Tests" begin
            println("Type: QCLSPowerModel")
            for case in testcases
                run_instance(case, QCLSPowerModel, false, false)
            end
        end
    end
    @testset "Quadratic Conic Models" begin
        @testset "SOCWR Model Tests" begin
            println("Type: SOCWRConicPowerModel")
            for case in testcases
                run_instance(case, SOCWRConicPowerModel, false)
            end
        end
        @testset "SOCBF Model Tests" begin
            println("Type: SOCBFConicPowerModel")
            for case in testcases
                run_instance(case, SOCBFConicPowerModel, false)
            end
        end
    end
    @testset "SDP Models" begin
        @testset "SDPWRM Model Tests" begin
            println("Type: SDPWRMPowerModel")
            for case in testcases
                run_instance(case, SDPWRMPowerModel, false)
            end
        end
        @testset "SparseSDPWRM Model Tests" begin
            println("Type: SparseSDPWRMPowerModel")
            for case in testcases
                run_instance(case, SparseSDPWRMPowerModel, false)
            end
        end
    end
end


text:
abstract type AbstractFpDataset{T, N} <: AbstractArray{T, N} end
mutable struct MetaData{T<:Real}
    startd::DateTime
    endd::DateTime
    times::Vector{DateTime}
    totparts::Int
    rellats::Vector{T}
    rellons::Vector{T}
end

mutable struct FlexpartOutput
    filename::String
    lons::Vector{Real}
    lats::Vector{Real}
    # ncdataset::Union{NCDatasets.NCDataset, Nothing}
    # dataset::Union{Array{Real}, Nothing}
    # ncvar::Union{NCDatasets.CFVariable, Nothing}
    # # remdim::Union{NamedTuple, Nothing}
    # selected::Union{Dict{Symbol, Any}, Nothing}
    metadata::MetaData
end

struct FpDataset{T, N} <: AbstractFpDataset{T, N}
    fpoutput::FlexpartOutput
    dataset::Array{T, N}
    varname::String
    selected::NamedTuple
end

# output_dir(fpdir::FlexpartDir) = joinpath(fpdir.path, fpdir[:output])
function FlexpartOutput(filename::String)
    filename = abspath(filename)
    lons, lats = mesh(filename)
    FlexpartOutput(filename, lons, lats, MetaData(filename))
end

FlexpartOutput(fpdir::FlexpartDir, name::String) = FlexpartOutput(joinpath(getdir(fpdir, :output), name))
FlexpartOutput(path::String, name::String) = FlexpartOutput(FlexpartDir(path), name)

function ncf_files(fpdir::FlexpartDir; onlynested=false)
    out_files = readdir(getdir(fpdir, :output))
    f = onlynested ? x -> occursin(".nc", x) && occursin("nest", x) :  x ->  occursin(".nc", x)
    files = filter(f, out_files)
    [joinpath(getdir(fpdir, :output), x) |> abspath for x in files]
end

ncf_files(path::String; onlynested=false) = ncf_files(FlexpartDir(path), onlynested=onlynested)

function MetaData(filename::String)
    NCDataset(filename, "r") do ds
        rellats, rellons = relloc(filename)
        parts = ds["RELPART"]
        MetaData(start_dt(filename), end_dt(filename), times_dt(filename), sum(parts), rellons, rellats)
    end
end

function opends(f::Function, args...)
    NCDataset(args[1], "r") do ds
        f(ds)
    end
end
is2d(ncvar::NCDatasets.CFVariable) = ncvar[:] isa Matrix

Base.size(fpds::AbstractFpDataset) = size(fpds.dataset)
Base.getindex(fpds::AbstractFpDataset, args...) = getindex(fpds.dataset, args...)
Base.setindex!(fpds::AbstractFpDataset, args...) = setindex!(fpds.dataset, args...)
function Base.show(io::IO, m::MIME"text/plain", fpds::AbstractFpDataset{T, N}) where {T, N}
    println(io, "Flexpart Output with variable: ", fpds.varname)
    println(io, "Selected dimensions: ", fpds.selected)
    print(io, "Dataset: ", typeof(fpds.dataset))
end
function Base.:+(fpds1::AbstractFpDataset{T, N}, fpds2::AbstractFpDataset{T, N}) where {T, N}
    addable(fpds1, fpds2)
    FpDataset(fpds1.fpoutput, fpds1.dataset + fpds2.dataset, fpds1.varname*"_added", fpds1.selected)
end
function Base.:-(fpds1::AbstractFpDataset{T, N}, fpds2::AbstractFpDataset{T, N}) where {T, N}
    addable(fpds1, fpds2)
    FpDataset(fpds1.fpoutput, fpds1.dataset - fpds2.dataset, fpds1.varname*"_diff", fpds1.selected)
end

function addable(fpds1::AbstractFpDataset{T, N}, fpds2::AbstractFpDataset{T, N}) where {T, N}
    fpds1.varname == fpds2.varname || error("Selected variables must be equal")
    Flexpart.completedim(fpds1) |> keys == Flexpart.completedim(fpds2) |> keys || error("Same dimensions must be selected")
    true
end

function variables(ncdataset::NCDataset)
    variables = keys(ncdataset)
end

function variables2d(ncdataset::NCDataset)
    vars = variables(ncdataset)
    lons, lats = mesh(ncdataset)
    inds = findall( x ->
        try
            size(ncdataset[x])[1:2] == (length(lons), length(lats))
        catch
            false
        end,
        vars
        )
    vars[inds]
end

function variables(fpoutput::FlexpartOutput)
    opends(fpoutput.filename) do ds
        variables(ds)
    end
end
function variables2d(fpoutput::FlexpartOutput)
    opends(fpoutput.filename) do ds
        variables2d(ds)
    end
end

function getvar(fpoutput::FlexpartOutput, varname)
    opends(fpoutput.filename) do ds
        ds[varname][:]
    end
end

function remdim(output::FlexpartOutput, varname::String)
    opends(output.filename) do ds
        ncvar = ds[varname]
        remdim(ncvar)
    end
end
function remdim(ncvar::NCDatasets.CFVariable)
    names = NCDatasets.dimnames(ncvar)
    sizes = NCDatasets.dimsize(ncvar)
    NamedTuple{Symbol.(names[3:end])}(sizes)
end

attrib(ds::NCDatasets.NCDataset) = Dict(ds.attrib)
attrib(ds::NCDatasets.CFVariable) = push!(Dict(ds.attrib), "name" => NCDatasets.name(ds))
function attrib(fpoutput::FlexpartOutput)
    opends(fpoutput.filename) do ds
        attrib(ds)
    end
end
function attrib(fpoutput::FlexpartOutput, varname::String)
    opends(fpoutput.filename) do ds
        attrib(ds[varname])
    end
end

function alldims(ds::NCDatasets.NCDataset, ncvar::NCDatasets.CFVariable) :: Dict{Symbol, Vector{Any}}
    dims = remdim(ncvar)
    d = Dict{Symbol,  Vector{Any}}()
    for k in keys(dims)
        try
            push!(d, k => ds[k][:])
        catch
            push!(d, k => 1:dims[k] |> collect)
        end
    end
    d
end
function alldims(fpoutput::FlexpartOutput, varname::String)
    opends(fpoutput.filename) do ds
        alldims(ds, ds[varname])
    end
end

function getspatial(ncds::NCDataset, ncvar::NCDatasets.CFVariable, dims::NamedTuple = NamedTuple())
    lons, lats = mesh(ncds)
    try
        @assert size(ncvar)[1:2] == (length(lons), length(lats))
    catch
        throw(ArgumentError("The selected field is not a spatial field"))
    end
    rd = remdim(ncvar)
    indices = []
    for k in keys(rd)
        if k in keys(dims)
            push!(indices, dims[k])
        else
            push!(indices, :)
        end
    end
    nt = NamedTuple{keys(rd)}(indices)
    ncvar[:,:,nt...]
end

function select(fpoutput::FlexpartOutput, varname::String, dims::NamedTuple = NamedTuple())
    opends(fpoutput.filename) do ds
        field = getspatial(ds, ds[varname], dims)
        FpDataset(fpoutput, field, varname, dims)
    end
end

function select(fpoutput::FlexpartOutput, varname::String, dims::Dict) :: FpDataset
    adims = alldims(fpoutput, varname)
    indices = []
    for (k,v) in dims
        v = try
            v = try
                parse(Float64, v)
            catch
                v
            end
            d = adims[k][1]
            v, d = try
                convert(Float64, v), convert(Float64, d)
            catch
                v, d
            end
            typev = typeof(v)
            typed = typeof(d)
            typev == typed ?
                v :
                parse_date([(v, "yyyymmddTHHMMSS"), (v[1:19], "yyyy-mm-ddTHH:MM:SS")])
        catch
            error("Wrong type to be selected")
        end

        inds = try
            findfirst(x -> x == v, adims[k])
        catch
            error("Some selected variables don't exist")
        end
        push!(indices, inds)
    end
    ndims = NamedTuple{Tuple(keys(dims) |> collect)}(indices)
    select(fpoutput, varname, ndims)
end

function parse_date(args::Vector{Tuple{String, String}})
    tried = []
    err = nothing
    for arg in args
        t = try
            DateTime(arg...)
        catch ex
            err = ex
            nothing
        end
        !isnothing(t) && push!(tried, t)
    end

    length(tried) == 0 && throw(err)
    tried[1]
end

function selected(fpds::FpDataset)
    sel = fpds.selected
    adims = alldims(fpds.fpoutput, fpds.varname)

    Dict(k => adims[k][v] for (k, v) in zip(keys(sel), sel))
end

mesh(fpoutput::FlexpartOutput) = mesh(fpoutput.filename)
function mesh(file)
    opends(file) do ds
        mesh(ds)
    end
end

function mesh(ncdataset::NCDataset)
    lon = ncdataset["longitude"]; lat = ncdataset["latitude"]
    lon[:], lat[:]
end

function deltamesh(lons, lats)
    dxs = lons[2:end] - lons[1:end-1]
    dys = lats[2:end] - lats[1:end-1]

    dx = unique(round.(dxs, digits=5))
    dy = unique(round.(dys, digits=5))

    if (length(dx) != 1) || (length(dy) != 1)
        error("mesh is not uniform")
    end

    dx[1], dy[1]
end

function areamesh(lons, lats)
    min_lon = minimum(lons)
    max_lon = maximum(lons)
    if min_lon > 180 || max_lon > 180
        min_lon -= 360
        max_lon -= 360
    end
    if min_lon < -180 || max_lon < -180
        min_lon += 360
        max_lon += 360
    end
    [maximum(lats), min_lon, minimum(lats), max_lon]
end

function filter_fields(lon, lat, field)
    if size(field) != (length(lon), length(lat)) error("dimension mismatch : size(field) == (length(lon), length(lat))") end
    mask = (!).(isapprox.(0., field))

    mg_lon = lon .* ones(length(lat))'
    mg_lat = lat' .* ones(length(lon))

    return mg_lon[mask], mg_lat[mask], field[mask]
end

function relloc(file)
    NCDataset(file, "r") do ds
        rellats = ds["RELLAT1"] |> unique
        rellons = ds["RELLNG1"] |> unique
        rellons, rellats
    end
end

function faverage(alltimes::Array{<:Real, 3})

    added = zeros(eltype(alltimes), size(alltimes)[1:2])
    N = length(alltimes[1, 1, :])
    for i=1:N
        added += alltimes[:, :, i]
    end

    return added ./ N
end

# function faverage(output::FlexpartOutput)
#     hasselection(output) && length(size(output.dataset)) == 3 || error("The dataset must be 3D")
#     faverage(output.dataset)
# end

function daily_average(alltimes::Array{<:Real, 3}, times::Vector)
    sdays = split_days(times)
    daily_av = zeros(eltype(alltimes), (size(alltimes)[1:2]..., length(sdays)))
    for (index, day_index) in enumerate(sdays)
        days = alltimes[:, :, day_index];
        daily_av[:, :, index] = faverage(days)
    end
    return daily_av
end

function daily_average(fpds::FpDataset)
    fpds.dataset isa Array{T, 3} where T || error("The dataset must be 3D")
    #TODO: throw error if selected is not time
    daily_average(fpds.dataset, alldims(fpds.fpoutput, fpds.varname)[:time])
end

completedim(fpds::FpDataset) = filter(e -> !(e.second isa Vector), selected(fpds))

# function write_daily_average(fpds::FpDataset; copy = true)
#     name = fpds.varname*"_daily_av"
#     filename = fpds.fpoutput.filename
#     if copy
#         newfn = split(basename(filename), ".")
#         newfn = joinpath(dirname(filename), newfn[1] * "_daily_av" * "." * newfn[2])
#         filename = cp(filename, newfn, force=true)
#     end
#     daily_ds, days = daily_average(fpds)

#     attr = attrib(fpds.fpoutput, fpds.varname)
#     toadd = completedim(fpds)

#     NCDataset(filename, "a") do ds
#         defVar(ds, "day", days, ("day",))

#         defVar(ds,name,convert.(Float64, daily_ds),("lon","lat","day"), attrib = merge(symb2string.((attr, toadd))...))
#     end
#     filename
# end

function write_daily_average(fpds::FpDataset; copy = true)
    name = fpds.varname*"_daily_av"
    filename = fpds.fpoutput.filename
    if copy
        newfn = split(basename(filename), ".")
        newfn = joinpath(dirname(filename), newfn[1] * "_daily_av" * "." * newfn[2])
        filename = cp(filename, newfn, force=true)
    end

    times = alldims(fpds.fpoutput, fpds.varname)[:time]
    s = size(fpds)
    s = (s[1:3]..., s[5:end]...)
    sdays = split_days(times)
    ndays = sdays|> length
    days = DateTime.(times[[s[1] for s in sdays]])
    new_dataset = zeros(s..., ndays)
    for height in 1:s[3], pointspec in 1:s[4], nageclass in 1:s[5]
        new_dataset[:,:,height, pointspec, nageclass, :] = daily_average(fpds[:,:,height, :, pointspec, nageclass], times)
    end

    attr = attrib(fpds.fpoutput, fpds.varname)
    toadd = completedim(fpds)

    NCDataset(filename, "a") do ds
        defVar(ds, "day", days, ("day",))

        defVar(ds,name,convert.(Float64, new_dataset),("lon","lat", "height", "pointspec", "nageclass", "day"), attrib = merge(symb2string.((attr, toadd))...))
    end
    filename
end

# macro

function write_worst_day(output::FlexpartOutput; copy=true)
    (hasselection(output) && length(size(output.dataset))) == 3 || error("The dataset must be 3D")
    name = "worst_day"
    filename = output.filename
    if copy
        newfn = split(output.filename, ".")
        newfn = newfn[1] * "_worst_day" * "_copy." * newfn[2]
        filename = cp(output.filename, newfn, force=true)
    end
    dsworst, indworst = findworst_day(output.dataset)
    dayworst = output.ncdataset["day"][:][indworst]

    attr = attrib(output.ncvar)
    toadd = Dict("day" => Dates.format(dayworst, "yyyymmdd"))
    copy || close(output.ncdataset)
    ds = NCDataset(filename, "a")

    defVar(ds,name,convert.(Float64, dsworst),("lon","lat"), attrib = merge(symb2string.((attr, toadd))...))
    close(ds)
    FlexpartOutput(filename)
end

symb2string(d::Dict) = Dict(String(k) => v for (k,v) in d)

function sum_abs(field)
    sum(abs, field)
end

function findworst_day(fields::AbstractArray)
    sums = [sum_abs(f) for f in eachslice(fields, dims=3)]
    _, ind = findmax(sums)
    fields[:, :, ind], ind
end

function worst_day(fields::AbstractArray)
    f, _ = findworst_day(fields)
    f
end

function start_dt(file)
    NCDataset(file, "r") do ds
        globatt = ds.attrib
        Dates.DateTime(globatt["ibdate"]*globatt["ibtime"], dateformat"yyyymmddHHMMSS")
    end
end

function end_dt(file)
    NCDataset(file, "r") do ds
        globatt = ds.attrib
        Dates.DateTime(globatt["iedate"]*globatt["ietime"], dateformat"yyyymmddHHMMSS")
    end
end

function split_days(times::Vector)
    split = Vector{Vector{Int}}()
    cur_day = Dates.day(times[1])
    day_inds = Vector{Int}()
    for (index, time) in enumerate(times)
        day = Dates.day(time)
        if day == cur_day
            push!(day_inds, index)
        else
            push!(split, day_inds)
            day_inds = Vector{Int}()
            push!(day_inds, index)
            cur_day = day
        end
    end
    return split
end

function times_dt(file)::Vector{DateTime}
    NCDataset(file, "r") do ds
        ds["time"][:]
    end
end


text:
# http://juliaimages.github.io/latest/function_reference.html#ImageCore.scalesigned

using ImageCore, TestImages
import ImageView
img = testimage("cameraman")
img1 = img[10:500, 10:500]
img2 = img[12:502, 10:500]
dimg = float(img1)-float(img2)
const ss = scalesigned(1)   # replace 1 with whatever scaling you want to encode
const cs = colorsigned()
imshow_now(dimg, scalei=x->cs(ss(x)))


text:
import GitLab
using Base.Test
myauth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode your access tokens!
println("Authentication successful")
options = Dict("private_token" => myauth.token)

myrepo = GitLab.repo_by_name("TestProject1"; headers=options)

branches = first(GitLab.branches(myrepo; params=options))

@test GitLab.name(branches[1]) == "branch1"
## @show branches

println("Branches - Tests Done !!!")


text:
function modalmethod(M,K,zm,F,nm,dt)

	# eigenvalues and eigenectors
	D,V = eig(K,M)

	# modal matrices
	Mm = V'*M*V
	Km = V'*K*V
	Fm = V'*F

	ndof,nt = size(F)

	qm=zeros(ndof,nt)
	q1m=zeros(ndof,nt)
	q2m=zeros(ndof,nt)
	for modo in 1:nm
		  m = Mm[modo,modo]
		  k = Km[modo,modo]
		  c = 2*sqrt(m*k)*zm[modo]
		  f = Fm[modo,:]

		  ym,vm,am = state1(m,c,k,f,dt,0.,0.)
		  qm[modo,:] = ym
		  q1m[modo,:] = vm
		  q2m[modo,:] = am
	end

	# deshacemos el cambio
	q = V*qm
	q1 = V*q1m
	q2 = V*q2m

	return q,q1,q2

end

# ========================================================
function modalmethod_test()
	include("FEMmodel.jl")
	E = 2.1e11
	rho = 7850.
	b = 2.
	h = 0.20
	L = 20.
	ne = 18

	I=1/12*b*h^3
	A=b*h

	k1=100*E*I
	k2=100*E*I

	(K,M) = FEMmodel(E,rho,A,I,L,ne,k1,k2)

	ndof = size(K,1)
	zm = 0.02*ones(ndof) # modal dampings
	F = randn(ndof,100)

	nm=3

	dt = 0.1

	q,q1,q2 = modalmethod(M,K,zm,F,nm,dt)

	return q,q1,q2

end


text:
#========================================================================

	Steady state solution for housing problem with non-convexities

	Translated julia code based on code from Ben Moll's Website:
		http://www.princeton.edu/%7Emoll/HACTproject.htm

		For Julia 1.0.0
========================================================================#
using LinearAlgebra, SparseArrays, Plots

# Model parameters
γ = 2 # CRRA utility
r = 0.035 #return on asset
α = 1/3 # parameter for cobb-douglas production function
η = 0.2 # for our benefit function
p = 1 #price
ϕ = 2 # coefficient of a in the budget constraint
hmin = 2.3 #lower threshold for houses
δ = 0.05 # captal depreciation
ρ = 0.05 # discount rate

# Parameters for productivity and the poisson process q
z_1 = .1
z_2 = .135
z = [z_1 z_2]
λ_1 = .5 # prob of transitioning from state 1 to 2
λ_2 = .5  # state 2 to 1
λ = [λ_1 λ_2]
z_ave = (z_1*λ_2 + z_2*λ_1)/(λ_1 + λ_2) # Average value of z

H =500 # Number of elements in grid space for a
a_min = 0.0
a_max = 3.0

# Grid spaces for variables
a = LinRange(a_min, a_max, H)
a = convert(Array,a)
da = (a_max-a_min)/(H-1)

aa = [a a]
zz = ones(H,1)*z #check later

# Construct the matrix to summarize the evolution of z over time
# this matrix will not change over the finite differences and is based on the poisson process

Aswitch = [-sparse(I,H,H)*λ[1] sparse(I,H,H)*λ[2] ; sparse(I,H,H)*λ[1] -sparse(I,H,H)*λ[2]]

# Paramters for the simulation
maxit = 120
crit = 10^(-10)
Δ = 1000 #step size for HJB

# Spaces for the finite difference terms
Vaf, Vab, c= [zeros(H,2) for i in 1:3]

dist =[]
V_n=[]

h = min.((α*η/(r*p))^(1/(1-α)) + hmin, ϕ*aa/p)
h = h.*(h.>=hmin)
f = η*(max.(h.-hmin,0)).^α - r*p*h

v0 = (zz + r*aa).^(1-γ)/(1-γ)/ρ

global v = v0

# The main loop, this makes sure markets clear

for n in 1:maxit
	V=v
	push!(V_n, V)
	# Forward difference
	Vaf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
	Vaf[H,:] = (z' .+ f[H,:] .+ r*a_max).^(-γ)

	# Backward difference
	Vab[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
	Vab[1,:] = (z' .+ f[1,:].+ r*a_min).^(-γ)

	# caculate consumption and savings
	# First the forward difference case
	cf = (max.(Vaf,10^(-10))).^(-1/γ)
	ssf = zz + f + r.*aa -cf
	Hf = cf.^(1-γ)/(1/γ) + Vaf.*ssf

	# consumption and savings for the backward difference
	cb = (max.(Vab,10^(-10))).^(-1/γ)
	ssb = zz + f + r.*aa -cb
	Hb = cb.^(1-γ)/(1-γ) + Vab.*ssb

	#Consumption and V' at the steady state
	c0 = zz + f + r.*aa
	Va0 = c0.^(-γ)

	# Upwind scheme chooses between the forward or backward difference
	Ieither = (1 .- (ssf.>0)) .* (1 .- (ssb.<0))
	Iunique = (ssb.<0).*(1 .- (ssf.>0)) + (1 .- (ssb.<0)).*(ssf.>0)
	Iboth = (ssb.<0).*(ssf.>0)
    If= Iunique.*(ssf .> 0) + Iboth.*(Hf.>=Hb)  #positive drift → forward difference
    Ib = Iunique.*(ssb .< 0) + Iboth.*(Hb.>=Hf) #negative drift → backward difference
    I0 = Ieither  #at steady state

	global c = cf.*If + cb.*Ib + c0.*I0
    u = (c.^(1-γ))/(1-γ)

	X = -min.(ssb,0)/da
	Y = -max.(ssf,0)/da + min.(ssb,0)/da
	Z = max.(ssf,0)/da

	A1 = spdiagm(0=>Y[:,1], -1 => X[2:H,1], 1=> Z[1:H-1,1])
	A2 = spdiagm(0=>Y[:,2], -1 => X[2:H,2], 1=> Z[1:H-1,2])
	AA = [A1 spzeros(H,H); spzeros(H,H) A2]

	global A= AA + Aswitch

	B = (1/Δ + ρ)*sparse(I,2*H,2*H)-A

	u_stacked = [u[:,1]; u[:,2]]
	V_stacked = [V[:,1]; V[:,2]]

	b = u_stacked + V_stacked/Δ

	V_stacked = B\b

	global V = [V_stacked[1:H] V_stacked[H+1:2*H]]

	V_change = V-v

	global v = V

	push!(dist, findmax(V_change)[1])
		if dist[n] < crit
			println("Value Function Converge Iteration = $(n)")
			break
		end
end

# Now solve the Fokker-Planck equation
AT = A'
gg =[] #create empty cell to store g

# Create an initial distribution
gg0 = ones(2*H,1)
g_sum = gg0'*ones(2*H,1)*da
gg0 = gg0./g_sum #normalize the initial distribution

push!(gg,gg0)
N = 1000
dt = 10
g_dist = zeros(N,1)

for n in 1:N
	push!(gg,(sparse(I,2*H,2*H)-AT*dt)\gg[n])
	g_dist[n] = findmax(abs.(gg[n+1]-gg[n]))[1]
end

g = [gg[N][1:H] gg[N][H+1:2*H]]
adot = zz + f + r.*aa - c
astar = p*hmin/ϕ

obj, index = findmin(abs.(astar.-a))

c2 = c - η*max.(h.-hmin,0).^α

amax1 = a_max
amin1 = -0.1

#Plot savings, output has some rough parts unlike matlab counter part...
plot(a, adot[:,1], label="\$s_{1}(a)\$", legend=:topright,
	xlims=(amin1,amax1),color=:blue,
	ylabel="Savings", xlabel="Wealth")
plot!(a,adot[:,2], label="\$s_{2}(a)\$", color=:red)
plot!(LinRange(amin1,amax1,H),zeros(H,1), label="",line=:dash, color=:black)
plot!(a_min.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)
plot!(astar.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)
png("savings_plot")

# Housing consumption, less smooth as well
plot(a, c2[:,1], label="\$c_{1}(a)\$",
	xlims=(amin1,amax1),color=:blue,
	ylabel="Housing Consumption", xlabel="Wealth", legend=:bottomright)
plot!(a,c2[:,2], label="\$c_{2}(a)\$", color=:red)
plot!(a_min.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)
plot!(astar.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)

# Housing
plot(a, h[:,1], label="\$h_{1}(a)\$",
	xlims=(amin1,amax1),color=:blue,
	ylabel="Housing", xlabel="Wealth", legend=:false)
plot!(a,h[:,2], label="\$h_{2}(a)\$", color=:red, line=:dashdot)
plot!(a_min.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)
plot!(astar.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)

# Benefit
plot(a, f[:,1],
	xlims=(amin1,amax1),color=:blue,
	ylabel="Pecuniary Benefit from Housing", xlabel="Wealth", legend=:false)
plot!(a,f[:,2], color=:red, line=:dashdot)
plot!(a_min.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)
plot!(astar.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)

# Value Function
plot(a, v[:,1], label="\$v_{1}(a)\$",
	xlims=(amin1,amax1),color=:blue,
	ylabel="Value Function", xlabel="Wealth", legend=:bottomright)
plot!(a,v[:,2], label="\$v_{2}(a)\$", color=:red, line=:dashdot)
plot!(a_min.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)
plot!(astar.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)

# Value Function
plot(a, g[:,1], label="\$g_{1}(a)\$",
	xlims=(amin1,amax1), ylim=(0,3.5),color=:blue,
	ylabel="Densities", xlabel="Wealth", legend=:bottomright)
plot!(a,g[:,2], label="\$g_{2}(a)\$", color=:red, line=:dashdot)
plot!(a_min.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)
plot!(astar.*ones(H,1), label="",line=:dash, color=:black, seriestype=:vline)


text:
"""
Convert result from evaluating single policy with eval_policy into a single row
"""
function res2row(res)
    row_dfs = []

    # process pp_params into df with single row
    row = DataFrame(res[:pp_params])
    push!(row_dfs, row)

    # process agg: move evry row from "description" into single dataframe of one row
    agg = select(res[:agg], [:variable, :mean, :median, :min, :max])
    for row in eachrow(agg)
        v = string(row.variable)
        row = DataFrame(row[2:end])
        rename!(s-> s*"_$v", row)
        push!(row_dfs, row)
    end

    # compact dfs into single row
    df = hcat(row_dfs...)

    # add method and its order
    method = string(first(res[:results]).name)*"_"*res[:method_info]
    insertcols!(df, 1, :method=>method)

    # add objective
    if haskey(res[:pp_params], :objective)
        objective = string(res[:pp_params][:objective])
    else
        objective = :revenue
    end
    if !("objective" in names(df))
#         df.objective = objective
        insertcols!(df, 2, :objective=>objective)
    end

    # add n_traces
    insertcols!(df, 2, :N=>res[:N])

    # add pp_params_string (without :objective)
    insertcols!(df,2, :pp_params_str=> string(delete!(Dict(res[:pp_params]), :objective)))

end

function folder_report(res_folder::String; raw_result_array=false)
    res_files = readdir(res_folder)
    rows = []
    raw = []
    for res_file in res_files
        if res_file[end-4:end]==".bson"
            res = BSON.load(joinpath(res_folder,res_file))
            row = res2row(res)
            push!(rows, row)
            raw_result_array && push!(raw, res)
        end
    end
    results = vcat(rows...)
    (;results=results, raw=raw)
end

function format_result_table(results::DataFrame;N=10)

    columns = [:method, :pp_params_str, :objective, :mean_r, :mean_u, :mean_bytes, :mean_time]

    df10 = filter(:N => N->N==N, results)
    gps = groupby(df10, [:objective])

    restable = outerjoin(
        [select(gr, columns) for gr in gps]...;
        on=[:method, :pp_params_str], makeunique=true, renamecols = "_obj_r" => "_obj_u")

    # restable = restable[:, [1,2,8,9,10,11,12,18,19,20,21,22]]
    sort!(restable, [:pp_params_str, order(:method, rev=true)])
end


text:
function testA()
  a = []
  for i in 1:1000
    push!(a, repeat("a", i))
  end
  return a
end

function testB()
  b = []
  for i in 1:10000
    push!(b, repeat("b", i))
  end
  return b
end

function tester()
  for i in 1:100
    println(length(testA()))
    println(length(testB()))
  end
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using ACE, JuLIP, BenchmarkTools

#---

fname = "/Users/ortner/Dropbox/Public/SHIPPots/Si_split_1.2_reg_test_v05.json"
D = load_dict(fname)
V = ACE.Import.import_pipot_v05(D)

at = bulk(:Si, cubic=true) * 3;
@show length(at)

#---

println("length(PIBasis) = ")
@show length(V.pibasis.inner[1])

#---

@info("Standard Evaluator")

@btime energy($V, $at)
@btime forces($V, $at)

print("time/atom (energy): ");
println(" $((@timed energy(V, at)).time / length(at) * 1000) ms")

print("time/atom (forces): ");
println(" $((@timed forces(V, at)).time / length(at) * 1000) ms")

#---

@info("Graph Evaluator")
println("(it takes a long time to generate the graph - sorry!")
Vgr = ACE.GraphPIPot(V)
println("Number of nodes")
@show length(Vgr.dags[1].nodes)

@btime energy($Vgr, $at)
@btime forces($Vgr, $at)
print("time/atom (energy): ");
println(" $((@timed energy(Vgr, at)).time / length(at) * 1000) ms")
print("time/atom (forces): ");
println(" $((@timed forces(Vgr, at)).time / length(at) * 1000) ms")


text:
#==============================================================================#
#
#               Routines for sampling trajectories from the
#               diffusion laws defined by this package.
#
#==============================================================================#
"""
    abstract type AbstractSDESolver end

Supertype of flags indicating ODE solvers
"""
abstract type AbstractSDESolver end

"""
    struct EulerMaruyama <: AbstractSDESolver end

Flag for indicating use of the Euler-Maruyama scheme for sampling diffusions.
IMPORTANT: this is the only diffusion path sampler implemented in this package.
There are no plans for implementing any other SDE solver in the forseeable
future.
"""
struct EulerMaruyama <: AbstractSDESolver end

"""
    struct Wiener{D}
    end

A struct defining the Wiener process. `D` indicates the dimension of the process
if it cannot be inferred from the DataType of the Trajectory. If the dimension
can be inferred then it takes precedence over the value of `D`.
"""
struct Wiener{D}
    Wiener{D}() where D = new{D}()
    Wiener(D::Integer) = new{D}()
    Wiener() = new{false}()
end
wiener() = Wiener()

ismutable(el) = ismutable(typeof(el))
ismutable(::Type) = Val(false)
ismutable(::Type{<:Array}) = Val(true)

"""
    Base.zero(K::Type, D, ::Val)

If `K` is a mutable type, then create `zeros` of dimension `D` and entries with
types `eltype(K)`. Otherwise, calls regular zero(`K`).
"""
Base.zero

Base.zero(K::Type, D, ::Val{true}) = zeros(eltype(K), D)
Base.zero(K::Type, D, ::Val{false}) = zero(K)

"""
    Trajectories.trajectory(tt, v::Type, D::Number)

Create a Trajectory with mutable states of dimension `D` along the time
collection `tt`.
"""
function Trajectories.trajectory(tt, v::Type, D::Number)
    trajectory(collect(tt), [zeros(eltype(v),D) for _ in tt])
end

"""
    Trajectories.trajectory(tt, v::Type)

Create a Trajectory with immutable states along the time collection `tt`.
"""
function Trajectories.trajectory(tt, v::Type)
    trajectory(collect(tt), zeros(v, length(tt)))
end
Trajectories.trajectory(tt, v::Type, D, ::Val{true}) = trajectory(tt, v, D)
Trajectories.trajectory(tt, v::Type, D, ::Val{false}) = trajectory(tt, v)

"""
    Random.rand!(
        path::Trajectory{T,Vector{K}},
        w::Wiener{D},
        y1=zero(K,D,ismutable(K))
    ) where {T,K,D}

Samples Wiener process started from `y1` and saves the data in `path`. Uses
a default random number generator.
"""
function Random.rand!(
        path::Trajectory{T,Vector{K}},
        w::Wiener{D},
        y1=zero(K,D,ismutable(K))
    ) where {T,K,D}
    rand!(Random.GLOBAL_RNG, path, w, y1)
end

"""
    Random.rand!(
        rng::Random.AbstractRNG,
        path::Trajectory{T,Vector{K}},
        ::Wiener,
        y1=zero(K)
    ) where {T,K}

Samples Wiener process with immutable states, started from `y1` and saves the
data in `path`. `rng` is used as a random number generator.
"""
function Random.rand!(
        rng::Random.AbstractRNG,
        path::Trajectory{T,Vector{K}},
        ::Wiener,
        y1=zero(K)
    ) where {T,K}
    yy, tt = path.x, path.t
    N = length(path)
    yy[1] = y1
    for i in 2:N
        rootdt = sqrt(tt[i] - tt[i-1])
        yy[i] = yy[i-1] + rootdt*randn(rng, K)
    end
    path
end

"""
    Random.rand!(
        rng::Random.AbstractRNG,
        path::Trajectory{T,Vector{Vector{K}}},
        ::Wiener{D},
        y1=zeros(K,D)
    ) where {K,T,D}

Samples Wiener process with mutable states, started from `y1` and saves the data
in `path`. `rng` is used as a random number generator.
"""
function Random.rand!(
        rng::Random.AbstractRNG,
        path::Trajectory{T,Vector{Vector{K}}},
        ::Wiener{D},
        y1=zeros(K,D)
    ) where {K,T,D}
    yy, tt = path.x, path.t
    N = length(path)
    yy[1] = y1
    for i in 2:N
        rootdt = sqrt(tt[i] - tt[i-1])
        for j in 1:D
            yy[i][j] = yy[i-1][j] + rootdt*randn(rng, K)
        end
    end
    path
end

"""
    Base.rand(tt, w::Wiener, y1)

Samples Wiener process on `tt`, started from `y1` and returns a new object with
a sampled trajectory.
"""
Base.rand(tt, w::Wiener, y1) = rand(Random.GLOBAL_RNG, tt, w, y1)

function Base.rand(
        rng::Random.AbstractRNG,
        tt,
        w::Wiener{D},
        y1::K
    ) where {D,K}
    path = trajectory(tt, K, D, ismutable(K))
    rand!(rng, path, w, y1)
end

"""
    solve!(XX, WW, P, y1)

Compute a trajectory, started from `y1` and following the diffusion law `P`,
from the sampled Wiener process `WW`. Save the sampled path in `XX`. Return
prematurely with a `false` massage if the numerical scheme has led to the solver
violating the state-space restrictions.
"""
solve!(XX, WW, P, y1) = solve!(EulerMaruyama(), XX, WW, P, y1)

"""
    solve!(XX, WW, P, y1, buffer)

Same as `solve!(XX, WW, P, y1)`, but additionally provides a pre-allocated
buffer for performing in-place computations.
"""
solve!(XX, WW, P, y1, buffer) = solve!(EulerMaruyama(), XX, WW, P, y1, buffer)

function solve!(
        ::EulerMaruyama,
        XX::Trajectory{T,Vector{KX}},
        WW::Trajectory{T,Vector{KW}},
        P,
        y1::KX,
    ) where {KX,KW,T}
    yy, ww, tt = XX.x, WW.x, XX.t
    N = length(XX)

    yy[1] = y1
    for i in 2:N
        dt = tt[i] - tt[i-1]
        dW = ww[i] - ww[i-1]
        yy[i] = yy[i-1] + _b((tt[i-1], i-1), yy[i-1], P)*dt + _σ((tt[i-1], i-1), yy[i-1], P)*dW
        bound_satisfied(P, yy[i]) || return false
    end
    true
end

function solve!(
        ::EulerMaruyama,
        XX::Trajectory{T,Vector{Vector{K}}},
        WW::Trajectory{T,Vector{Vector{K}}},
        P,
        y1::Vector{K},
        buffer=StandardEulerBuffer{K}(dimensions(P)...)
    ) where {K,T}
    yy, ww, tt = XX.x, WW.x, XX.t
    N = length(XX)
    dim_proc, dim_wiener = dimension(P)

    yy[1] = y1
    dt = tt[2] - tt[1]
    y = yy[1]

    for i in 2:N
        dt = tt[i] - tt[i-1]
        y = yy[i-1]
        _b!(buffer, (tt[i-1], i-1), y, P)
        _σ!(buffer, (tt[i-1], i-1), y, P)
        for j in 1:dim_wiener
            buffer.dW[j] = ww[i][j] - ww[i-1][j]
        end
        mul!(yy[i], buffer.σ, buffer.dW)
        for j in 1:dim_proc
            yy[i][j] += y[j] + buffer.b[j]*dt
        end
        bound_satisfied(P, yy[i]) || return false
    end
    true
end


text:
module ImageSmooth

using ImageFiltering, FFTW
using ImageBase
using ImageBase.ImageCore: GenericGrayImage, GenericImage
using OffsetArrays

#TODO: port SmoothAPI to ImagesAPI
include("SmoothAPI/SmoothAPI.jl")
import .SmoothAPI: AbstractImageSmoothAlgorithm,
                   smooth, smooth!

# Smooth algorithms

include("algorithms/l0_smooth.jl") # L0Smooth

export
    # generic API
    smooth, smooth!,

    # Algorithms
    L0Smooth

end # module ImageSmooth


text:
using Test
using CliOptions

include("testutils.jl")

@testset "CliOptionSpec" begin
    @testset "ctor(); program = $(repr(v[1]))" for v in [
        ("", "PROGRAM"),
        ("a", "a"),
    ]
        program, expected = v
        spec = CliOptionSpec(
            FlagOption("-a"),
            program = program,
            onerror = error,
        )
        @test spec.program == expected
    end

    @testset "show(io, x)" begin
        let spec = CliOptionSpec(
            FlagOption("-a"),
            CounterOption("-b"),
            Option("-c"),
            Positional("d"),
            onerror = error,
        )
            buf = IOBuffer()
            show(buf, spec)
            @test String(take!(buf)) == "CliOptionSpec(" * join([
                "FlagOption(:a)",
                "CounterOption(:b)",
                "Option(:c)",
                "Positional(:d)",
            ], ',') * ")"
        end
    end

    @testset "show(io, x)" begin
        let spec = CliOptionSpec(
            FlagOption("-a"),
            CounterOption("-b"),
            Option("-c"),
            Positional("d"),
            onerror = error,
        )
            buf = IOBuffer()
            redirect_stdout(buf) do
                show(spec)
            end
            @test String(take!(buf)) == "CliOptionSpec(" * join([
                "FlagOption(:a)",
                "CounterOption(:b)",
                "Option(:c)",
                "Positional(:d)",
            ], ',') * ")"
        end
    end
end


text:
using Neuroimaging, Test

@testset "Sensors" begin

    fname = joinpath(dirname(@__FILE__), "..", "..", "data", "test.sfp")
    s = read_sfp(fname)

    @testset "Read file" begin

        s = read_sfp(fname)
        @test length(s) == 3
    end

    @testset "Show" begin

        show(s)
        show(s[1])
    end

    @testset "Info" begin

        @test label(s) == ["Fp1", "Fpz", "Fp2"]
        @test label(s[1]) == "Fp1"
        @test labels(s) == ["Fp1", "Fpz", "Fp2"]
        @test labels(s[1]) == "Fp1"

        @test Neuroimaging.x(s) == [-27.747648u"m", -0.085967u"m", 27.676888u"m"]
        @test Neuroimaging.x(s[1]) == -27.747648u"m"

        @test Neuroimaging.y(s) == [98.803864u"m", 103.555275u"m", 99.133354u"m"]
        @test Neuroimaging.y(s[1]) == 98.803864u"m"

        @test Neuroimaging.z(s) == [34.338360u"m", 34.357265u"m", 34.457005u"m"]
        @test Neuroimaging.z(s[1]) == 34.338360u"m"
    end

    @testset "Matching" begin

        a, b = match_sensors(s, ["Fp1", "Fp2"])

    end

    @testset "Leadfield matching" begin

        L = ones(4, 3, 5)
        L[:, :, 2] = L[:, :, 2] * 2
        L[:, :, 3] = L[:, :, 3] * 3
        L[:, :, 4] = L[:, :, 4] * 4
        L[:, :, 5] = L[:, :, 5] * 5

        L, idx = match_sensors(L, ["meh", "Cz", "Fp1", "Fp2", "eh"], ["Fp2", "Cz"])

        @test size(L) == (4, 3, 2)
        @test L[:, :, 1] == ones(4, 3) * 4
        @test L[:, :, 2] == ones(4, 3) * 2
        @test idx == [4, 2]

    end

    @testset "Standard sets" begin

        @test length(EEG_64_10_20) == 64
        @test length(EEG_Vanvooren_2014) == 18
        @test length(EEG_Vanvooren_2014_Left) == 9
        @test length(EEG_Vanvooren_2014_Right) == 9

    end

    @testset "Optodes" begin

        label = "source1"
        coord = Talairach(68.3, -26.9, 8.3)
        s = Source(label, coord, Dict())
        @test isa(s, Neuroimaging.Sensor)
        @test isa(s, Neuroimaging.Optode)
        @test isa(s, Neuroimaging.Source)
        @test !isa(s, Neuroimaging.Detector)
        @test !isa(s, Neuroimaging.Electrode)
        d = Detector(label, coord, Dict())
        @test isa(d, Neuroimaging.Sensor)
        @test isa(d, Neuroimaging.Optode)
        @test isa(d, Neuroimaging.Detector)
    end

end


text:
module CUPTI

using ..APIUtils

using ..CUDA
using ..CUDA: libcupti
using ..CUDA: CUcontext, CUstream, CUdevice, CUdevice_attribute,
              CUgraph, CUgraphNode, CUgraphNodeType, CUgraphExec

# TODO: move to APIUtils
struct CUuuid
    bytes::NTuple{16,Int8}
end

using CEnum

# core library
include("libcupti_common.jl")
include("error.jl")
include("libcupti.jl")

include("wrappers.jl")

end


text:
module Brackets

export read_tree

import ...Tree

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp

struct Open end
struct Close end
struct Node; val end
struct Leaf; val end

function token!(stack, node::Node)
    push!(stack, Tree(node.val))
    return stack
end
function token!(stack, ::Close)
    node = pop!(stack)
    push!(last(stack).children, node)
    return stack
end
function token!(stack, leaf::Leaf)
    push!(last(stack).children, leaf.val)
    return stack
end

function node_string(data, from, to)
    !isvalid(data, from) && (from += 1)
    !isvalid(data, to)   && (to -= 1)
    return data[from:to]
end
const leaf_string = node_string

struct TreeReader
    tokenizer
    read
end

function TreeReader()
    close      = re"\)"
    node       = re"\([ \t\n]*[^ \)\(\t\n]*"
    leaf       = re"[^ \)\(\t\n]+"
    whitespace = re"[ \t\n]*"

    tokenizer = Automa.compile(
        node       => :(token!(stack, Node(read_node(node_string(data, ts+1, te))))),
        close      => :(token!(stack, Close())),
        leaf       => :(token!(stack, Leaf(read_leaf(leaf_string(data, ts, te))))),
        whitespace => :()
    )

    context = Automa.CodeGenContext()
    reader = @eval function (data; read_node=identity, read_leaf=identity)
        stack = Tree[Tree()]

        $(Automa.generate_init_code(context, tokenizer))
        p_end = p_eof = sizeof(data)
        while p ≤ p_eof && cs > 0
            $(Automa.generate_exec_code(context, tokenizer))
        end
        if cs != 0
            error("parsing error at position $p: code ", cs)
        end
        return last(last(stack).children)
    end

    return TreeReader(tokenizer, reader)
end

const READER = TreeReader()

"""
    read_tree(str)
    read_tree(reader, string)

Read a constituency parse tree from bracketed format.
"""
function read_tree end
read_tree(str; ks...) = READER.read(str; ks...)
read_tree(reader, str; ks...) = reader.read(str; ks...)

(reader::TreeReader)(data; ks...) = reader.read(data; ks...)

end # module


text:
using ImagineAnalyses
using Test

using ImagineInterface, Unitful, Interpolations, AxisArrays
import Unitful:s, μm

###get_cycles

##whole number of cycles
#raw array
cyc = [1.0:1.0:10.0...]
ncycs = 5
cyc_mat = repeat(cyc, 1, ncycs)
cyc_vec = reshape(cyc_mat, (ncycs * length(cyc),))
cyc_mat2 = ImagineAnalyses.get_cycles(cyc_vec, length(cyc))
@test cyc_mat == cyc_mat2 #test special case with whole number of cycles
#AxisArray
#(with sample-aligned cycles)
cyc_vec_aa = AxisArray(cyc_vec, Axis{:time}(range(1.0s, stop=ncycs*10.0s, length=ncycs*length(cyc))))
@test all(ImagineAnalyses.get_cycles(cyc_vec_aa, 10.0s) .== cyc_mat)

#(without sample-aligned cycles)
#cycle still takes 10 seconds, but we sample at 1.22 samples per second which works out evenly to 5 cycles in 61 samples
cyc_vec_itp = extrapolate(interpolate(cyc_vec, BSpline(Linear())), Flat())
#cyc_vec3 = cyc_vec_itp(range(1, stop=length(cyc_vec), length=61))
cyc_vec3 = similar(cyc_vec, (61,))
#for (i, t) in enumerate(range((1/1.22), stop=length(cyc_vec), length=61))
#cyc_vec3 = cyc_vec_itp(range(1/1.22, stop=length(cyc_vec), length=61)) #not sure why this doesn't work
for (i, t) in enumerate(range(1/1.22, stop=length(cyc_vec), length=61))
    cyc_vec3[i] = cyc_vec_itp(t)
end
cyc_vec3_aa = AxisArray(cyc_vec3, Axis{:time}(range((1.0s/1.22), stop=50.0s, length=length(cyc_vec3))))
cycs_out = ImagineAnalyses.get_cycles(cyc_vec3_aa, 10.0s)
@test size(cycs_out) == (13,5) #12.2 samples per cycle but it rounds up
same_subcyc = cycs_out[3:11,:] #due to interpolation at edges only the middle portion of the cycle vector is consistent
for i = 1:size(same_subcyc,1)
    val = same_subcyc[i,1]
    @test all(isapprox.(val,same_subcyc[i,:]))
end

##partial cycle
#raw array
cyc = [1.0:1.0:10.0...]
ncycs = 5
cyc_mat = repeat(cyc, 1, ncycs)
cyc_vec = reshape(cyc_mat, (ncycs * length(cyc),))
cyc_vec2 = push!(copy(cyc_vec), 1.0)
@test all(ImagineAnalyses.get_cycles(cyc_vec2, length(cyc)) .== cyc_mat) #should ignore the last partial cycle
#AxisArray
#(with sample-aligned cycles)
cyc_vec_aa = AxisArray(push!(cyc_vec,1.0), Axis{:time}(range(1.0s, stop=51.0s, length=51)))
@test all(ImagineAnalyses.get_cycles(cyc_vec_aa, 10.0s) .== cyc_mat)

#(without sample-aligned cycles)

#cycle still takes 10 seconds, but we sample at 1.02 samples per second which does not work out evenly
cyc_vec_itp = extrapolate(interpolate(cyc_vec, BSpline(Linear())), Flat())
#cyc_vec3 = cyc_vec_itp(range(1, stop=length(cyc_vec), length=61))
cyc_vec3 = similar(cyc_vec, (51,))
#for (i, t) in enumerate(range((1/1.22), length(cyc_vec), 61))
for (i, t) in enumerate(range(1/1.02, stop=51*(1/1.02), length=51))
    cyc_vec3[i] = cyc_vec_itp(t)
end
cyc_vec3_aa = AxisArray(cyc_vec3, Axis{:time}(range((1.0s/1.02), stop=51.0*(1s/1.02), length=length(cyc_vec3))))
cycs_out = ImagineAnalyses.get_cycles(cyc_vec3_aa, 10.0s)
@test size(cycs_out) == (11,5)
same_subcyc = cycs_out[2:10,:] #due to interpolation at edges only the middle portion of the cycle vector is consistent
for i = 1:size(same_subcyc,1)
    val = same_subcyc[i,1]
    @test all(isapprox.(val,same_subcyc[i,:]))
end

#ImagineSignal
ocpi2 = rigtemplate("ocpi-2"; sample_rate=10200*inv(Unitful.s)) #1.02 * 10k
pos = first(getpositioners(ocpi2))
append!(pos, "test", cyc_vec3*Unitful.μm)
samps_vec = get_samples(pos) #returns an AxisArray
sig_cycs = ImagineAnalyses.get_cycles(pos, (10.0s/10000))
@test all(sig_cycs.==ImagineAnalyses.get_cycles(samps_vec, 10.0s/10000))

#largest_cycle_diff
cyc = [1.0:1.0:10.0...]
ncycs = 5
cyc_mat1 = repeat(cyc, 1, ncycs)
cyc_mat2 = repeat(cyc, 1, ncycs)
cyc_mat2[:,3] .+= 1.0

cyc_vec1 = reshape(cyc_mat1, (ncycs * length(cyc),))
cyc_vec2 = reshape(cyc_mat2, (ncycs * length(cyc),))
@test all(largest_cycle_diff(cyc_vec1, 10).==0.0)
@test all(largest_cycle_diff(cyc_vec2, 10).==1.0)
@test all(largest_cycle_diff(cyc_vec2, 10; ignore_ncycs=3).==0.0) #ignore first 3

#AxisArray
cyc_vec2_aa = AxisArray(cyc_vec2, Axis{:time}(range(1.0s, stop=ncycs*10.0s, length=ncycs*length(cyc))))
@test all(largest_cycle_diff(cyc_vec2_aa, 10.0s).==1.0)
@test all(largest_cycle_diff(cyc_vec2_aa, 10.0s; ignore_ncycs=3).==0.0) #ignore first 3

#ImagineSignal
ocpi2 = rigtemplate("ocpi-2"; sample_rate=1*inv(Unitful.s))
pos = first(getpositioners(ocpi2))
append!(pos, "test", cyc_vec2*Unitful.μm)
myapprox = (x,y)->isapprox(ustrip(x),ustrip(y); rtol=sqrt(1/typemax(Int16))) #account for Int16-limited encoding of piezo signal
@test all(map(x->myapprox(x,1.0μm),  largest_cycle_diff(pos, 10.0s)))
@test all(map(x->myapprox(x,0.0μm),  largest_cycle_diff(pos, 10.0s; ignore_ncycs=3)))

samps_vec = get_samples(pos)
ocpi2 = rigtemplate("ocpi-2"; sample_rate=1*inv(Unitful.s))
pos = first(getpositioners(ocpi2))
samps_vec = vcat(fill(0.0*Unitful.μm, 111), samps_vec)
append!(pos, "test", samps_vec)
@test all(map(x->myapprox(x,1.0μm),  largest_cycle_diff(pos, 10.0s, 111))) #111 sample delay period
@test all(map(x->myapprox(x,1.0μm),  largest_cycle_diff(pos, 10.0s, 111/samprate(pos)))) #time units


text:
mutable struct TypeAst
    name ::Symbol
    args ::Vector{Union{TypeAst, String}}
    TypeAst(name::Symbol) = new(name, TypeAst[])
end

Base.push!(a::TypeAst, arg::TypeAst) = Base.push!(a.args, deepcopy(arg))
Base.push!(a::TypeAst, arg::AbstractString) = Base.push!(a.args, String(arg))


text:
@testset "STORE" begin

  mktempdir() do path
    for i in 0:4
      # test directly save
      file = joinpath(path, "builtin_type_p$i.jpkl")
      store(file, builtin_type_samples; proto=i)
      @test check_bts(file)

      # test saving to string
      sbts = stores(builtin_type_samples; proto=i)
      sfile = joinpath(path, "string_builtin_type_p$i")
      open(sfile, "w+") do f
        write(f, sbts)
      end
      @test scheck_bts(sfile)
    end
  end

end


text:
include("timedomain/timedomain.jl")
include("frequencydomain/frequencydomain.jl")


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 2069c2c0-01b4-11eb-3cc9-a7c60467a92a
# Load packages
begin
	using Pkg
	Pkg.activate(joinpath(Pkg.envdir(), "dev-env"))
	using FractalTools
	using Plots
	plotly()
	default(:label, "")
	default(:lw, 2);
end

# ╔═╡ 24b8a340-01b5-11eb-0551-318932fe4a10
@bind d html"<input type=range min=0 max=1 step=1e-1>"

# ╔═╡ 2caf47f8-047f-11eb-1f15-db8e661dc9ce
@bind niter html"<input type=range min=1 max=10 step=1>"

# ╔═╡ 69dbf05e-047f-11eb-3d6e-892350bec24f
niter, d

# ╔═╡ 402f7b62-01b5-11eb-2789-7bdb04c0d694
begin
	f(xi) = sin(2π * xi)

	# Consrucct data
	xi, dx, xf = 0, 0.1, 1
	x = collect(xi : dx : xf)
	y = sin.(2π * x)
	dvar = d * ones(length(x) - 1)

	# Interplate
	interp = interpolate(x, y, dvar, niter=niter)

	# Plots
	xd = collect(xi : 0.1dx : xf)
	plot(xd, f.(xd))
	yd = interp.(xd)
	plot!(xd, yd)
	scatter!(x, y, marker=(:circle, 3))
end

# ╔═╡ Cell order:
# ╠═2069c2c0-01b4-11eb-3cc9-a7c60467a92a
# ╠═24b8a340-01b5-11eb-0551-318932fe4a10
# ╠═2caf47f8-047f-11eb-1f15-db8e661dc9ce
# ╠═69dbf05e-047f-11eb-3d6e-892350bec24f
# ╠═402f7b62-01b5-11eb-2789-7bdb04c0d694


text:
mutable struct Sphere
    x::Float64
    y::Float64
    z::Float64
    r::Float64
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
abstract type AbstractLGSSM end

"""
    LGSSM{Ttransitions<:GaussMarkovModel, Temissions<:StructArray} <: AbstractLGSSM

A linear-Gaussian state-space model. Represented in terms of a Gauss-Markov model
`transitions` and collection of emission dynamics `emissions`.
"""
struct LGSSM{Ttransitions<:GaussMarkovModel, Temissions<:StructArray} <: AbstractLGSSM
    transitions::Ttransitions
    emissions::Temissions
end

@inline function transitions(model::LGSSM)
    return Zygote.literal_getfield(model, Val(:transitions))
end

@inline function emissions(model::LGSSM)
    return Zygote.literal_getfield(model, Val(:emissions))
end

@inline ordering(model::LGSSM) = ordering(transitions(model))

Zygote._pullback(::AContext, ::typeof(ordering), model) = ordering(model), nograd_pullback

function Base.:(==)(x::LGSSM, y::LGSSM)
    return (transitions(x) == transitions(y)) && (emissions(x) == emissions(y))
end

Base.length(model::LGSSM) = length(transitions(model))

Base.eachindex(model::LGSSM) = eachindex(transitions(model))

storage_type(model::LGSSM) = storage_type(transitions(model))

Zygote.@nograd storage_type

function is_of_storage_type(model::LGSSM, s::StorageType)
    return is_of_storage_type((transitions(model), emissions(model)), s)
end

# Get the (Gaussian) distribution over the latent process at time t=0.
x0(model::LGSSM) = x0(transitions(model))

emission_type(model::LGSSM) = eltype(emissions(model))

# Functionality for indexing into an LGSSM.

struct ElementOfLGSSM{Tordering, Ttransition, Temission}
    ordering::Tordering
    transition::Ttransition
    emission::Temission
end

@inline function ordering(x::ElementOfLGSSM)
    return Zygote.literal_getfield(x, Val(:ordering))
end

@inline function transition_dynamics(x::ElementOfLGSSM)
    return Zygote.literal_getfield(x, Val(:transition))
end

@inline function emission_dynamics(x::ElementOfLGSSM)
    return Zygote.literal_getfield(x, Val(:emission))
end

@inline function Base.getindex(model::LGSSM, n::Int)
    return ElementOfLGSSM(ordering(model), model.transitions[n], model.emissions[n])
end

# Draw a sample from the model.

function AbstractGPs.rand(rng::AbstractRNG, model::LGSSM)
    iterable = zip(ε_randn(rng, model), model)
    init = rand(rng, x0(model))
    return scan_emit(step_rand, iterable, init, eachindex(model))[1]
end

# Generate randomness used only once so that checkpointing works.
function ε_randn(rng::AbstractRNG, model::LGSSM)
    return zip(
        map(n -> ε_randn(rng, model.transitions[n]), eachindex(model.transitions)),
        map(f -> ε_randn(rng, f), model.emissions),
    )
end

step_rand(x::AbstractVector, (rng, model)) = step_rand(ordering(model), x, (rng, model))

function step_rand(::Forward, x::AbstractVector, ((ε_t, ε_e), model))
    x_next = conditional_rand(ε_t, transition_dynamics(model), x)
    y = conditional_rand(ε_e, emission_dynamics(model), x_next)
    return y, x_next
end

function step_rand(::Reverse, x::AbstractVector, ((ε_t, ε_e), model))
    y = conditional_rand(ε_e, emission_dynamics(model), x)
    x_next = conditional_rand(ε_t, transition_dynamics(model), x)
    return y, x_next
end

"""
    marginals(model::LGSSM)

Compute the complete marginals at each point in time. These are returned as a `Vector` of
length `length(model)`, each element of which is a dense `Gaussian`.
"""
function AbstractGPs.marginals(model::LGSSM)
    return scan_emit(step_marginals, model, x0(model), eachindex(model))[1]
end

step_marginals(x::Gaussian, model) = step_marginals(ordering(model), x, model)

function step_marginals(::Forward, x::Gaussian, model)
    xp = predict(x, transition_dynamics(model))
    y = predict(xp, emission_dynamics(model))
    return y, xp
end

function step_marginals(::Reverse, x::Gaussian, model)
    y = predict(x, emission_dynamics(model))
    xp = predict(x, transition_dynamics(model))
    return y, xp
end

"""
    marginals_diag(model::LGSSM)

Compute the diagonal of the marginals at each point in time. These are returned as a
`Vector` of length `length(model)`, each element of which is a diagonal `Gaussian`.
"""
function marginals_diag(model::LGSSM)
    return scan_emit(step_marginals_diag, model, x0(model), eachindex(model))[1]
end

step_marginals_diag(x::Gaussian, model) = step_marginals_diag(ordering(model), x, model)

function step_marginals_diag(::Forward, x::Gaussian, model)
    xp = predict(x, transition_dynamics(model))
    y = predict_marginals(xp, emission_dynamics(model))
    return y, xp
end

function step_marginals_diag(::Reverse, x::Gaussian, model)
    y = predict_marginals(x, emission_dynamics(model))
    xp = predict(x, transition_dynamics(model))
    return y, xp
end

# Compute the log marginal likelihood of the observations `y`.

function AbstractGPs.logpdf(
    model::LGSSM, y::AbstractVector{<:Union{AbstractVector, <:Real}},
)
    return sum(scan_emit(step_logpdf, zip(model, y), x0(model), eachindex(model))[1])
end

step_logpdf(x::Gaussian, (model, y)) = step_logpdf(ordering(model), x, (model, y))

function step_logpdf(::Forward, x::Gaussian, (model, y))
    xp = predict(x, transition_dynamics(model))
    xf, lml = posterior_and_lml(xp, emission_dynamics(model), y)
    return lml, xf
end

function step_logpdf(::Reverse, x::Gaussian, (model, y))
    xf, lml = posterior_and_lml(x, emission_dynamics(model), y)
    xp = predict(xf, transition_dynamics(model))
    return lml, xp
end

# Compute the filtering distributions.

function _filter(model::LGSSM, y::AbstractVector)
    return scan_emit(step_filter, zip(model, y), x0(model), eachindex(model))[1]
end

step_filter(x::Gaussian, (model, y)) = step_filter(ordering(model), x, (model, y))

function step_filter(::Forward, x::Gaussian, (model, y))
    xp = predict(x, transition_dynamics(model))
    xf, lml = posterior_and_lml(xp, emission_dynamics(model), y)
    return xf, xf
end

function step_filter(::Reverse, x::Gaussian, (model, y))
    xf, lml = posterior_and_lml(x, emission_dynamics(model), y)
    xp = predict(xf, transition_dynamics(model))
    return xf, xp
end

# Construct the posterior model.

function posterior(prior::LGSSM, y::AbstractVector)
    _check_inputs(prior, y)
    new_trans, xf = _a_bit_of_posterior(prior, y)
    A = zygote_friendly_map(x -> Zygote.literal_getfield(x, Val(:A)), new_trans)
    a = zygote_friendly_map(x -> Zygote.literal_getfield(x, Val(:a)), new_trans)
    Q = zygote_friendly_map(x -> Zygote.literal_getfield(x, Val(:Q)), new_trans)
    ems = Zygote.literal_getfield(prior, Val(:emissions))
    return LGSSM(GaussMarkovModel(reverse(ordering(prior)), A, a, Q, xf), ems)
end

function _check_inputs(prior, y)
    if length(prior) != length(y)
        lp = length(prior)
        ly = length(y)
        throw(error("Dimension mismatch. length(prior) is $lp, but length(y) is $ly"))
    end
end

ChainRulesCore.@non_differentiable _check_inputs(::Any, ::Any)

function _a_bit_of_posterior(prior, y)
    return scan_emit(step_posterior, zip(prior, y), x0(prior), eachindex(prior))
end
step_posterior(xf::Gaussian, (prior, y)) = step_posterior(ordering(prior), xf, (prior, y))

function step_posterior(::Forward, xf::Gaussian, (prior, y))
    t = transition_dynamics(prior)
    xp = predict(xf, t)
    new_dynamics = invert_dynamics(xf, xp, t)
    xf, _ = posterior_and_lml(xp, emission_dynamics(prior), y)
    return new_dynamics, xf
end

function step_posterior(::Reverse, x::Gaussian, (prior, y))
    xf, _ = posterior_and_lml(x, emission_dynamics(prior), y)
    t = transition_dynamics(prior)
    xp = predict(xf, t)
    return invert_dynamics(xp, xf, t), xp
end

# inlining for the benefit of type inference. Needed in at least julia-1.5.3.
@inline function invert_dynamics(xf::Gaussian, xp::Gaussian, prior::SmallOutputLGC)
    A, _, _ = get_fields(prior)
    mf, Pf = get_fields(xf)
    mp, Pp = get_fields(xp)
    U = cholesky(Symmetric(Pp + ident_eps(1e-10))).U
    Gt = U \ (U' \ (A * Pf))
    return SmallOutputLGC(_collect(Gt'), mf - Gt'mp, _compute_Pf(Pf, U * Gt))
end

_compute_Pf(Pp::AbstractMatrix, B::AbstractMatrix) = Pp - B'B

ident_eps(xf::Gaussian) = ident_eps(xf, 1e-12)

ident_eps(xf, ε) = UniformScaling(convert(eltype(xf), ε))

ident_eps(ε::Real) = UniformScaling(ε)

ident_eps(x::ColVecs, ε::Real) = UniformScaling(convert(eltype(x.X), ε))

function Zygote._pullback(::NoContext, ::typeof(ident_eps), args...)
    return ident_eps(args...), nograd_pullback
end

_collect(U::Adjoint{<:Any, <:Matrix}) = collect(U)
_collect(U::SMatrix) = U

# AD stuff. No need to understand this unless you're really plumbing the depths...

function get_adjoint_storage(
    x::LGSSM, n::Int, Δx::NamedTuple{(:ordering, :transition, :emission)},
)
    return (
        transitions = get_adjoint_storage(x.transitions, n, Δx.transition),
        emissions = get_adjoint_storage(x.emissions, n, Δx.emission)
    )
end

function _accum_at(
    Δxs::NamedTuple{(:transitions, :emissions)},
    n::Int,
    Δx::NamedTuple{(:ordering, :transition, :emission)},
)
    return (
        transitions = _accum_at(Δxs.transitions, n, Δx.transition),
        emissions = _accum_at(Δxs.emissions, n, Δx.emission),
    )
end


text:
using Statistics,Serialization,StatsPlots, DataFrames

filename_prefix = "../../figures/working_paper_2/"

path_to_data_prefix = "../../data/working_paper_2/"

mkpath(filename_prefix)

# define period of plotting for economic plots
start_week = 2 # cut of first week
end_week = 102
# set periods for covid plots
inf_start_day = 1
inf_end_day = 600 # cut of after period 400
inf_start_week = 1
inf_end_week = 86 # 58*7 = 406 and 42*7 = 294

include("../../data/working_paper_2/baseline_GER/covid_par_ini.jl") # necessary for some basic parameter values (nhh)

datapoint = fld(T,datat)

# parameter values
paras = ["beta_l5","beta_l10","beta_l30","baseline_GER"]

# scatterplot
gdploss_means = []
totcas_means = []
pubacc_means = []

for (i,p) in enumerate(paras)
    local path_to_data = string(path_to_data_prefix, p , "/")
    # store data scatterplot
    local worker_results = deserialize(open("$(path_to_data)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    local pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means,gdploss)
    push!(totcas_means,totcas)
    push!(pubacc_means,pubacc)
    if p == "beta_l5"
        global point_A_x = gdploss
        global point_A_y = totcas
    end
    if p == "baseline_GER"
            global point_D_x = gdploss
            global point_D_y = totcas
        end
end
plfull = plot(gdploss_means, totcas_means, linewidth = 2, linecolor = [:black], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means, totcas_means, linewidth = 2, seriestype = :scatter, color = [:black], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

gdploss_means = reverse(gdploss_means)
totcas_means = reverse(totcas_means)
ll = length(paras)
plarrows = quiver(gdploss_means[1:ll-1],totcas_means[1:ll-1],quiver=(gdploss_means[2:ll]-gdploss_means[1:ll-1],totcas_means[2:ll]-totcas_means[1:ll-1]), color = [:black], linewidth = 2, xlabel = "average GDP loss [%]", ylabel = "mortality [%]")

scatter!(gdploss_means, totcas_means, markersize = 3, markercolor = [:black], legend=false)

paras = ["beta_l5_alpha_l025","beta_l5_alpha_l025_alpha_o025"]

# scatterplot
gdploss_means1 = []
totcas_means1 = []
pubacc_means1 = []

for (i,p) in enumerate(paras)
    local path_to_data = string(path_to_data_prefix, p , "/")
    # store data scatterplot
    local worker_results = deserialize(open("$(path_to_data)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    local pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means1,gdploss)
    push!(totcas_means1,totcas)
    push!(pubacc_means1,pubacc)
    if p == "beta_l5_alpha_l025_alpha_o025"
        global point_E_x = gdploss
        global point_E_y = totcas
    end
end
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linecolor = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, seriestype = [:scatter], color = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, gdploss_means1[1:ll-1],totcas_means1[1:ll-1],quiver=(gdploss_means1[2:ll]-gdploss_means1[1:ll-1],totcas_means1[2:ll]-totcas_means1[1:ll-1]), linecolor = [:blue], linewidth = 2)
scatter!(gdploss_means1, totcas_means1, markersize = 3, markercolor = [:blue], legend=false)

paras = ["baseline_GER","alpha_o025","alpha_o050","alpha_o075"]

# scatterplot
gdploss_means1 = []
totcas_means1 = []
pubacc_means1 = []

for (i,p) in enumerate(paras)
    local path_to_data = string(path_to_data_prefix, p , "/")
    # store data scatterplot
    local worker_results = deserialize(open("$(path_to_data)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    local pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means1,gdploss)
    push!(totcas_means1,totcas)
    push!(pubacc_means1,pubacc)
    if p == "alpha_o050"
        global point_C_x = gdploss
        global point_C_y = totcas
    end
end
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linestyle= [:dot], linecolor = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linestyle= [:dot], seriestype = [:scatter], color = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, gdploss_means1[1:ll-1],totcas_means1[1:ll-1],quiver=(gdploss_means1[2:ll]-gdploss_means1[1:ll-1],totcas_means1[2:ll]-totcas_means1[1:ll-1]), linecolor = [:blue], linestyle= [:dot], linewidth = 2)
scatter!(gdploss_means1, totcas_means1, markersize = 3, markercolor = [:blue], linestyle= [:dot], legend=false)

paras = ["beta_l5_alpha_l025","beta_l5_alpha_l050","beta_l5_alpha_l075","beta_l5","beta_l5_alpha_l125"]

# scatterplot
gdploss_means2 = []
totcas_means2 = []
pubacc_means2 = []

for (i,p) in enumerate(paras)
    local path_to_data = string(path_to_data_prefix, p , "/")
    # store data scatterplot
    local worker_results = deserialize(open("$(path_to_data)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    local pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means2,gdploss)
    push!(totcas_means2,totcas)
    push!(pubacc_means2,pubacc)
    if p == "beta_l5_alpha_l025"
        global point_B_x = gdploss
        global point_B_y = totcas
    end
end
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, linecolor = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, seriestype= [:scatter], color = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, [gdploss_means2[ll-1]],[totcas_means2[ll-1]],quiver=([gdploss_means2[ll]-gdploss_means2[ll-1]],[totcas_means2[ll]-totcas_means2[ll-1]]), linecolor = [:red], linewidth = 2)
gdploss_means2 = reverse(gdploss_means2)
totcas_means2 = reverse(totcas_means2)
quiver!(plarrows, gdploss_means2[2:ll-1],totcas_means2[2:ll-1],quiver=(gdploss_means2[3:ll]-gdploss_means2[2:ll-1],totcas_means2[3:ll]-totcas_means2[2:ll-1]), linecolor = [:red], linewidth = 2)
scatter!(gdploss_means2, totcas_means2, markersize = 3, markercolor = [:red], legend=false)

paras = ["alpha_l025","alpha_l050","alpha_l075","baseline_GER","alpha_l125"]

# scatterplot
gdploss_means2 = []
totcas_means2 = []
pubacc_means2 = []

for (i,p) in enumerate(paras)
    local path_to_data = string(path_to_data_prefix, p , "/")
    # store data scatterplot
    local worker_results = deserialize(open("$(path_to_data)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    local pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means2,gdploss)
    push!(totcas_means2,totcas)
    push!(pubacc_means2,pubacc)
end
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, linestyle= [:dot], linecolor = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, linestyle= [:dot], seriestype= [:scatter], color = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, [gdploss_means2[ll-1]],[totcas_means2[ll-1]],quiver=([gdploss_means2[ll]-gdploss_means2[ll-1]],[totcas_means2[ll]-totcas_means2[ll-1]]), linecolor = [:red], linestyle= [:dot], linewidth = 2)
gdploss_means2 = reverse(gdploss_means2)
totcas_means2 = reverse(totcas_means2)
quiver!(plarrows, gdploss_means2[2:ll-1],totcas_means2[2:ll-1],quiver=(gdploss_means2[3:ll]-gdploss_means2[2:ll-1],totcas_means2[3:ll]-totcas_means2[2:ll-1]), linecolor = [:red], linestyle= [:dot], linewidth = 2)
scatter!(gdploss_means2, totcas_means2, markersize = 3, markercolor = [:red], linestyle= [:dot], legend=false)

### add specific points A B C
scatter!([point_A_x], [point_A_y], markersize = 4, markercolor = [:green])
scatter!([point_A_x], [point_A_y+0.0008], markersize = 0, markercolor = [:white], series_annotations = text.(["B"], :bottom))

scatter!([point_B_x], [point_B_y], markersize = 4, markercolor = [:green])
scatter!([point_B_x], [point_B_y+0.0008], markersize = 0, markercolor = [:white], series_annotations = text.(["C"], :bottom))

scatter!([point_C_x], [point_C_y], markersize = 4, markercolor = [:green])
scatter!([point_C_x-0.25], [point_C_y], markersize = 0, markercolor = [:white], series_annotations = text.(["D"], :right))

scatter!([point_D_x], [point_D_y], markersize = 4, markercolor = [:green])
scatter!([point_D_x], [point_D_y+0.001], markersize = 0, markercolor = [:white], series_annotations = text.(["A"], :bottom))

scatter!([point_E_x], [point_E_y], markersize = 4, markercolor = [:green])
scatter!([point_E_x-0.3], [point_E_y-0.002], markersize = 0, markercolor = [:white], series_annotations = text.(["E"], :left))

ylims!(plfull,0.005,0.07)
xlims!(plfull,0,9)

savefig(plarrows,"$(filename_prefix)policy_comparison.pdf")


text:
module Optimizer

#=====
params
======#
abstract type Optimizer_params end
mutable struct SGD_params <: Optimizer_params
    learning_rate::Float64
end
mutable struct Adam_params <: Optimizer_params
    learning_rate::Float64
    p1
    p2
    e
    vs #momentam
    ss #RMSProp
end

#=========
コンストラクタ
==========#
function SGD(params::Array)
    return SGD_params(params[1])
end
function Adam(params::Array)
    return Adam_params(params[1],params[2],params[3],params[4],params[5],params[6])
end

#=========
実行
==========#
function update(this::SGD_params,model)
    for i in 1:length(model.params)
        # layers[i].W = layers[i].W -this.learning_rate*layers[i].dW
        model.params[i] .-= this.learning_rate*model.grads[i]
    end
end
function update(this::Adam_params,model) #layer->model
    for i in 1:length(model.params)
        v = this.p1 .* this.vs[i] .+ (1-this.p1).*model.grads[i]
        s = this.p2 .* this.ss[i] .+ (1-this.p2).*(model.grads[i] .^2)
        model.params[i] .-= this.learning_rate*(v./sqrt.(this.e .+s))
        this.vs[i] = v
        this.ss[i] = s
    end
end

end


text:
function wald(obj::EconometricModel)
    β = coef(obj)
    V = vcov(obj)
	k = length(β) - 1
	rdf = dof_residual(obj)
	R = hasintercept(obj) ? hcat(zeros(k), I) : I
	Bread = R * β
	Meat = inv(bunchkaufman!(Hermitian(R * V * R')))
	W = (Bread' * Meat * Bread) / size(R, 1)
	if size(R, 1) > 0
		F = FDist(size(R, 1), rdf)
		p = ccdf(F, W)
	else
		F = FDist(1, 1)
		p = NaN
	end
	W, F, p
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function fem_advection(u,v,m,n,h)
    fem_advection_ = load_op_and_grad("./build/libFemAdvection","fem_advection", multiple=true)
    u,v,m_,n_,h = convert_to_tensor(Any[u,v,m,n,h], [Float64,Float64,Int64,Int64,Float64])
    ii, jj, vv = fem_advection_(u,v,m_,n_,h)
    SparseTensor(ii+1, jj+1, vv, (m+1)*(n+1), (m+1)*(n+1))
end

m = 10
n = 10
h = 0.1
u = rand(4*m*n)
v = rand(4*m*n)

# TODO: specify your input parameters
U = fem_advection(u,v,m,n,h)
sess = Session(); init(sess)
# @show run(sess, U)

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    return sum(values(fem_advection(u,x,m,n,h))^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(4*m*n))
v_ = rand(4*m*n)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
# Toplevel interface functions for recurisve broadcast expressions

"""
    slab(data::AbstractData, h::Integer)

A "pancake" view into an underlying
data layout `data` at location `h`.
"""
function slab end

# generic fallback
@inline slab(x, inds...) = x
@inline slab(tup::Tuple, inds...) = slab_args(tup, inds...)

# Recursively call slab() on broadcast arguments in a way that is statically reducible by the optimizer
# see Base.Broadcast.preprocess_args
@inline slab_args(args::Tuple, inds...) =
    (slab(args[1], inds...), slab_args(Base.tail(args), inds...)...)
@inline slab_args(args::Tuple{Any}, inds...) = (slab(args[1], inds...),)
@inline slab_args(args::Tuple{}, inds...) = ()

"""
    column(data::AbstractData, i::Integer)

A contiguous "column" view into an underlying
data layout `data` at nodal point index `i`.
"""
function column end

# generic fallback
@inline column(x, inds...) = x
@inline column(tup::Tuple, inds...) = column_args(tup, inds...)

# Recursively call column() on broadcast arguments in a way that is statically reducible by the optimizer
# see Base.Broadcast.preprocess_args
@inline column_args(args::Tuple, inds...) =
    (column(args[1], inds...), column_args(Base.tail(args), inds...)...)
@inline column_args(args::Tuple{Any}, inds...) = (column(args[1], inds...),)
@inline column_args(args::Tuple{}, inds...) = ()

function level end


text:
using Printf

"""
    printmatrix(mat [, Val(F), colnames])

Display a matrix in full while controlling the formatting of each value and
optionally showing the column names.
  * `Val(F)` - display each number in the given format `F`. The format is in the
    form of a decimal point number where the whole part indicates the total
    width and the fractional part is the number of digits printed after the
    decimal point. Default is `Val(12.7)`
  * `colnames` - a list of names to display in the first row. The names are
    displayed as given, possibly with padding to match the width given in the
    `Val` argument. If any names are longer than that, they will not be
    truncated and so the display will not be aligned properly. Sorry about that!

"""
printmatrix(mat, args...) = printmatrix(mat, Val(12.7), args...)
@generated function printmatrix(mat, ::Val{N}, cols = nothing) where N
    fmts = "% $(N)s "
    fmtn = "% $(N)f "
    return quote
        m, n = size(mat)
        if cols !== nothing
            s = ""
            for j ∈ 1:n
                s *= Printf.@sprintf($fmts, cols[j])
            end
            println(s)
        end
        for i in 1:m
            s = ""
            for j in 1:n
                s *= Printf.@sprintf($fmtn, mat[i,j])
            end
            println(s)
        end
        return nothing
    end
end
export printmatrix


text:
# Autogenerated wrapper script for WCS_jll for armv7l-linux-musleabihf
export libwcs

JLLWrappers.@generate_wrapper_header("WCS")
JLLWrappers.@declare_library_product(libwcs, "libwcs.so.7")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libwcs,
        "lib/libwcs.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, Plots, Parameters

Nx = 500
xm = 10.
x = LinRange(-xm,xm,Nx); y = x

# create some test data
psi = @. 10*exp(-.1*(x^2+y'^2)) + 15*exp(-((x-5)^2+y'^2)) # + randn()+im*randn()
heatmap(x,y,abs2.(psi))

abstract type Basis end
struct Oscillator <: Basis
    ω::Float64
    n::Int64
end
Oscillator(n::Int64) = Oscillator(1.0,n)
function (b::Oscillator)(x)
    @unpack ω,n = b
    T = x * ones(1, n) |> zero
    T[:,1] = @. exp(-(√ω * x)^2 / 2) * (ω / π)^(1 / 4)
    n >= 2 && (T[:,2] = @. sqrt(2 * ω) * x * T[:,1])
    for j ∈ index(b)[3:end]
        m = qnumbers(b)[j]
        T[:,j] = @. sqrt(2 / m) * (√ω * x) * T[:,j-1] - sqrt((m-1) / m) * T[:,j-2]
    end
    return T
end
index(b::Oscillator) = 1:b.n
spectrum(b::Oscillator) = (0:b.n-1)*b.ω
qnumbers(b::Oscillator) = 0:b.n-1

function hermite(x,n)
    @assert n >= 1
    h1 = exp(-x^2 / 2) * (1 / π)^(1 / 4)
    n == 1 && (return h1)
    h2 = sqrt(2) * x * h1
    n == 2 && (return h2)
    h3 = 0.0
    for j ∈ 3:n
        m = j - 1
        h3 = sqrt(2/m) * x * h2 - sqrt((m-1)/m) * h1
        h1 = h2; h2 = h3
    end
    n > 2 && (return h3)
end
hermite(x,n,ω) = hermite(√ω * x,n) * ω ^ ( 1 / 4 )

# test that we can create the basis, or any one of oscillator modes
Nv = 1:20
@time hermite.(x,Nv')
plot(x,hermite.(x,Nv'),legend=false,grid=false)
n = 20
b1 = Oscillator(n)
@time b1(x)
plot(x,b1(x),legend=false,grid=false)

plot(x,hermite.(x,1,[1. .5 .2 .1]))

# filtering cartesian data
function filterH(psi,x,N)
    H = Oscillator(1.0,N)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T
    psiF = zero(psi)
        for m = 1:N, n = 1:(N+1-m)
            @. psiF += H[:,m]*F̂[m,n]*H[:,n]'
        end
    return psiF
end

n = 30
psiF = filterH(psi,x,n)
heatmap(x,y,abs2.(psiF))

# coversion to polar coords
function polar(psi,x,N)
    Nx = length(x)
    H = Oscillator(1.0,N)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T
    r = LinRange(0,last(x),Nx/2 |> Int)
    θ = LinRange(0,2*pi,2*Nx)'
    psiP = zero(r*θ)
    for m = 1:N, n = 1:(N + 1 - m)
        @. psiP += hermite(r*cos(θ),n) * F̂[m,n] * hermite(r*sin(θ),m)
    end
    return psiP
end

r = LinRange(0,last(x),Nx/2 |> Int)
θ = LinRange(0,2*pi,2*Nx)'
psiP = polar(psi,x,30)

# plot in polar coordinates
heatmap(θ',r,abs2.(psiP))
dθ = θ[2]-θ[1]
dr = r[2]-r[1]
dx = x[2]-x[1]
psiPth = sum(psiP,dims=2)*dθ
plot(r,psiPth.*r)

# check total integral preserved (< 1%)
sum(psiPth.*r)*dr
sum(psiF)*dx^2

# ======== ideas for cfields
# CField:

#= something like this?
function transform(b::Oscillator)
    x,w = gausshermite(b.n)
    T = b(x)
    return x,w,T
end

xk,wk,Tk = transform(b1)
c = randn(n)
f = Tk*c
@test sum(@. wk*exp(xk^2)*abs2(f)) ≈ sum(abs2.(c))

abstract type CField end
struct YField{N} <: CField
    basis::Basis
    psi::Array{Complex{Float64},N}
end

f1 = YField(b1,randn(b1.n)+im*randn(b1.n))

=#


text:
using Test, SafeTestsets

@safetestset "Multithreading" begin include("multithreading.jl") end

@safetestset "SemObs" begin include("data_input_formats.jl") end