text:
# save real space correlations to HDF5 file
function save_χ!(
    file     :: HDF5.File,
    Λ        :: Float64,
    symmetry :: String,
    m        :: Mesh,
    χ        :: Vector{Matrix{Float64}}
    )        :: Nothing

    # save symmetry group
    if haskey(file, "symmetry") == false
        file["symmetry"] = symmetry
    end

    # save frequency mesh
    file["χ/$(Λ)/mesh"] = m.χ

    if symmetry == "su2"
        file["χ/$(Λ)/diag"] = χ[1]
    elseif symmetry == "u1-dm"
        file["χ/$(Λ)/xx"] = χ[1]
        file["χ/$(Λ)/zz"] = χ[2]
        file["χ/$(Λ)/xy"] = χ[3]
    end

    return nothing
end

"""
    read_χ_labels(
        file :: HDF5.File
        )    :: Vector{String}

Read labels of available real space correlations from HDF5 file (*_obs).
"""
function read_χ_labels(
    file :: HDF5.File
    )    :: Vector{String}

    ref    = keys(file["χ"])[1]
    labels = String[]

    for key in keys(file["χ/$(ref)"])
        if key != "mesh"
            push!(labels, key)
        end
    end

    return labels
end

"""
    read_χ(
        file    :: HDF5.File,
        Λ       :: Float64,
        label   :: String
        ;
        verbose :: Bool = true
        )       :: Tuple{Vector{Float64}, Matrix{Float64}}

Read real space correlations with name `label` and the associated frequency mesh from HDF5 file (*_obs) at cutoff Λ.
"""
function read_χ(
    file    :: HDF5.File,
    Λ       :: Float64,
    label   :: String
    ;
    verbose :: Bool = true
    )       :: Tuple{Vector{Float64}, Matrix{Float64}}

    # filter out nearest available cutoff
    list    = keys(file["χ"])
    cutoffs = parse.(Float64, list)
    index   = argmin(abs.(cutoffs .- Λ))

    if verbose
        println("Λ was adjusted to $(cutoffs[index]).")
    end

    # read frequency mesh
    m = read(file, "χ/$(cutoffs[index])/mesh")

    # read correlations with requested label
    χ = read(file, "χ/$(cutoffs[index])/" * label)

    return m, χ
end

"""
    read_χ_all(
        file    :: HDF5.File,
        Λ       :: Float64
        ;
        verbose :: Bool = true
        )       :: Tuple{Vector{Float64}, Vector{Matrix{Float64}}}

Read all available real space correlations and the associated frequency mesh from HDF5 file (*_obs) at cutoff Λ.
"""
function read_χ_all(
    file    :: HDF5.File,
    Λ       :: Float64
    ;
    verbose :: Bool = true
    )       :: Tuple{Vector{Float64}, Vector{Matrix{Float64}}}

    # filter out nearest available cutoff
    list    = keys(file["χ"])
    cutoffs = parse.(Float64, list)
    index   = argmin(abs.(cutoffs .- Λ))

    if verbose
        println("Λ was adjusted to $(cutoffs[index]).")
    end

    # read symmetry group
    symmetry = read(file, "symmetry")

    # read frequency mesh
    m = read(file, "χ/$(cutoffs[index])/mesh")

    # read correlations
    χ = Matrix{Float64}[]

    for label in read_χ_labels(file)
        push!(χ, read(file, "χ/$(cutoffs[index])/" * label))
    end

    return m, χ
end

"""
    read_χ_flow_at_site(
        file  :: HDF5.File,
        site  :: Int64,
        label :: String
        )     :: NTuple{2, Vector{Float64}}

Read flow of static real space correlations with name `label` from HDF5 file (*_obs) at irreducible site.
"""
function read_χ_flow_at_site(
    file  :: HDF5.File,
    site  :: Int64,
    label :: String
    )     :: NTuple{2, Vector{Float64}}

    # filter out a sorted list of cutoffs
    list    = keys(file["χ"])
    cutoffs = sort(parse.(Float64, list), rev = true)

    # allocate array to store values
    χ = zeros(Float64, length(cutoffs))

    # fill array with values at given site
    for i in eachindex(cutoffs)
        χ[i] = read(file, "χ/$(cutoffs[i])/" * label)[site, 1]
    end

    return cutoffs, χ
end

"""
    compute_structure_factor_flow!(
        file_in  :: HDF5.File,
        file_out :: HDF5.File,
        k        :: Matrix{Float64},
        label    :: String
        )        :: Nothing

Compute the flow of the structure factor from real space correlations with name `label` in file_in (*_obs) and save the result to file_out.
The momentum space discretization k should be formatted such that k[:, n] is the n-th momentum.
"""
function compute_structure_factor_flow!(
    file_in  :: HDF5.File,
    file_out :: HDF5.File,
    k        :: Matrix{Float64},
    label    :: String
    )        :: Nothing

    # filter out a sorted list of cutoffs
    list    = keys(file_in["χ"])
    cutoffs = sort(parse.(Float64, list), rev = true)

    # read lattice and reduced lattice
    l, r = read_lattice(file_in)

    # save momenta
    if haskey(file_out, "k") == false
        file_out["k"] = k
    end

    println("Computing structure factor flow ...")

    # compute and save structure factors
    for Λ in cutoffs
        # read correlations
        m, χ = read_χ(file_in, Λ, label, verbose = false)

        # save frequency mesh
        file_out["s/$(Λ)/mesh"] = m

        # allocate output matrix
        smat = zeros(Float64, size(k, 2), length(m))

        # compute structure factors for all Matsubara frequencies
        for w in eachindex(m)
            smat[:, w] .= compute_structure_factor(χ[:, w], k, l, r)
        end

        # save structure factors
        file_out["s/$(Λ)/" * label] = smat
    end

    println("Done.")

    return nothing
end

"""
    compute_structure_factor_flow_all!(
        file_in  :: HDF5.File,
        file_out :: HDF5.File,
        k        :: Matrix{Float64}
        )        :: Nothing

Compute the flows of the structure factors for all available real space correlations in file_in (*_obs) and save the result to file_out.
The momentum space discretization k should be formatted such that k[:, n] is the n-th momentum.
"""
function compute_structure_factor_flow_all!(
    file_in  :: HDF5.File,
    file_out :: HDF5.File,
    k        :: Matrix{Float64}
    )        :: Nothing

    # filter out a sorted list of cutoffs
    list    = keys(file_in["χ"])
    cutoffs = sort(parse.(Float64, list), rev = true)

    # read lattice and reduced lattice
    l, r = read_lattice(file_in)

    # save momenta
    if haskey(file_out, "k") == false
        file_out["k"] = k
    end

    println("Computing structure factor flows ...")

    # read available labels
    for label in read_χ_labels(file_in)
        # compute and save structure factors
        for Λ in cutoffs
            # read correlations
            m, χ = read_χ(file_in, Λ, label, verbose = false)

            # save frequency mesh
            if haskey(file_out, "s/$(Λ)/mesh") == false
                file_out["s/$(Λ)/mesh"] = m
            end

            # allocate output matrix
            smat = zeros(Float64, size(k, 2), length(m))

            # compute structure factors for all Matsubara frequencies
            for w in eachindex(m)
                smat[:, w] .= compute_structure_factor(χ[:, w], k, l, r)
            end

            # save structure factors
            file_out["s/$(Λ)/" * label] = smat
        end
    end

    println("Done.")

    return nothing
end

"""
    read_structure_factor_labels(
        file :: HDF5.File
        )    :: Vector{String}

Read labels of available structure factors from HDF5 file (*_obs).
"""
function read_structure_factor_labels(
    file :: HDF5.File
    )    :: Vector{String}

    ref    = keys(file["s"])[1]
    labels = String[]

    for key in keys(file["s/$(ref)"])
        if key != "mesh"
            push!(labels, key)
        end
    end

    return labels
end

"""
    read_structure_factor(
        file    :: HDF5.File,
        Λ       :: Float64,
        label   :: String
        ;
        verbose :: Bool = true
        )       :: Tuple{Vector{Float64}, Matrix{Float64}}

Read structure factor with name `label` and the associated frequency mesh from HDF5 file at cutoff Λ.
"""
function read_structure_factor(
    file    :: HDF5.File,
    Λ       :: Float64,
    label   :: String
    ;
    verbose :: Bool = true
    )       :: Tuple{Vector{Float64}, Matrix{Float64}}

    # filter out nearest available cutoff
    list    = keys(file["s"])
    cutoffs = parse.(Float64, list)
    index   = argmin(abs.(cutoffs .- Λ))

    if verbose
        println("Λ was adjusted to $(cutoffs[index]).")
    end

    # read frequency mesh
    m = read(file, "s/$(cutoffs[index])/mesh")

    # read structure factor with requested label
    s = read(file, "s/$(cutoffs[index])/" * label)

    return m, s
end

"""
    read_structure_factor_all(
        file    :: HDF5.File,
        Λ       :: Float64
        ;
        verbose :: Bool = true
        )       :: Tuple{Vector{Float64}, Vector{Matrix{Float64}}}

Read all available structure factors and the associated frequency mesh from HDF5 file at cutoff Λ.
"""
function read_structure_factor_all(
    file    :: HDF5.File,
    Λ       :: Float64
    ;
    verbose :: Bool = true
    )       :: Tuple{Vector{Float64}, Vector{Matrix{Float64}}}

    # filter out nearest available cutoff
    list    = keys(file["s"])
    cutoffs = parse.(Float64, list)
    index   = argmin(abs.(cutoffs .- Λ))

    if verbose
        println("Λ was adjusted to $(cutoffs[index]).")
    end

    # read frequency mesh
    m = read(file, "s/$(cutoffs[index])/mesh")

    # read structure factors
    s = Matrix{Float64}[]

    for label in read_structure_factor_labels(file)
        push!(s, read(file, "s/$(cutoffs[index])/" * label))
    end

    return m, s
end

"""
    read_structure_factor_flow_at_momentum(
        file    :: HDF5.File,
        p       :: Vector{Float64},
        label   :: String
        ;
        verbose :: Bool = true
        )       :: NTuple{2, Vector{Float64}}

Read flow of static structure factor with name `label` from HDF5 file at momentum p.
"""
function read_structure_factor_flow_at_momentum(
    file    :: HDF5.File,
    p       :: Vector{Float64},
    label   :: String
    ;
    verbose :: Bool = true
    )       :: NTuple{2, Vector{Float64}}

    # filter out a sorted list of cutoffs
    list    = keys(file["s"])
    cutoffs = sort(parse.(Float64, list), rev = true)

    # locate closest momentum
    k     = read(file, "k")
    dists = Float64[norm(k[:, i] .- p) for i in 1 : size(k, 2)]
    index = argmin(dists)

    if verbose
        println("Momentum was adjusted to k = $(k[:, index]).")
    end

    # allocate array to store values
    s = zeros(Float64, length(cutoffs))

    # fill array with values at given momentum
    for i in eachindex(cutoffs)
        s[i] = read(file, "s/$(cutoffs[i])/" * label)[index, 1]
    end

    return cutoffs, s
end

"""
    read_reference_momentum(
        file  :: HDF5.File,
        Λ     :: Float64,
        label :: String
        )     :: Vector{Float64}

Read the momentum with maximum static structure factor value with name `label` at cutoff Λ from HDF5 file.
"""
function read_reference_momentum(
    file  :: HDF5.File,
    Λ     :: Float64,
    label :: String
    )     :: Vector{Float64}

    # read struture factor
    m, s = read_structure_factor(file, Λ, label)

    # determine momentum with maximum amplitude
    k = read(file, "k")
    p = k[:, argmax(s[:, 1])]

    return p
end


text:
struct precipMesh
    sMesh::eFEM.ScalarMesh
    fMesh::eFEM.FluidMesh
end

struct precipSolution
    ψA::eFEM.ScalarSolution
    ψB::eFEM.ScalarSolution
    ψC::eFEM.ScalarSolution
    θs::eFEM.ScalarSolution
    θf::eFEM.ScalarSolution
    fluid::eFEM.FluidSolution
end

struct precipProblem
    ψA::eFEM.Problem
    ψB::eFEM.Problem
    ψC::eFEM.Problem
    fluid::eFEM.Problem
end

@with_kw struct precipParam
    # --------------------
    # Physical Parameters
    # --------------------
    μ0::Float64  = 1.0  # viscosity
    μ1::Float64  = 1.0  # brinkman viscosity
    α0::Float64  = 0.0  # initial value of thetam/thetaf
    G::Float64   = 1e0  # pressure inflow
    ξϵ::Float64  = 1e0  # regularization of friction coefficient

    PeA::Float64 = 30    # Peclet number for chemical A
    PeB::Float64 = 30    # Peclet number for chemical B
    PeC::Float64 = 1     # Peclet number for chemical C

    ρf::Float64 = 1.0     # mass density of fluid
    ρs::Float64 = 4.1     # mass density of solid

    # computed parameters
    α::Float64 = 1
    h::Float64 = 1

    # ---------------------
    # Chemistry Parameters aA + bB -> cC
    # ---------------------
    a::Int = 1
    b::Int = 1
    c::Int = 1
    Ma::Float64 = 58.693            # molar mass of A
    Mb::Float64 = 17.008            # molar mass of A
    Mc::Float64 = 0.0               # molar mass of C

    # ----------------------------
    # Time-integration parameters
    # ----------------------------
    Δt::Float64        = 5e-2
    Tfinal::Float64    = 80.0
    integrator::Symbol = :cranknicholson   # :cranknicholson, :backwardeuler

    # ---------------
    # vtk Parameters
    # ---------------
    ψAname::String = "aqueous_A"
    ψBname::String = "aqueous_B"
    ψCname::String = "aqueous_C"
    θsname::String = "membrane_vf"
    θfname::String = "fluid_vf"
    printSkipIndex::Int = 10
    printBool::Bool     = true
    SaveFolder::String  = string(homedir(),"/data/6-7-2020-y-dom")
    SaveName::String    = "y-dom"

    # ----------------
    # Mesh Parameters
    # ----------------
    MeshFolder::String="/home/peastham/Dropbox/projects/project-precipitation---numerics/PrecipSim.jl/meshes/y-dom"
    MeshFile::String = "y-dom6.msh"
    MeshTransformFunc::Function = none
    MeshTransformParam::Float64 = 1.0

    # -------------
    # run ID (1-7)
    # -------------
    runID::Int = 1
  end

ext: jl


text:
# https://github.com/JuliaGraphs/JuliaGraphsTutorials
using Test
using LightGraphs
using GraphPlot
using Compose

## Getting Started

# https://juliagraphs.org/LightGraphs.jl/latest/#Basic-library-examples
function learn_basic_library_examples(io::IO)
    # A path that connect 1 vertex to the next
    g = path_graph(6)
    @test 6 == nv(g)
    @test 5 == ne(g)

    add_edge!(g, 1, 6) # connect first and last
    @test 6 == nv(g)
    @test 6 == ne(g)
end

## Making and Modifying Graphs

# https://juliagraphs.org/LightGraphs.jl/latest/generators/#Graph-Generators
function learn_graph_generators(io::IO)
    # 1 - 3 - 5
    # |   |   |
    # 2 - 4 - 6
    g = grid([2,3], periodic=false)
    @test 6 == nv(g)
    @test 7 == ne(g)

    plt = gplot(g, nodelabel=1:nv(g), edgelabel=1:ne(g))
    display(plt)

    # Numbering follow the 1st dimension first then 2nd dimension
    g = grid([10,5], periodic=false)
    plt = gplot(g, nodelabel=1:nv(g))
    display(plt)
end
function learn_modifying_graphs()
    # 1 - 3 - 5
    # |   |   |
    # 2 - 4 - 6

    # Remove an edge
    g = grid([2,3], periodic=false)
    @test 7 == ne(g)
    @test has_edge(g, 5,6)

    @test rem_edge!(g, 5,6)
    @test 6 == ne(g)
    @test !has_edge(g, 5,6)

    # Add existing edge does not change the graph
    g = grid([2,3], periodic=false)
    @test 7 == ne(g)
    @test !add_edge!(g, 5,6)
    @test 7 == ne(g)
end

## Reading/Writing Graphs

## Operators

## Plotting Graphs

## Path and Traversal
function learn_bfs_tree()
    # 1 - 3 - 5
    # |   |   |
    # 2 - 4 - 6
    g = grid([2,3], periodic=false)
    #display(gplot(g, nodelabel=1:nv(g), edgelabel=1:ne(g)))

    bfs = bfs_tree(g, 1)
    # 6 <- 4 <- 2 <- 1 -> 3 -> 5
    #display(gplot(bfs, nodelabel=1:nv(bfs), edgelabel=1:ne(bfs)))
    draw(SVG("bfs_tree.svg", 16cm, 16cm), gplot(bfs, nodelabel=1:nv(bfs), edgelabel=1:ne(bfs)))

    @test 6 == nv(bfs)
    @test 5 == ne(bfs)

    edge = first(edges(bfs))
    @test edge isa LightGraphs.SimpleGraphs.SimpleEdge
    @test 1 == src(edge)
    @test 2 == dst(edge)

    # 1 - 3   5
    # |   |   |
    # 2 - 4   6
    g = grid([2,3], periodic=false)

    rem_edge!(g, 3,5)
    rem_edge!(g, 4,6)

    @test 6 == nv(bfs)
    @test 5 == ne(bfs)

    bfs = bfs_tree(g, 2)
    @show collect(vertices(bfs))
    @show collect(edges(bfs))
    @test 6 == nv(bfs)
    @test 3 == ne(bfs)
end

function learn_dijkstra_shortest_paths(dp)
    g = grid([2,2], periodic=false)

    # weight matrix
    w = [
        0 2 2 0;    # distance of v1 to v1-v4
        2 0 0 100;  # distance of v2 to v1-v4
        2 0 0 50;
        0 100 50 0;
    ]

    edgelabel = Vector{String}(undef, ne(g))
    for (i,e) in enumerate(edges(g))
        #@show i, e
        edgelabel[i] = string(w[e.src, e.dst])
    end
    ds = dijkstra_shortest_paths(g, 1, w)
    @show ds
    @show ds.dists

    plt = gplot(g, nodelabel=1:nv(g), edgelabel=edgelabel)
    #display(plt)

    return ds
end

function learn_dijkstra_shortest_paths_2(dp)
    g = grid([2,2], periodic=false)

    # Assign permeability to each grid block (vertex)
    k = Float64.([100 100 50 1]) .* 0.01

    # Calculate transmissibility
    R = Matrix{Float64}(undef, nv(g), nv(g)) # resistance
    for edge in edges(g)
        src = edge.src
        dst = edge.dst
        r = 0.5*(1/k[src] + 1/k[dst])
        R[src,dst] = r
        R[dst,src] = r # symmetric
    end

    edgelabel = Vector{String}(undef, ne(g))
    for (i,e) in enumerate(edges(g))
        edgelabel[i] = string(R[e.src, e.dst])
    end

    ds = dijkstra_shortest_paths(g, 1, R)
    @info "Distance from Vertex 1 to: $(ds.dists)"

    if dp
        plt = gplot(g, nodelabel=1:nv(g), edgelabel=edgelabel)
        display(plt)
    end
end

## Distance

function learn()
    g = path_graph(0)

    @show g
    @show nv(g)
    @show ne(g)

    add_vertex!(g)
    add_vertex!(g)

    for v in vertices(g)
        @show v, typeof(v)
    end
end

# https://nbviewer.jupyter.org/github/JuliaGraphs/JuliaGraphsTutorials/blob/master/Basics.ipynb
function learn_Basics()
    # Graph -> SimpleGraph
    G = Graph(3) # graph with 3 vertices

    # make a triangle
    add_edge!(G, 1, 2)
    add_edge!(G, 1, 3)
    add_edge!(G, 2, 3)

    plt = gplot(G, nodelabel=1:3)
    display(plt)

    # Adjacency matrix
    A = [
        0 1 1 # node 1 connects to 2, 3
        1 0 1
        1 1 0
    ]

    G2 = Graph(A)

    @test G == G2
end

# https://nbviewer.jupyter.org/github/JuliaGraphs/JuliaGraphsTutorials/blob/master/Basics.ipynb#Graph-properties
function learn_graph_properties(io::IO)
    G = smallgraph("house")
    nvertices = nv(G)
    nedges = ne(G)

    # Get vertices
    for v in vertices(G)
        println(io, "vertex $v")
    end

    # Get edges
    for e in edges(G)
        u = src(e)
        v = dst(e)
        println(io, "edge $u - $v")
    end

    plt = gplot(G, nodelabel=1:nvertices, edgelabel=1:nedges)
    display(plt)
end

# https://nbviewer.jupyter.org/github/JuliaGraphs/JuliaGraphsTutorials/blob/master/Basics.ipynb#Basic-operations
function learn_basic_operations(io::IO)
    G = Graph()
    add_vertices!(G, 3)
    add_vertices!(G, 2)

    add_edge!(G, 1, 2)
    add_edge!(G, 1, 3)
    add_edge!(G, 2, 4)
    add_edge!(G, 3, 4)
    add_edge!(G, 3, 5)
    add_edge!(G, 4, 5)

    plt = gplot(G, nodelabel=1:nv(G), edgelabel=1:ne(G))
    display(plt)

    # add vertex
    add_vertex!(G)
    add_edge!(G, 5, 6)

    plt = gplot(G, nodelabel=1:nv(G), edgelabel=1:ne(G))
    display(plt)

    # remove vertex
    rem_vertex!(G, 1)
    plt = gplot(G, nodelabel=1:nv(G), edgelabel=1:ne(G))
    display(plt)
end

function learn_set_interface(io::IO)
    skeleton = Graph(11)
    add_edge!(skeleton, 1, 2)
    add_edge!(skeleton, 2, 3)
    add_edge!(skeleton, 3, 4)
    add_edge!(skeleton, 4, 5)
    add_edge!(skeleton, 3, 6)
    add_edge!(skeleton, 3, 7)
    add_edge!(skeleton, 3, 8)
    add_edge!(skeleton, 3, 9)
    add_edge!(skeleton, 9, 10)
    add_edge!(skeleton, 9, 11)

    plt = gplot(skeleton)
    display(plt)

    □ = Graph(4)
    add_edge!(□, 1, 2)
    add_edge!(□, 1, 3)
    add_edge!(□, 2, 4)
    add_edge!(□, 3, 4)

    plt = gplot(cartesian_product(□, skeleton))
    display(plt)
end

# https://nbviewer.jupyter.org/github/JuliaGraphs/JuliaGraphsTutorials/blob/master/DAG-Julia-Pkgs.ipynb
function learn_dag(io)
    G = smallgraph("house")

    plt = gplot(G, nodelabel=1:nv(G), edgelabel=1:ne(G))
    display(plt)

    @show gdistances(G, 5)
end

function learn_MetaGraphs(io)
end

#default_logger = global_logger()
#global_logger(ConsoleLogger(stdout, Logging.Info))

io = stdout
#io = devnull

# Display Plot
dp = true
dp = false

@testset "Base" begin
    #learn_basic_library_examples(io)
    #learn_graph_generators(io)
    learn_modifying_graphs()
    ds = learn_dijkstra_shortest_paths(dp)
    learn_bfs_tree()

    #learn_dijkstra_shortest_paths_2(dp)
    #learn()
    #pl = learn_Basics()
    #learn_graph_properties(io)
    #learn_basic_operations(io)
    #learn_set_interface(io)
    #learn_dag(io)
end
#global_logger(default_logger) # Restore global logger
nothing


text:
module LEAPfunctions
using DelimitedFiles, PyCall, DataFrames, CSV

export visible, outputtoleap, calculateleap, energyinvestment

include("./IOlib.jl")
using .IOlib

"""
    visible(state::Bool)

Hide or show LEAP by setting visibility.
"""
function visible(state::Bool)
	LEAP = connecttoleap()
	LEAP.Visible = state
	disconnectfromleap(LEAP)
end

"""
    outputtoleap(file::String, indices::Array)

First obtain LEAP branch info from the YAML config file and then send Macro model results to LEAP.
"""
function outputtoleap(file::String, indices::Array, run::Int64)
    # LEAP parameters from config file
    if run == 0
        params = IOlib.parse_input_file(file, force = true)
    else
        params = IOlib.parse_input_file(file)
    end
    base_year = params["years"]["start"]
    final_year = params["years"]["end"]

    # connects program to LEAP
    LEAP = connecttoleap()

    # Set ActiveView
    LEAP.ActiveView = "Analysis"

	branch_data = Dict(:branch => String[], :variable => String[], :last_historical_year => Int64[], :col => Int64[])
	col = 0
    if haskey(params, "GDP-branch")
        col += 1
        append!(branch_data[:branch], [params["GDP-branch"]["branch"]])
        append!(branch_data[:variable], [params["GDP-branch"]["variable"]])
        append!(branch_data[:last_historical_year], [params["GDP-branch"]["last_historical_year"]])
        append!(branch_data[:col], [col])
    end
    if haskey(params, "Employment-branch")
        col += 1
        append!(branch_data[:branch], [params["Employment-branch"]["branch"]])
        append!(branch_data[:variable], [params["Employment-branch"]["variable"]])
        append!(branch_data[:last_historical_year], [params["Employment-branch"]["last_historical_year"]])
        append!(branch_data[:col], [col])
    end

	for leap_sector in params["LEAP-sectors"]
		col += 1
		for branch in leap_sector["branches"]
			append!(branch_data[:branch], [branch["branch"]])
			append!(branch_data[:variable], [branch["variable"]])
			append!(branch_data[:last_historical_year], [branch["last_historical_year"]])
			append!(branch_data[:col], [col])
		end
	end

	branch_df = DataFrame(branch_data)

    # send results to LEAP
    ndxrows = final_year - base_year + 1
    try
        for i = 1:size(branch_df, 1) # loops through each branch path
            branch = branch_df[i,:branch]
            variable = branch_df[i,:variable]
            lasthistoricalyear = branch_df[i,:last_historical_year]
            col = branch_df[i,:col]
            start_ndx = (1+(col*ndxrows))
            end_ndx = (col+1)*ndxrows

            if lasthistoricalyear > base_year
                newexpression = interp_expression(base_year, indices[start_ndx:end_ndx], lasthistoricalyear=lasthistoricalyear)
            else
                newexpression = interp_expression(base_year, indices[start_ndx:end_ndx])
            end
            setbranchvar_expression(LEAP, branch, variable, newexpression, scenario=params["LEAP-info"]["input_scenario"])

        end
        LEAP.SaveArea()
    finally
	    disconnectfromleap(LEAP)
    end
end

"""
    connecttoleap()

Connect to the currently running instance of LEAP, if one exists; otherwise starts an instance of LEAP.

Return a `PyObject` corresponding to the instance.
If LEAP cannot be started, return `missing`
"""
function connecttoleap()
	try
		LEAPPyObj = pyimport("win32com.client").Dispatch("Leap.LEAPApplication")
        max_loops = 5
        while !LEAPPyObj.ProgramStarted & max_loops > 0
            sleep(5)
            max_loops -= 1
        end
        if !LEAPPyObj.ProgramStarted
            error("LEAP is not responding.")
            return missing
        else
            return LEAPPyObj
        end
	catch
        error("Cannot connect to LEAP. Is it installed?")
		return missing
	end
end  # connecttoleap

"""
    disconnectfromleap(LEAPPyObj)

Wrapper for PyCall's pydecref(obj), after saving
"""
function disconnectfromleap(LEAPPyObj)
	pydecref(LEAPPyObj)
end

"""
    interp_expression(base_year::Int64, newdata::Array; lasthistoricalyear::Int64=0)

Create LEAP Interp expression from an array of values.
"""
function interp_expression(base_year::Int64, newdata::Array; lasthistoricalyear::Int64=0)
    # Creates start of expression. Includes historical data if available
    if lasthistoricalyear > 0
        newexpression = string("If(year <= ", lasthistoricalyear, ", ScenarioValue(Current Accounts), Value(", base_year,") * Interp(")
        diff = lasthistoricalyear - base_year + 2
        year = lasthistoricalyear + 1
    else
        newexpression = string("(Value(", base_year,") * Interp(")
        diff = 2
        year = base_year + 1
    end

    # Incorporates IO model data into the rest of the expression
    for i = diff:size(newdata,1)
        if isnan(newdata[i]) == false
            newexpression = string(newexpression, year, ", ", newdata[i], ", ")
        end
        if i == size(newdata,1)
            newexpression = newexpression[1:(lastindex(newexpression)-2)]
            newexpression = string(newexpression, "))")
        end
        year = year + 1
    end
    return newexpression
end

"""
    setbranchvar_expression(leapapplication::PyObject, branch::String, variable::String, newexpression::String; region::String = "", scenario::String = "")

Set a LEAP branch-variable expression.

The region and scenario arguments can be omitted by leaving them as empty strings.
Note that performance is best if neither region nor scenario is specified.
"""
function setbranchvar_expression(leapapplication::PyObject, branch::String, variable::String, newexpression::String; region::String = "", scenario::String = "")
    # Set ActiveRegion and ActiveScenario as Julia doesn't allow a function call (ExpressionRS) to be set to a value
    if region != ""
        leapapplication.ActiveRegion = region
    end

    if scenario != ""
        leapapplication.ActiveScenario = scenario
    end

    # Set expression
    leapapplication.Branch(branch).Variable(variable).Expression = newexpression

    # Refresh LEAP display
    leapapplication.Refresh()
end  # setbranchvarexpression

"""
    calculateleap(scen_name::String)

Calculate the LEAP model, returning results for the specified scenario.
"""
function calculateleap(scen_name::String)
    # connects program to LEAP
    LEAP = connecttoleap()
    try
        LEAP.Scenario(scen_name).ResultsShown = true
        LEAP.Calculate()
        LEAP.SaveArea()
    finally
	    disconnectfromleap(LEAP)
    end
end

"""
    energyinvestment(file::String, run::Int64)

Obtain energy investment data from the LEAP model.
"""
function energyinvestment(file::String, run::Int64)
    # LEAP parameters from config file
    if run == 0
        params = IOlib.parse_input_file(file, force = true)
    else
        params = IOlib.parse_input_file(file)
    end
    base_year = params["years"]["start"]
    final_year = params["years"]["end"]

    # connects program to LEAP
    LEAP = connecttoleap()

    # Set ActiveView and ActiveScenario
    LEAP.ActiveView = "Results"
    LEAP.ActiveScenario = params["LEAP-info"]["result_scenario"]

    nrows = (final_year - base_year) + 1
    I_en = Array{Float64}(undef, nrows)
    I_en_temp = Array{Float64}(undef, nrows)

    n_energy = 0
    try
        for b in LEAP.Branches
            if b.BranchType == 2 && b.Level == 2 && b.VariableExists("Investment Costs")
                for y = base_year:final_year
                    I_en_temp[(y-base_year+1)] = b.Variable("Investment Costs").Value(y, params["LEAP-info"]["inv_costs_unit"]) / params["LEAP-info"]["inv_costs_scale"]
                end
                I_en = hcat(I_en, I_en_temp)
                n_energy += 1
            end
        end
    finally
    	disconnectfromleap(LEAP)
    end

    if n_energy > 0
        I_en = sum(I_en[:,2:size(I_en, 2)], dims=2)
    else
        I_en .= 0.0
    end
    writedlm(joinpath(params["results_path"],string("I_en_",run,".csv")), I_en, ',')

    return I_en
end

end

max_stars_repo_name: sei-international/LEAPMacro

max_issues_repo_name: sei-international/LEAPMacro

max_forks_repo_name: sei-international/LEAPMacro


text:
#
# tlibs2 julia module test
# @author Tobias Weber <tobias.weber@tum.de>, <tweber@ill.fr>
# @date 2017 -- 2018
# @license GPLv3, see <<LICENSE>> file
# @desc Forked on 7-Nov-2018 from my privately and TUM-PhD-developed <<tlibs>> project (https://github.com/t-weber/tlibs).
#
# ----------------------------------------------------------------------------
# tlibs
# Copyright (C) 2017-2021  Tobias WEBER (Institut Laue-Langevin (ILL),
#                          Grenoble, France).
# Copyright (C) 2015-2017  Tobias WEBER (Technische Universitaet Muenchen
#                          (TUM), Garching, Germany).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 3 of the License.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
# ----------------------------------------------------------------------------
#

include("tl2.jl")

# load data file
#a = tl2.loadinstr("/Users/tw/Measurements/mnsi/mira-18/data/13511_00035079.dat")
a = tl2.loadinstr("/home/tw/Measurements/mnsi/mira_18/data/13511_00035079.dat")
cols = a[1]
data = a[2]
#println(cols)
#println(typeof(data))

idxx = findfirst(isequal("h"), cols)
idxCtr = findfirst(isequal("ctr1"), cols)

xs = data[:,idxx]
cts = data[:,idxCtr]
cts_err = tl2.poisson_err.(cts)
#println(typeof(xs))
#println(typeof(cts))
#println(typeof(cts_err))

# print data
println(lpad(cols[idxx], 12), " ", lpad(cols[idxCtr], 12))
for (x, ct, err) in zip(xs, cts, cts_err)
	println(lpad(x, 12), " ", lpad(ct, 12), " ", lpad(err, 12))
end
println()

# fit gaussian
fitresult = tl2.fit(tl2.gauss_model_amp, xs, cts, cts_err, fixed = [ "offs" ],
	values = Dict("x0" => 1., "sigma" => 0.003, "amp" => 20000., "offs" => 5.),
	errors = Dict("x0" => 0.05, "sigma" => 0.002, "amp" => 2000., "offs" => 2.5))
println(fitresult)
println()

data_fine = tl2.eval_func(fitresult, minimum(xs), maximum(xs))

# write results
f = open("tst_pts.dat", "w")
for (x,y,yerr) in zip(xs, cts, cts_err)
	println(f, rpad(x, 20), " ", rpad(y, 20), " ", rpad(yerr, 20))
end
close(f)

f = open("tst_fit.dat", "w")
for (x,y) in zip(data_fine[1], data_fine[2])
	println(f, rpad(x, 20), " ", rpad(y, 20))
end
close(f)

# plot results
run(`gnuplot -p -e "plot \"tst_pts.dat\" w yerrorbars pt 7 lc rgb \"#ff0000\",
	\"tst_fit.dat\" with lines lw 2 lc rgb \"#0000ff\""`)


text:
using SimJulia, Distributions, BenchmarkTools

function exp_source(sim::Simulation, lambd::Float64, server::Resource, mu::Float64)
  while true
    dt = rand(Exponential(1/lambd))
    yield(timeout(sim, dt))
    @oldprocess customer(sim, server, mu)
  end
end

function customer(sim::Simulation, server::Resource, mu::Float64)
  yield(request(server))
  dt = rand(Exponential(1/mu))
  yield(timeout(sim, dt))
  yield(release(server))
end

function customer2(sim::Simulation, server::Resource, mu::Float64)
  request(server) do req
    yield(req)
    dt = rand(Exponential(1/mu))
    yield(timeout(sim, dt))
  end
end

function test_mm1(n::Float64)
  sim = Simulation()
  server = Resource(sim)
  @oldprocess exp_source(sim, 1.0, server, 1.1)
  run(sim, n)
end

test_mm1(100.0)
@btime test_mm1(100.0)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    Strata(record)

Stratigraphic model from geological `record`.
"""
struct Strata
  horizons::Vector{<:Matrix}
end

function Strata(record::Record{LandState}, stack=:erosional)
  # land maps for all states in record
  horizons = [state.land for state in record]

  if stack == :erosional
    # erode land maps backward in time
    for t=length(horizons):-1:2
      Lt = horizons[t]
      Lp = horizons[t-1]
      erosion = Lp .> Lt
      Lp[erosion] = Lt[erosion]
    end
  end

  if stack == :depositional
    # deposit sediments forward in time
    for t=1:length(horizons)-1
      Lt = horizons[t]
      Lf = horizons[t+1]
      nodeposit = Lf .≤ Lt
      Lf[nodeposit] = Lt[nodeposit]
    end
  end

  Strata(horizons)
end

function voxelize(strata::Strata, nz::Int; fillbase=NaN, filltop=NaN)
  # retrieve horizons
  horizons = strata.horizons

  # initial land map
  init   = copy(horizons[1])
  init .-= minimum(init[.!isnan.(init)])
  nx, ny = size(init)

  # estimate bulk sediment
  thickness = diff(horizons)

  # handle NaNs
  init[isnan.(init)] .= 0.
  for t=1:length(thickness)
    thick = thickness[t]
    thick[isnan.(thick)] .= 0.
  end

  # estimate maximum z coordinate
  zmax = maximum(init + sum(thickness))

  # current elevation
  elevation = floor.(Int, (init/zmax)*nz)

  # voxel model
  model = fill(NaN, nx, ny, nz)

  # fill model base
  if !isnan(fillbase)
    for j=1:ny, i=1:nx
      for k=1:elevation[i,j]
        model[i,j,k] = fillbase
      end
    end
  end

  # add layers to model
  for t=1:length(thickness)
    thick = thickness[t]
    sediments = floor.(Int, (thick/zmax)*nz)
    for j=1:ny, i=1:nx
      for k=1:sediments[i,j]
        model[i,j,elevation[i,j]+k] = t
      end
    end
    elevation .+= sediments
  end

  # fill model top
  if !isnan(filltop)
    for j=1:ny, i=1:nx
      for k=elevation[i,j]+1:nz
        model[i,j,k] = filltop
      end
    end
  end

  model
end

@recipe function f(strata::Strata)
  aspect_ratio --> :equal
  seriestype --> :surface
  colorbar --> false
  for horizon in strata.horizons
    @series begin
      horizon
    end
  end
end


text:
function apply_boundary_conditions_on_φ_axis!(  rbpot::Array{T, 4}, iz::Int, ir::Int, rbi::Int, nrbi::Int, ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T},
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[iz,   1, ir, rbi] = rbpot[ iz, end - 1, ir, rbi] # cycling boundary
    rbpot[iz, end, ir, rbi] = rbpot[ iz,       2, ir, rbi] # cycling boundary
    nothing
end

function apply_boundary_conditions_on_φ_axis!(  rbpot::Array{T, 4}, iz::Int, ir::Int, rbi::Int, nrbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{false})::Nothing where {T}
    rbpot[iz,   1, ir, rbi] = rbpot[ iz,       3, ir, rbi] # cycling boundary
    rbpot[iz, end, ir, rbi] = rbpot[ iz, end - 2, ir, rbi] # cycling boundary
    nothing
end
function apply_boundary_conditions_on_φ_axis!(  rbpot::Array{T, 4}, iz::Int, ir::Int, rbi::Int, nrbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{true})::Nothing where {T}
    rbpot[iz,   1, ir, nrbi] = rbpot[ iz, 2, ir, rbi] # cycling boundary
    rbpot[iz, end, ir, nrbi] = rbpot[ iz, 2, ir, rbi] # cycling boundary
    nothing
end
function apply_boundary_conditions_on_r_axis!(  rbpot::Array{T, 4}, iz::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :r0, :infinite}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[iz, iφ, end, rbi] = grid_boundary_factors[2] * rbpot[iz, iφ, end - 2, rbi] # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_r_axis!(  rbpot::Array{T, 4}, iz::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :r0, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[iz, iφ, end, rbi] = rbpot[iz, iφ, end - 2, rbi] # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_r_axis!(  rbpot::Array{T, 4}, iz::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :r0, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    # rbpot[iz, iφ, end, rbi] = rbpot[iz, iφ, end - 2, rbi] # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_cyl_z_axis!(  rbpot::Array{T, 4}, ir::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :infinite}, int::Interval{:closed, :closed, T},
                                                    grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iφ, ir, rbi] = grid_boundary_factors[1] * rbpot[ 2, iφ, ir, rbi]  # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, iφ, ir, rbi] = grid_boundary_factors[2] * rbpot[ end - 1, iφ, ir, rbi]    # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_cyl_z_axis!(  rbpot::Array{T, 4}, ir::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :fixed}, int::Interval{:closed, :closed, T},
                                                    grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iφ, ir, rbi] = grid_boundary_factors[1] * rbpot[ 2, iφ, ir, rbi]  # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_cyl_z_axis!(  rbpot::Array{T, 4}, ir::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :infinite}, int::Interval{:closed, :closed, T},
                                                    grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[end, iφ, ir, rbi] = grid_boundary_factors[2] * rbpot[ end - 1, iφ, ir, rbi]    # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_cyl_z_axis!(  rbpot::Array{T, 4}, ir::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T},
                                                    grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iφ, ir, rbi] = rbpot[ 2, iφ, ir, rbi]  # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, iφ, ir, rbi] = rbpot[ end - 1, iφ, ir, rbi]    # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_cyl_z_axis!(  rbpot::Array{T, 4}, ir::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :reflecting}, int::Interval{:closed, :closed, T},
                                                    grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[end, iφ, ir, rbi] = rbpot[ end - 1, iφ, ir, rbi]    # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_cyl_z_axis!(  rbpot::Array{T, 4}, ir::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :fixed}, int::Interval{:closed, :closed, T},
                                                    grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iφ, ir, rbi] = rbpot[ 2, iφ, ir, rbi]  # infinity boundaries in z (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_cyl_z_axis!(  rbpot::Array{T, 4}, ir::Int, iφ::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :fixed}, int::Interval{:closed, :closed, T},
                                                    grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    nothing
end

function apply_boundary_conditions!(fssrb::PotentialSimulationSetupRB{T, 3, 4, :cylindrical}, update_even_points::Val{even_points}, only2d::Val{only_2d}) where {T, even_points, only_2d}
    rbi::Int = even_points ? rb_even::Int : rb_odd::Int
    nrbi::Int = even_points ? rb_odd::Int : rb_even::Int
    ax1::typeof(fssrb.grid.axes[1]) = fssrb.grid.axes[1]
    ax2::typeof(fssrb.grid.axes[2]) = fssrb.grid.axes[2]
    ax3::typeof(fssrb.grid.axes[3]) = fssrb.grid.axes[3]
    int1::typeof(ax1.interval) = ax1.interval
    int2::typeof(ax2.interval) = ax2.interval
    int3::typeof(ax3.interval) = ax3.interval
    if only_2d
        iφ::Int = 2
        @inbounds for iz in axes(fssrb.potential, 1)
            # for ir in axes(fssrb.potential, 3)
            #     apply_boundary_conditions_on_φ_axis!( fssrb.potential, iz, ir, rbi, nrbi, ax2, int2, fssrb.grid_boundary_factors[2], only2d)
            # end
            apply_boundary_conditions_on_r_axis!( fssrb.potential, iz, iφ, rbi, ax1, int1, fssrb.grid_boundary_factors[1])
        end
        @inbounds for ir in axes(fssrb.potential, 3)
            apply_boundary_conditions_on_cyl_z_axis!( fssrb.potential, ir, iφ, rbi, ax3, int3, fssrb.grid_boundary_factors[3])
        end
    else
        @inbounds for iz in axes(fssrb.potential, 1)
            for iφ in axes(fssrb.potential, 2)
                apply_boundary_conditions_on_r_axis!( fssrb.potential, iz, iφ, rbi, ax1, int1, fssrb.grid_boundary_factors[1])
            end
            for ir in axes(fssrb.potential, 3)
                apply_boundary_conditions_on_φ_axis!( fssrb.potential, iz, ir, rbi, nrbi, ax2, int2, fssrb.grid_boundary_factors[2], only2d)
            end
        end
        @inbounds for iφ in axes(fssrb.potential, 2) # z boundaries
            for ir in axes(fssrb.potential, 3)
                apply_boundary_conditions_on_cyl_z_axis!( fssrb.potential, ir, iφ, rbi, ax3, int3, fssrb.grid_boundary_factors[3])
            end
        end
        begin # r = 0 handling
            nφ::Int = size(fssrb.potential, 2) - 1
            gw_φ::Array{T, 2} = fssrb.geom_weights[2].weights
            @inbounds for inz in 1:(size(fssrb.ϵ, 3) - 1)
                m::T = 0
                l::T = 0
                for inφ in 1:nφ
                    if even_points ? isodd(inz + inφ) : iseven(inz + inφ)
                        l += gw_φ[3, inφ]
                        m += fssrb.potential[rbidx(inz), inφ + 1, 2, rbi] * gw_φ[3, inφ]
                    end
                end
                m *= inv(l)
                for inφ in 1:nφ
                    if even_points ? isodd(inz + inφ) : iseven(inz + inφ)
                        fssrb.potential[rbidx(inz), inφ + 1, 2, rbi]::T = m
                    end
                end
            end
        end
    end
    nothing
end


text:
function get_vecs!((phi, q), M, V, AV, ni)
  F = eigen(Hermitian(M))
  lambda = F.values[1]
  u = F.vectors[:, 1]
  phi .= u[1] .* V[1]
  q .= u[1] .* AV[1]
  for n in 2:ni
    phi .+= u[n] .* V[n]
    q .+= u[n] .* AV[n]
  end
  q .-= lambda .* phi
  #Fix sign
  if real(u[1]) < 0.0
    phi .*= -1
    q .*= -1
  end
  return lambda
end

function orthogonalize!(q::ITensor, V, ni)
  q0 = copy(q)
  for k in 1:ni
    Vq0k = dot(V[k], q0)
    q .+= -Vq0k .* V[k]
  end
  qnrm = norm(q)
  if qnrm < 1E-10 #orthog failure, try randomizing
    randn!(q)
    qnrm = norm(q)
  end
  q .*= 1.0 / qnrm
  return nothing
end

function expand_krylov_space(M::Matrix{elT}, V, AV, ni) where {elT}
  newM = fill(zero(elT), (ni + 1, ni + 1))
  newM[1:ni, 1:ni] = M
  for k in 1:(ni + 1)
    newM[k, ni + 1] = dot(V[k], AV[ni + 1])
    newM[ni + 1, k] = conj(newM[k, ni + 1])
  end
  return newM
end

function davidson(A, phi0::ITensorT; kwargs...) where {ITensorT<:ITensor}
  elTA = eltype(A)
  elTphi = eltype(phi0)

  # if the matrix is complex and the starting vector is real,
  # that's not going to last long
  # and we're going to need phi to be complex
  # down when we reuse the storage in get_vecs!
  if !(elTA <: Real) && (elTphi <: Real)
    phi = complex(phi0)
  else
    phi = copy(phi0)
  end

  maxiter = get(kwargs, :maxiter, 2)
  miniter = get(kwargs, :miniter, 1)
  errgoal = get(kwargs, :errgoal, 1E-14)
  Northo_pass = get(kwargs, :Northo_pass, 1)

  approx0 = 1E-12

  nrm = norm(phi)
  if nrm < 1E-18
    phi_ = similar(phi)
    randn!(phi_)
    phi = phi_
    nrm = norm(phi)
  end
  phi .*= 1.0 / nrm

  maxsize = size(A)[1]
  actual_maxiter = min(maxiter, maxsize - 1)

  if dim(inds(phi)) != maxsize
    error("linear size of A and dimension of phi should match in davidson")
  end

  Aphi = A(phi)

  V = ITensorT[copy(phi)]
  AV = ITensorT[Aphi]

  last_lambda = NaN
  lambda::Float64 = real(dot(V[1], AV[1]))
  q = AV[1] - lambda * V[1]

  M = fill(elTA(lambda), (1, 1))

  for ni in 1:actual_maxiter
    qnorm = norm(q)

    errgoal_reached = (qnorm < errgoal && abs(lambda - last_lambda) < errgoal)
    small_qnorm = (qnorm < max(approx0, errgoal * 1E-3))
    converged = errgoal_reached || small_qnorm

    if (qnorm < 1E-20) || (converged && ni > miniter) #|| (ni >= actual_maxiter)
      #@printf "  done with davidson, ni=%d, qnorm=%.3E\n" ni qnorm
      break
    end

    last_lambda = lambda

    for pass in 1:Northo_pass
      orthogonalize!(q, V, ni)
    end

    Aq = A(q)

    push!(V, copy(q))
    push!(AV, Aq)

    M = expand_krylov_space(M, V, AV, ni)

    lambda = get_vecs!((phi, q), M, V, AV, ni + 1)
  end #for ni=1:actual_maxiter+1

  return lambda, phi
end


text:
using Test
using TwoPointFunctions

import TwoPointFunctions: ∂θ₁∂θ₂cosχ, ∂θ₂∂θ₁cosχ, ∂ϕ₁∂ϕ₂cosχ, ∂ϕ₂∂ϕ₁cosχ,
∂ϕ₂∂ϕ₂cosχ, ∂ϕ₁∂ϕ₁cosχ, ∂θ₁∂ϕ₁cosχ, ∂θ₂∂ϕ₂cosχ, ∂ϕ₁∂θ₁cosχ, ∂ϕ₂∂θ₂cosχ

isapproxzero(x) = isapprox(x, zero(x), atol=1e-14, rtol=1e-8)

@testset "cosχ" begin
	# cos(θ₁)cos(θ₂) + sin(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂)
    @test cosχ((0, 0), (0, 0)) == 1
    @test cosχ((0, 0), (0, pi/3)) == 1
    @test cosχ((0, 0), (pi, 0)) == -1
    @test cosχ((0, 0), (pi, pi/3)) == -1

    @test cosχ((0, 0), (0,0)) == 1
    @test cosχ((0, pi/3), (0,0)) == 1
    @test cosχ((pi, 0), (0,0)) == -1
    @test cosχ((pi, pi/3), (0,0)) == -1

    @test isapproxzero(cosχ((0, 0), (pi/2, 0)))
    @test isapproxzero(cosχ((0, 0), (pi/2, pi/3)))

    @test isapproxzero(cosχ((pi/2, 0), (0, 0)))
    @test isapproxzero(cosχ((pi/2, pi/3), (0, 0)))

    @test cosχ((0, 0), (pi/3, 0)) ≈ cos(pi/3)
    @test cosχ((0, 0), (pi/3, pi/3)) ≈ cos(pi/3)
end

@testset "first derivatives" begin
    @testset "∂θ₁cosχ" begin
		# -sin(θ₁)cos(θ₂) + cos(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂)
	    @test isapproxzero(∂θ₁cosχ((0, 0), (0, 0)))
	    @test isapproxzero(∂θ₁cosχ((0, 0), (0, pi/3)))
	    @test isapproxzero(∂θ₁cosχ((0, 0), (pi, 0)))
	    @test isapproxzero(∂θ₁cosχ((0, 0), (pi, pi/3)))
	    @test ∂θ₁cosχ((0, 0), (pi/2, 0)) ≈ 1
	    @test ∂θ₁cosχ((0, 0), (pi/2, pi/3)) ≈ 0.5
	    @test ∂θ₁cosχ((pi/2, 0), (0, 0)) ≈ -1
	    @test ∂θ₁cosχ((pi/2, pi/3), (0, 0)) ≈ -1
	end

	@testset "∂θ₂cosχ" begin
		# -cos(θ₁)sin(θ₂) + sin(θ₁)cos(θ₂)cos(ϕ₁ - ϕ₂)
	    @test isapproxzero(∂θ₂cosχ((0, 0), (0, 0)))
	    @test isapproxzero(∂θ₂cosχ((0, 0), (0, pi/3)))
	    @test isapproxzero(∂θ₂cosχ((0, 0), (pi, 0)))
	    @test isapproxzero(∂θ₂cosχ((0, 0), (pi, pi/3)))
	    @test ∂θ₂cosχ((0, 0), (pi/2, 0)) ≈ -1
	    @test ∂θ₂cosχ((0, 0), (pi/2, pi/3)) ≈ -1
	    @test ∂θ₂cosχ((pi/2, 0), (0, 0)) ≈ 1
	    @test ∂θ₂cosχ((pi/2, pi/3), (0, 0)) ≈ 0.5
	end

	@testset "∂ϕ₁cosχ" begin
		# -sin(θ₁)sin(θ₂)sin(ϕ₁ - ϕ₂)
	    @test isapproxzero(∂ϕ₁cosχ((0, 0), (0, 0)))
	    @test isapproxzero(∂ϕ₁cosχ((0, 0), (0, pi/3)))
	    @test isapproxzero(∂ϕ₁cosχ((0, 0), (pi, 0)))
	    @test isapproxzero(∂ϕ₁cosχ((0, 0), (pi, pi/3)))

	    @test isapproxzero(∂ϕ₁cosχ((0, 0), (0, 0)))
	    @test isapproxzero(∂ϕ₁cosχ((0, pi/3), (0, 0)))
	    @test isapproxzero(∂ϕ₁cosχ((pi, 0), (0, 0)))
	    @test isapproxzero(∂ϕ₁cosχ((pi, pi/3), (0, 0)))

	    @test isapproxzero(∂ϕ₁cosχ((pi/2, 0), (pi/2, 0)))
	    @test ∂ϕ₁cosχ((pi/2, 0), (pi/2, pi/6)) ≈ 0.5
	    @test ∂ϕ₁cosχ((pi/2, pi/6), (pi/2, 0)) ≈ -0.5
	    @test ∂ϕ₁cosχ((pi/2, pi/2), (pi/2, 0)) ≈ -1
	end

	@testset "∇ϕ₁cosχ" begin
	    # -sin(θ₂)sin(ϕ₁ - ϕ₂)
	    @test isapprox(∇ϕ₁cosχ((0, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₁cosχ((0, 0), (0, pi/3)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₁cosχ((0, 0), (pi, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₁cosχ((0, 0), (pi, pi/3)), 0, atol=1e-14)

	    @test isapprox(∇ϕ₁cosχ((0, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₁cosχ((0, pi/3), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₁cosχ((pi, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₁cosχ((pi, pi/3), (0, 0)), 0, atol=1e-14)

	    @test isapprox(∇ϕ₁cosχ((pi/2, 0), (pi/2, 0)), 0, atol=1e-14)
	    @test ∇ϕ₁cosχ((pi/2, 0), (pi/2, pi/6)) ≈ 0.5
	    @test ∇ϕ₁cosχ((pi/2, pi/6), (pi/2, 0)) ≈ -0.5
	    @test ∇ϕ₁cosχ((pi/2, pi/2), (pi/2, 0)) ≈ -1
	end

	@testset "∂ϕ₂cosχ" begin
		# sin(θ₁)sin(θ₂)sin(ϕ₁ - ϕ₂)
	    @test isapprox(∂ϕ₂cosχ((0, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∂ϕ₂cosχ((0, 0), (0, pi/3)), 0, atol=1e-14)
	    @test isapprox(∂ϕ₂cosχ((0, 0), (pi, 0)), 0, atol=1e-14)
	    @test isapprox(∂ϕ₂cosχ((0, 0), (pi, pi/3)), 0, atol=1e-14)

	    @test isapprox(∂ϕ₂cosχ((0, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∂ϕ₂cosχ((0, pi/3), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∂ϕ₂cosχ((pi, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∂ϕ₂cosχ((pi, pi/3), (0, 0)), 0, atol=1e-14)

	    @test isapprox(∂ϕ₂cosχ((pi/2, 0), (pi/2, 0)), 0, atol=1e-14)
	    @test ∂ϕ₂cosχ((pi/2, 0), (pi/2, pi/6)) ≈ -0.5
	    @test ∂ϕ₂cosχ((pi/2, pi/6), (pi/2, 0)) ≈ 0.5
	    @test ∂ϕ₂cosχ((pi/2, pi/2), (pi/2, 0)) ≈ 1
	end

	@testset "∇ϕ₂cosχ" begin
	    # sin(θ₁)sin(ϕ₁ - ϕ₂)
	    @test isapprox(∇ϕ₂cosχ((0, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₂cosχ((0, 0), (0, pi/3)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₂cosχ((0, 0), (pi, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₂cosχ((0, 0), (pi, pi/3)), 0, atol=1e-14)

	    @test isapprox(∇ϕ₂cosχ((0, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₂cosχ((0, pi/3), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₂cosχ((pi, 0), (0, 0)), 0, atol=1e-14)
	    @test isapprox(∇ϕ₂cosχ((pi, pi/3), (0, 0)), 0, atol=1e-14)

	    @test isapprox(∇ϕ₂cosχ((pi/2, 0), (pi/2, 0)), 0, atol=1e-14)
	    @test ∇ϕ₂cosχ((pi/2, 0), (pi/2, pi/6)) ≈ -0.5
	    @test ∇ϕ₂cosχ((pi/2, pi/6), (pi/2, 0)) ≈ 0.5
	    @test ∇ϕ₂cosχ((pi/2, pi/2), (pi/2, 0)) ≈ 1
	end
end

@testset "second derivatives" begin
    @testset "∂²θ₁cosχ" begin
        # -cos(θ₁)cos(θ₂) -sin(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂) == -cosχ
        @test isapproxzero(∂²θ₁cosχ((pi/2, 0), (0, 0)))
        @test ∂²θ₁cosχ((pi/2, 0), (pi/2, 0)) == -1
        @test ∂²θ₁cosχ((pi/2, 0), (pi/2, pi/3)) ≈ -0.5

        @test ∂²θ₁cosχ((0, 0), (0, 0)) ≈ -cosχ((0, 0), (0, 0))
        @test ∂²θ₁cosχ((0, 0), (pi, 0)) ≈ -cosχ((0, 0), (pi, 0))
        @test ∂²θ₁cosχ((0, 0), (pi/2, 0)) ≈ -cosχ((0, 0), (pi/2, 0))
        @test ∂²θ₁cosχ((0, pi/3), (pi/2, 0)) ≈ -cosχ((0, pi/3), (pi/2, 0))
    end

    @testset "∂²θ₂cosχ" begin
        # -cos(θ₁)cos(θ₂) - sin(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂) = -cosχ
        @test isapproxzero(∂²θ₂cosχ((pi/2, 0), (0, 0)))
        @test ∂²θ₂cosχ((pi/2, 0), (pi/2, 0)) == -1
        @test ∂²θ₂cosχ((pi/2, 0), (pi/2, pi/3)) ≈ -0.5

        @test ∂²θ₂cosχ((0, 0), (0, 0)) ≈ -cosχ((0, 0), (0, 0))
        @test ∂²θ₂cosχ((0, 0), (pi, 0)) ≈ -cosχ((0, 0), (pi, 0))
        @test ∂²θ₂cosχ((0, 0), (pi/2, 0)) ≈ -cosχ((0, 0), (pi/2, 0))
        @test ∂²θ₂cosχ((0, pi/3), (pi/2, 0)) ≈ -cosχ((0, pi/3), (pi/2, 0))
    end

    @testset "∂θ₁∂θ₂cosχ" begin
    	# sin(θ₁)sin(θ₂) + cos(θ₁)cos(θ₂)cos(ϕ₁ - ϕ₂)

    	@test ∂θ₁∂θ₂cosχ((0, 0), (0, 0)) ≈ 1
    	@test ∂θ₁∂θ₂cosχ((0, pi/3), (0, 0)) ≈ cos(pi/3)
    	@test ∂θ₁∂θ₂cosχ((0, pi/3), (pi, 0)) ≈ -cos(pi/3)

    	@test isapproxzero(∂θ₁∂θ₂cosχ((pi/2, pi/3), (0, 0)))
    	@test isapproxzero(∂θ₁∂θ₂cosχ((pi/2, pi/3), (0, pi/3)))
    	@test isapproxzero(∂θ₁∂θ₂cosχ((pi/2, pi/3), (pi, 0)))
    	@test isapproxzero(∂θ₁∂θ₂cosχ((pi/2, 0), (pi, 0)))
    	@test isapproxzero(∂θ₁∂θ₂cosχ((pi/2, 0), (pi, pi/3)))

    	@test ∂θ₁∂θ₂cosχ((pi/2, 0), (pi/2, 0)) ≈ 1
    	@test ∂θ₁∂θ₂cosχ((pi/2, 0), (pi/2, pi/3)) ≈ 1

        @testset "∂θ₂∂θ₁cosχ" begin
        	for θ₁ in LinRange(0,π,10), ϕ₁ in LinRange(0,2π,20),
            	θ₂ in LinRange(0,π,10), ϕ₂ in LinRange(0,2π,20)

            	@test isapprox(∂θ₂∂θ₁cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)) , ∂θ₁∂θ₂cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)), atol=1e-14, rtol=1e-8)
        	end
        end
    end

    @testset "∂²ϕ₁cosχ" begin
        # -sin(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂)

        @test isapproxzero(∂²ϕ₁cosχ((0, 0), (0, 0)))
        @test isapproxzero(∂²ϕ₁cosχ((0, 0), (pi,0)))
        @test isapproxzero(∂²ϕ₁cosχ((pi, π/4), (pi,0)))
        @test isapproxzero(∂²ϕ₁cosχ((0, pi/3), (0, 0)))
        @test isapproxzero(∂²ϕ₁cosχ((0, pi/3), (0, pi/2)))
        @test isapproxzero(∂²ϕ₁cosχ((pi/4, pi/2), (pi/4, 0)))

        @test ∂²ϕ₁cosχ((pi/4, 0), (pi/4, 0)) ≈ -sin(π/4)^2
        @test ∂²ϕ₁cosχ((pi/4, pi/6), (pi/4, pi/6)) ≈ -sin(π/4)^2

        @test ∂²ϕ₁cosχ((pi/4, pi/3), (pi/4, pi/6)) ≈ -sin(π/4)^2*cos(π/6)

        @testset "∂ϕ₁∂ϕ₁cosχ" begin
        	for θ₁ in LinRange(0,π,10), ϕ₁ in LinRange(0,2π,20),
            	θ₂ in LinRange(0,π,10), ϕ₂ in LinRange(0,2π,20)

            	@test isapprox(∂²ϕ₁cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)) , ∂ϕ₁∂ϕ₁cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)), atol=1e-14, rtol=1e-8)
        	end
        end
    end

    @testset "∂²ϕ₂cosχ" begin
        # -sin(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂)

        @test isapproxzero(∂²ϕ₂cosχ((0, 0), (0, 0)))
        @test isapproxzero(∂²ϕ₂cosχ((0, 0), (pi,0)))
        @test isapproxzero(∂²ϕ₂cosχ((pi, π/4), (pi,0)))
        @test isapproxzero(∂²ϕ₂cosχ((0, pi/3), (0, 0)))
        @test isapproxzero(∂²ϕ₂cosχ((0, pi/3), (0, pi/2)))
        @test isapproxzero(∂²ϕ₂cosχ((pi/4, pi/2), (pi/4, 0)))

        @test ∂²ϕ₂cosχ((pi/4, 0), (pi/4, 0)) ≈ -sin(π/4)^2
        @test ∂²ϕ₂cosχ((pi/4, pi/6), (pi/4, pi/6)) ≈ -sin(π/4)^2

        @test ∂²ϕ₂cosχ((pi/4, pi/3), (pi/4, pi/6)) ≈ -sin(π/4)^2*cos(π/6)

        @testset "∂ϕ₂∂ϕ₂cosχ" begin
        	for θ₁ in LinRange(0,π,10), ϕ₁ in LinRange(0,2π,20),
            	θ₂ in LinRange(0,π,10), ϕ₂ in LinRange(0,2π,20)

            	@test isapprox(∂²ϕ₂cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)) , ∂ϕ₂∂ϕ₂cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)), atol=1e-14, rtol=1e-8)
        	end
        end
    end

    @testset "∂ϕ₁∂ϕ₂cosχ" begin
        # sin(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂)

        @test isapproxzero(∂ϕ₁∂ϕ₂cosχ((0, 0), (0, 0)))
        @test isapproxzero(∂ϕ₁∂ϕ₂cosχ((0, 0), (pi,0)))
        @test isapproxzero(∂ϕ₁∂ϕ₂cosχ((pi, π/4), (pi,0)))
        @test isapproxzero(∂ϕ₁∂ϕ₂cosχ((0, pi/3), (0, 0)))
        @test isapproxzero(∂ϕ₁∂ϕ₂cosχ((0, pi/3), (0, pi/2)))
        @test isapproxzero(∂ϕ₁∂ϕ₂cosχ((pi/4, pi/2), (pi/4, 0)))

        @test ∂ϕ₁∂ϕ₂cosχ((pi/4, 0), (pi/4, 0)) ≈ sin(π/4)^2
        @test ∂ϕ₁∂ϕ₂cosχ((pi/4, pi/6), (pi/4, pi/6)) ≈ sin(π/4)^2

        @test ∂ϕ₁∂ϕ₂cosχ((pi/4, pi/3), (pi/4, pi/6)) ≈ sin(π/4)^2*cos(π/6)

        @testset "∂ϕ₂∂ϕ₁cosχ" begin
        	for θ₁ in LinRange(0,π,10), ϕ₁ in LinRange(0,2π,20),
            	θ₂ in LinRange(0,π,10), ϕ₂ in LinRange(0,2π,20)

            	@test isapprox(∂ϕ₂∂ϕ₁cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)) , ∂ϕ₁∂ϕ₂cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)), atol=1e-14, rtol=1e-8)
        	end
        end
    end

    @testset "∂θ₁∂ϕ₁cosχ" begin

        @testset "∂ϕ₁∂θ₁cosχ" begin
        	for θ₁ in LinRange(0,π,10), ϕ₁ in LinRange(0,2π,20),
            	θ₂ in LinRange(0,π,10), ϕ₂ in LinRange(0,2π,20)

            	@test isapprox(∂θ₁∂ϕ₁cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)) , ∂ϕ₁∂θ₁cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)), atol=1e-14, rtol=1e-8)
        	end
        end
    end

    @testset "∂θ₂∂ϕ₂cosχ" begin

        @testset "∂ϕ₂∂θ₂cosχ" begin
        	for θ₁ in LinRange(0,π,10), ϕ₁ in LinRange(0,2π,20),
            	θ₂ in LinRange(0,π,10), ϕ₂ in LinRange(0,2π,20)

            	@test isapprox(∂θ₂∂ϕ₂cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)) , ∂ϕ₂∂θ₂cosχ((θ₁, ϕ₁), (θ₂, ϕ₂)), atol=1e-14, rtol=1e-8)
        	end
        end
    end
end


text:
function include_parameters(canonical::Canonical,
                            data::Array,
                            param_reference::UpdateRef,
                            axs...)
    if !model_has_parameters(canonical)
        error("Operational Model doesn't have parameters enabled. Include the keyword use_parameters=true")
    end
    param = _add_param_container!(canonical, param_reference, axs...)
    Cidx = CartesianIndices(length.(axs))
    for idx in Cidx
        param.data[idx] = PJ.add_parameter(canonical.JuMPmodel, data[idx])
    end
    return param
end

function include_parameters(canonical::Canonical,
                            ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                            param_reference::UpdateRef,
                            expression_name::Symbol,
                            multiplier::Float64 = 1.0)
    if !model_has_parameters(canonical)
        error("Operational Model doesn't have parameters enabled. Include the keyword use_parameters=true")
    end
    time_steps = model_time_steps(canonical)
    param = _add_param_container!(canonical, param_reference, (r[1] for r in ts_data), time_steps)
    expr = get_expression(canonical, expression_name)
    for t in time_steps, r in ts_data
        param[r[1], t] = PJ.add_parameter(canonical.JuMPmodel, r[4][t]);
        _add_to_expression!(expr, r[2], t, param[r[1], t], r[3] * multiplier)
    end
    return param
end


text:
mutable struct FirstReactionMethod <: CoupledSearchAlgorithm
  rates::Vector{Float64}
  total_rate::Float64
end

RatesCache(::FirstReactionMethod) = HasRates()

@inline search_jumps(algorithm::FirstReactionMethod) = search_jump_times(RatesCache(algorithm), algorithm.rates)

@inline update!(algorithm::FirstReactionMethod, state, model, j) = update_jump_rates!(algorithm, state, model)


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_storage_size_constraints(m, p, b; _n="")
    # TODO add formal types for storage (i.e. "b")

	# Constraint (4b)-1: Lower bound on Storage Energy Capacity
	@constraint(m,
        m[Symbol("dvStorageEnergy"*_n)][b] >= p.storage.min_kwh[b]
    )

	# Constraint (4b)-2: Upper bound on Storage Energy Capacity
	@constraint(m,
        m[Symbol("dvStorageEnergy"*_n)][b] <= p.storage.max_kwh[b]
    )

	# Constraint (4c)-1: Lower bound on Storage Power Capacity
	@constraint(m,
        m[Symbol("dvStoragePower"*_n)][b] >= p.storage.min_kw[b]
    )

	# Constraint (4c)-2: Upper bound on Storage Power Capacity
	@constraint(m,
        m[Symbol("dvStoragePower"*_n)][b] <= p.storage.max_kw[b]
    )
end

function add_storage_dispatch_constraints(m, p, b; _n="")
    # Constraint (4a): initial state of charge
	@constraint(m,
        m[Symbol("dvStoredEnergy"*_n)][b, 0] == p.storage.soc_init_pct[b] * m[Symbol("dvStorageEnergy"*_n)][b]
    )

	# Constraint (4g): state-of-charge for electrical storage - with grid
	@constraint(m, [ts in p.time_steps_with_grid],
        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_timestep * (
            sum(p.storage.charge_efficiency[b] * m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.elec_techs)
            + p.storage.grid_charge_efficiency[b] * m[Symbol("dvGridToStorage"*_n)][b, ts]
            - m[Symbol("dvDischargeFromStorage"*_n)][b,ts] / p.storage.discharge_efficiency[b]
        )
	)

	# Constraint (4h): state-of-charge for electrical storage - no grid
	@constraint(m, [ts in p.time_steps_without_grid],
        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_timestep * (
            sum(p.storage.charge_efficiency[b] * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.elec_techs)
            - m[Symbol("dvDischargeFromStorage"*_n)][b, ts] / p.storage.discharge_efficiency[b]
        )
    )

	# Constraint (4j): Minimum state of charge
	@constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoredEnergy"*_n)][b, ts] >= p.storage.soc_min_pct[b] * m[Symbol("dvStorageEnergy"*_n)][b]
    )

	# Constraint (4i)-1: Dispatch to electrical storage is no greater than power capacity
	@constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoragePower"*_n)][b] >=
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.elec_techs) + m[Symbol("dvGridToStorage"*_n)][b, ts]
    )

	#Constraint (4j): Dispatch from storage is no greater than power capacity
	@constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts]
    )

	#Constraint (4k)-alt: Dispatch to and from electrical storage is no greater than power capacity
	@constraint(m, [ts in p.time_steps_with_grid],
        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts] +
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.elec_techs) + m[Symbol("dvGridToStorage"*_n)][b, ts]
    )

	#Constraint (4l)-alt: Dispatch from electrical storage is no greater than power capacity
	@constraint(m, [ts in p.time_steps_without_grid],
        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] +
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.elec_techs)
    )

	#Constraint (4n): State of charge upper bound is storage system size
	@constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoredEnergy"*_n)][b,ts] <= m[Symbol("dvStorageEnergy"*_n)][b]
    )

    for b in setdiff(p.storage.types, p.storage.can_grid_charge)
        for ts in p.time_steps_with_grid
            fix(m[Symbol("dvGridToStorage"*_n)][b, ts], 0.0, force=true)
        end
	end
end

function add_storage_sum_constraints(m, p; _n="")

	##Constraint (8c): Grid-to-storage no greater than grid purchases
	@constraint(m, [ts in p.time_steps_with_grid],
        m[Symbol("dvGridPurchase"*_n)][ts] >= sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.storage.types)
    )
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
"""
    s2(array, phase[; len][, plans][, directions,] periodic = false)
    s2(array, SeparableIndicator(χ₁, χ₂)[; len][, plans][,directions,] periodic = false)
    s2(array, InseparableIndicator(χ)[; len][,directions,] periodic = false)

Calculate `S₂` (two point) correlation function for one-, two- or
three-dimensional multiphase system.

`S₂(x)` equals to probability that corner elements of a line segment
with the length `x` cut from the array belong to the same phase. This
implementation calculates `S₂(x)` for all `x`es in the range from `1`
to `len` which defaults to half of the minimal dimenstion of the
array.

More generally, you can provide indicator function `χ` instead of
`phase`. In this case `S₂` function calculates probability of `χ(x,
y)` returing `true` where `x` and `y` are two corners of a line
segment. Indicator functions must be wrapped in either
`SeparableIndicator` or `InseparableIndicator`. Some computations for
separable indicator functions are optimized.

An argument `plans` can be used to support precomputed FFT plans which
can be helpful if you call `s2` often with the array of the same
size. Plans can be computed with `S2FTPlans` constructor.

# Examples
```jldoctest
julia> s2([1,1,1,0,1,1], 1; len = 6)[:x]
6-element Array{Float64,1}:
 0.8333333333333334
 0.6
 0.5
 0.6666666666666666
 1.0
 1.0
```

See also: [`direction1Dp`](@ref), [`direction2Dp`](@ref),
[`direction3Dp`](@ref), [`SeparableIndicator`](@ref),
[`InseparableIndicator`](@ref), [`S2FTPlans`](@ref).
"""
function s2 end

maybe_pad_with_zeros(slice :: AbstractVector, :: Torus) = slice
maybe_pad_with_zeros(slice :: AbstractVector{T}, :: Plane) where T =
    let z = zeros(T, length(slice)); vcat(z, slice, z) end

expand_coefficient(:: Plane) = 3
expand_coefficient(:: Torus) = 1

struct S2FTPlans{F, I}
    forward :: Dict{Int, F}
    inverse :: Dict{Int, I}
end

Base.in(key :: Int, plans :: S2FTPlans) = key ∈ keys(plans.forward)

"""
    S2FTPlans(array, Plane())
    S2FTPlans(array, Torus())

Create FFT plans for faster calculation of two-point, surface-surface
and surface-void  correlation functions. `Plane()` corresponds to CW
boundary conditions and `Torus()` corresponds to periodic boundary
conditions.

See also: [`s2`](@ref), [`surfsurf`](@ref), [`surfvoid`](@ref).
"""
function S2FTPlans(array    :: AbstractArray,
                   topology :: Topology)
    m = expand_coefficient(topology)
    fft_plans  = Dict(s => zeros(Float64, m*s) |> plan_rfft for s in size(array))
    ifft_plans = Dict(s => plan_irfft(fft_plans[s] * zeros(Float64, m*s), m*s) for s in size(array))
    return S2FTPlans(fft_plans, ifft_plans)
end

"""
    S2FTPlans(array, periodic)

Create FFT plans for faster calculation of two-point, surface-surface
and surface-void  correlation functions. Periodic boundary conditions
are used when `periodic` is `true`, otherwise CW boundary conditions
are used.
"""
S2FTPlans(array :: AbstractArray, periodic :: Bool) =
    S2FTPlans(array, periodic ? Torus() : Plane())

function s2(array      :: AbstractArray,
            indicator  :: SeparableIndicator;
            len        :: Integer        = (array |> size |> minimum) ÷ 2,
            directions :: Vector{Symbol} = array |> default_directions,
            periodic   :: Bool           = false,
            plans      :: S2FTPlans      = S2FTPlans(array, periodic))
    cd = CorrelationData(len, check_directions(directions, size(array), periodic))
    topology = periodic ? Torus() : Plane()
    χ1, χ2 = indicator_function(indicator)

    for direction in directions
        success = cd.success[direction]
        total = cd.total[direction]
        slicer = slice_generators(array, periodic, Val(direction))

        for slice in slicer
            # Get plans for FFT and inverse FFT
            slen = length(slice)
            local fft, ifft

            ind1 = maybe_pad_with_zeros(map(χ1, slice), topology)
            ind2 = maybe_pad_with_zeros(map(χ2, slice), topology)

            if slen ∈ plans
                fft  = plans.forward[slen]
                ifft = plans.inverse[slen]
            else
                l = length(ind1)
                fft  = plan_rfft(ind1)
                ifft = plan_irfft(zeros(ComplexF64, l >> 1 + 1), l)
            end

            fft1 = fft * ind1
            fft2 = (χ1 === χ2) ? fft1 : fft * ind2
            s2 = ifft * (fft1 .* conj.(fft2))

            # Number of correlation lengths
            shifts = min(len, slen)

            # Update success and total
            success[1:shifts] .+= s2[1:shifts]
            if periodic
                total[1:shifts] .+= slen
            else
                update_runs!(total, slen, shifts)
            end
        end
    end

    return cd
end

function s2(array      :: AbstractArray,
            indicator  :: InseparableIndicator;
            len        :: Integer        = (array |> size |> minimum) ÷ 2,
            directions :: Vector{Symbol} = array |> default_directions,
            periodic   :: Bool           = false)
    cd = CorrelationData(len, check_directions(directions, size(array), periodic))
    χ = indicator_function(indicator)

    for direction in directions
        slicer = slice_generators(array, periodic, Val(direction))

        for slice in slicer
            slen = length(slice)
            # Number of shifts (distances between two points for this slice)
            shifts = min(len, slen)

            # Calculate slices where χ(slice[x]) && χ(slice[x+y]) for
            # all y's from 1 to len.
            cd.success[direction][1:shifts] .+= Iterators.map(1:shifts) do shift
                # Periodic slice, if needed
                pslice = periodic ? vcat(slice, slice[1:shift-1]) : slice
                plen = periodic ? slen+shift-1 : slen
                mapreduce(χ, +, pslice, view(pslice, shift:plen))
            end

            # Calculate total number of slices with lengths from 1 to len
            if periodic
                cd.total[direction][1:shifts] .+= slen
            else
                update_runs!(cd.total[direction], slen, shifts)
            end
        end
    end

    return cd
end

s2(array      :: AbstractArray,
   phase;
   len        :: Integer        = (array |> size |> minimum) ÷ 2,
   directions :: Vector{Symbol} = array |> default_directions,
   periodic   :: Bool           = false,
   plans      :: S2FTPlans      = S2FTPlans(array, periodic)) =
       s2(array, SeparableIndicator(x -> x == phase);
          len        = len,
          directions = directions,
          periodic   = periodic,
          plans      = plans)


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

using Markdown
using InteractiveUtils

# ╔═╡ c10a143e-3226-11ec-06d1-13d8fedcc6e3
# hideall

using Genie;

# ╔═╡ e608c9a3-3875-48f3-b28c-991e05534072
md"""
# Adding your existing Julia code into Genie apps

If you have existing Julia code (modules and libraries) which you'd like to quickly expose on the web without building an app from scratch, Genie provides an easy way to add and load your code into an app.

## Adding your Julia code to a Genie app

If you have an existing Julia application or standalone codebase which you'd like to expose over the web through a Genie app, the simplest thing is to add the files to the `lib/` folder. The `lib/` folder is automatically added by Genie to the `LOAD_PATH`, _recursively_.

This means that you can also add folders under `lib/`, and they will be recursively added to the `LOAD_PATH`. Beware though that this only happens when the Genie app is initially loaded. Hence, an app restart might be required if you add nested folders after the app is started.
"""

# ╔═╡ 2589c03d-fd17-4ba3-bdf9-90c4d779a136
md"""
Creating a simple Genie web service:
"""

# ╔═╡ 03fce09f-0f05-47bd-960b-eed9ef1087dd
# hideall

Genie.newapp_webservice("MyGenieApp");

# ╔═╡ c7ca2f4a-09d4-48d3-bb96-1fea6195b9cf
md"""
**HEADS UP**

In most cases, Genie won't create the `lib/` folder by default. If the `lib/` folder is not present in the root of the app, just create it yourself:
"""

# ╔═╡ 47574366-5c9a-47be-ac60-6a93a4ab9d81
md"""

```julia
julia> mkdir("lib")
```
"""

# ╔═╡ 4f7496dd-fa87-4180-8d6e-dfe663d66cc4
# hideall

begin
	incPath= true
	mkdir("lib")
end

# ╔═╡ b439eaa1-f4e9-40a5-9671-8684d5e7bd2c
# hideall

begin
	incPath;
	include("routes.jl");
end;

# ╔═╡ 31191d7b-45c9-434f-8be5-9f0088859ce5
# hideall

begin
	incPath
	@assert ispath("lib") == true
end

# ╔═╡ 22c1473d-c623-475b-9a18-41f5977c02c9
md"""
Once your code is added to the `lib/` folder, it will become available in your app's environment. For example, say we have a file called `lib/MyLib.jl`:

```julia
# lib/MyLib.jl
module MyLib

using Dates

function isitfriday()
  Dates.dayofweek(Dates.now()) == Dates.Friday
end

end
```

Then we can reference it in `routes.jl` and expose it on the web as follows:

"""

# ╔═╡ 07e003a6-297e-486a-b8dc-acd18cf15d0c
# hideall

begin
	incPath;
	touch("lib/MyLib.jl");
end;

# ╔═╡ c6427379-f6fb-4a3e-8be7-e3d1b9e3ad52
# hideall

begin
	incPath;
	@assert isfile("lib/MyLib.jl") == true;
end;

# ╔═╡ de394ad3-d9f3-44da-89e0-411d38c72deb
# hideall

begin
	incPath;
	write("lib/MyLib.jl", """
		# lib/MyLib.jl
	module MyLib

	using Dates

	function isitfriday()
	  Dates.dayofweek(Dates.now()) == Dates.Friday
	end

	end""");
end;

# ╔═╡ 8f6ffc7c-0102-4682-8779-ccba834982c9
md"""

```julia
# routes.jl
using Genie.Router
using MyLib

route("/friday") do
  MyLib.isitfriday() ? "Yes, it's Friday!" : "No, not yet :("
end
```
"""

# ╔═╡ 394e18be-97d6-4456-8f7c-3f5c13d47d12
# hideall

Genie.Router.routes();

# ╔═╡ d9a2ee12-f679-44c7-bd9b-0b83192e37a7
md"""
Use the `lib/` folder to host your Julia code so that Genie knows where to look in order to load it and make it available throughout the application.
"""

# ╔═╡ de8ba170-f425-4cbf-a21c-8eb926fbe683
# hideall

begin
	incPath;
	@assert isfile("routes.jl") == true;
end;

# ╔═╡ 0e49172c-a2e1-4499-a2bf-672daa13f9b9
md"""
Adding above code snippet to `routes.jl`
"""

# ╔═╡ dd5a7ecd-a950-478c-b458-4e138f41d89d
# hideall

write("routes.jl", """
	# routes.jl
using Genie.Router
using MyLib

route("/friday") do
  MyLib.isitfriday() ? "Yes, it's Friday!" : "No, not yet :("
end""");

# ╔═╡ be30e358-da0f-42d5-b984-d3e80fe16c65
# hideall
down();

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"

[compat]
Genie = "~4.1.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgParse]]
deps = ["Logging", "TextWrap"]
git-tree-sha1 = "3102bce13da501c9104df33549f511cd25264d7d"
uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
version = "1.1.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSTParser]]
deps = ["Tokenize"]
git-tree-sha1 = "f9a6389348207faf5e5c62cbc7e89d19688d338a"
uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f"
version = "3.3.0"

[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "393ac9df4eb085c2ab12005fc496dae2e1da344e"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.3"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "618835ab81e4a40acf215c98768978d82abc5d97"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.16"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

[[Genie]]
deps = ["ArgParse", "Dates", "Distributed", "EzXML", "FilePathsBase", "HTTP", "HttpCommon", "Inflector", "JSON3", "JuliaFormatter", "Logging", "Markdown", "MbedTLS", "Millboard", "Nettle", "OrderedCollections", "Pkg", "REPL", "Reexport", "Revise", "SHA", "Serialization", "Sockets", "UUIDs", "Unicode", "VersionCheck", "YAML"]
git-tree-sha1 = "1cb249c696d4c81d1dde38277602fbad6dc16594"
uuid = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
version = "4.1.0"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[HttpCommon]]
deps = ["Dates", "Nullables", "Test", "URIParser"]
git-tree-sha1 = "46313284237aa6ca67a6bce6d6fbd323d19cff59"
uuid = "77172c1b-203f-54ac-aa54-3f1198fe9f90"
version = "0.5.0"

[[Inflector]]
deps = ["Unicode"]
git-tree-sha1 = "8555b54ddf27806b070ce1d1cf623e1feb13750c"
uuid = "6d011eab-0732-4556-8808-e463c76bf3b6"
version = "1.0.1"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[JuliaFormatter]]
deps = ["CSTParser", "CommonMark", "DataStructures", "Pkg", "Tokenize"]
git-tree-sha1 = "10c95cebcfa37c1f510a726c90886db4745e1238"
uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
version = "0.15.11"

[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Millboard]]
git-tree-sha1 = "ea6a5b7e56e76d8051023faaa11d91d1d881dac3"
uuid = "39ec1447-df44-5f4c-beaa-866f30b4d3b2"
version = "0.2.5"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[Nettle]]
deps = ["Libdl", "Nettle_jll"]
git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
version = "0.5.1"

[[Nettle_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
version = "3.7.2+0"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "41deb3df28ecf75307b6e492a738821b031f8425"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.1.20"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StringEncodings]]
deps = ["Libiconv_jll"]
git-tree-sha1 = "50ccd5ddb00d19392577902f0079267a72c5ab04"
uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
version = "0.3.5"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TextWrap]]
git-tree-sha1 = "9250ef9b01b66667380cf3275b3f7488d0e25faf"
uuid = "b718987f-49a8-5099-9789-dcd902bef87d"
version = "1.0.1"

[[Tokenize]]
git-tree-sha1 = "0952c9cee34988092d73a5708780b3917166a0dd"
uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624"
version = "0.5.21"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[UrlDownload]]
deps = ["HTTP", "ProgressMeter"]
git-tree-sha1 = "05f86730c7a53c9da603bd506a4fc9ad0851171c"
uuid = "856ac37a-3032-4c1c-9122-f86d88358c8b"
version = "1.0.0"

[[VersionCheck]]
deps = ["Dates", "JSON3", "Logging", "Pkg", "Random", "Scratch", "UrlDownload"]
git-tree-sha1 = "89ef2431dd59344ebaf052d0737205854ded0c62"
uuid = "a637dc6b-bca1-447e-a4fa-35264c9d0580"
version = "0.2.0"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[YAML]]
deps = ["Base64", "Dates", "Printf", "StringEncodings"]
git-tree-sha1 = "3c6e8b9f5cdaaa21340f841653942e1a6b6561e5"
uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
version = "0.4.7"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─e608c9a3-3875-48f3-b28c-991e05534072
# ╠═c10a143e-3226-11ec-06d1-13d8fedcc6e3
# ╟─2589c03d-fd17-4ba3-bdf9-90c4d779a136
# ╠═03fce09f-0f05-47bd-960b-eed9ef1087dd
# ╟─c7ca2f4a-09d4-48d3-bb96-1fea6195b9cf
# ╟─47574366-5c9a-47be-ac60-6a93a4ab9d81
# ╠═4f7496dd-fa87-4180-8d6e-dfe663d66cc4
# ╠═31191d7b-45c9-434f-8be5-9f0088859ce5
# ╟─22c1473d-c623-475b-9a18-41f5977c02c9
# ╠═07e003a6-297e-486a-b8dc-acd18cf15d0c
# ╠═c6427379-f6fb-4a3e-8be7-e3d1b9e3ad52
# ╠═de394ad3-d9f3-44da-89e0-411d38c72deb
# ╟─8f6ffc7c-0102-4682-8779-ccba834982c9
# ╠═394e18be-97d6-4456-8f7c-3f5c13d47d12
# ╠═d9a2ee12-f679-44c7-bd9b-0b83192e37a7
# ╠═de8ba170-f425-4cbf-a21c-8eb926fbe683
# ╟─0e49172c-a2e1-4499-a2bf-672daa13f9b9
# ╠═dd5a7ecd-a950-478c-b458-4e138f41d89d
# ╠═b439eaa1-f4e9-40a5-9671-8684d5e7bd2c
# ╠═be30e358-da0f-42d5-b984-d3e80fe16c65
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
    Shuffle( X, Y )

Shuffles the rows of the `X` and `Y` data without replacement. It returns a 2-Tuple of the shuffled set.
"""
function Shuffle( X, Y )
    Inds = StatsBase.sample( collect(1 : size(X)[1] ), size(X)[1], replace = false )
    return( X[Inds,:], Y[Inds,:] )
end

"""
    Shuffle!( X, Y )

Shuffles the rows of the `X` and `Y` data without replacement in place. In place, means that this function
alters the order of the data in memory and this function does not return anything.
"""
function Shuffle!( X, Y )
    Inds = StatsBase.sample( collect(1 : size(X)[1] ), size(X)[1], replace = false )
    X .= X[Inds,:]
    Y .= Y[Inds,:]
end

struct KFoldsValidation
    K::Int
    FoldSize::Int
    observations::Int
    remainder::Int
    X
    Y
end

"""
    KFoldsValidation(K::Int, x, y)

Returns a KFoldsValidation iterator with `K` folds. Because it's an iterator it can be used in for loops,
see the tutorials for pragmatic examples. The iterator returns a 2-Tuple of 2-Tuples which have the
 following form: ``` ((TrainX,TrainY),(ValidateX,ValidateY) ```.

"""
function KFoldsValidation(K::Int, x, y)
    observations = size( x )[ 1 ]
    if floor(size( x )[ 1 ] / K) == 1
        K = observations
        println("Warning: Hold out size is 1, defaulting to Leave One Out Validation")
    end
    Remainder = observations % K
    KFoldsValidation( K, floor(size( x )[ 1 ] / K), observations, Remainder, x, forceMatrix(y))
end

"""
    LeaveOneOut(x, y)

Returns a KFoldsValidation iterator with leave one out folds. Because it's an iterator it can be used in for loops,
see the tutorials for pragmatic examples. The iterator returns a 2-Tuple of 2-Tuples which have the
 following form: ``` ((TrainX,TrainY),(ValidateX,ValidateY) ```.
"""
LeaveOneOut(x, y) = KFoldsValidation( size( x )[ 1 ], 1, size( x )[ 1 ], 0, x, forceMatrix(y) )

function Base.iterate( iter::KFoldsValidation, state = (0,0) )
    ( i, o ) = state;
    predictInds = []; trainInds = []; holdoutInds = []
    r = ( ( iter.K - i) <= iter.remainder) ? 1 : 0

    if i == 0
        holdoutInds = collect( 1 : iter.FoldSize )
        trainInds = collect( ( iter.FoldSize + 1 ) : iter.observations )
    elseif i < iter.K
        HOBegin = ( i * iter.FoldSize ) + 1 + (o)
        HOEnd = ( i + 1 ) * iter.FoldSize + o + r
        TBegin = collect( 1 : (HOBegin - 1) )
        TEnd = collect( ( HOEnd + 1 ) : iter.observations )
        holdoutInds = collect( HOBegin : HOEnd )
        trainInds = [ TBegin; TEnd]
    else
        return nothing
    end
    return (    (   ( iter.X[trainInds,:]     , iter.Y[trainInds,:] ),
                    ( ( iter.X[holdoutInds,:]   , iter.Y[holdoutInds,:] ) )
                ),
                ( i + 1 , o + r) )
end


text:
using ProjectFlow

println("Testing Exports !")

p = Project(
    id="xyz",
    name="My Fancy? *Project1 2 ",
    template="jl",
    profile="default"
)

datalake, iviz, idata = initiate(p);


text:
using CircuitSimulator: ParsedCircuit, MNA, MNAbuilder,
                        parse_netlist,nodedict_type, Resistor,
                        Inductor,Capacitor, CurrentSource, VoltageSource
@testset "parser Tests" begin
    @test typeof(ParsedCircuit{Float64}()) == ParsedCircuit{Float64}
    @test typeof(MNAbuilder{Float64}()) == MNAbuilder{Float64}
    netlist1=IOBuffer(
    """Hello World
    R1 a b 10k a=1
    R2 a 0 10 bgh s + 1
    Lone 1 2 23.7Meg b = 4 c = 5
    C10 1 a 2
    V1 b 0 2
    I1 2 b 3
    Z4 z x 45 \$ not a valid part
    """)
    pc = parse_netlist(netlist1) # add @test_warn here someday
    @test pc.titleline == "Hello World"
    @test pc.nodedict==nodedict_type(Symbol("a")=>1,Symbol("b")=>2,Symbol(1)=>3,Symbol(2)=>4,Symbol(0)=>0)
    @test pc.max_node == maximum(collect(values(pc.nodedict)))
    @test pc.netlist[1] == Resistor(:R1,(1,2),10.0e3)
    @test pc.netlist[2] == Resistor(:R2,(1,0),10.0)
    @test pc.netlist[3] == Inductor(:Lone,(3,4),23.7e6)
    @test pc.netlist[4] == Capacitor(:C10,(3,1),2.0)
    @test pc.netlist[5] == VoltageSource(:V1,(2,0),2.0)
    @test pc.netlist[6] == CurrentSource(:I1,(4,2),3.0)
    @test pc.max_element == length(pc.netlist)
    netlist2=IOBuffer(
    """netlist2
    R1 a b 10k + V(a) + V(b,1)
    R1 a b 10k + V(a) + V(b,1) a = 1 b = 34
    R1 a b 10k + V(a) + V(b,1) a 1 b 34
    R1 c d I(R1) - 7.9megOhm
    """)
    pc2 = parse_netlist(netlist2)
    eval(:(@inline netvoltage(::Val{:a}) = 5e3))
    eval(:(@inline netvoltage(::Val{:b}) = 20e3))
    eval(:(@inline netvoltage(::Val{Symbol("1")}) = 40e3))
    eval(:(@inline branchcurrent(::Val{:R1}) = 4.3e6))
    eval(pc2.netlist[1].value)
    @test eval(pc2.netlist[1].value) ≈ -5.00e3
    @test eval(pc2.netlist[2].value) ≈ -5.00e3
    @test eval(pc2.netlist[3].value) ≈ -5.00e3
    @test eval(pc2.netlist[4].value) ≈ -3.6e6
end


text:
#EasyPlot Defaults control
#-------------------------------------------------------------------------------

#=Design philosophy:
Defaults are initialized from ENV[] variables, so that they can be set from
~/.julia/config/startup.jl without first loading the EasyPlot module.
=#

#==Types
===============================================================================#
#Default width/height for rendering inline plots (preferably low res):
const DEFAULT_RENDERW = 300
const DEFAULT_RENDERH = round(Int, DEFAULT_RENDERW / MathConstants.φ) #φ: golden ratio

#==Types
===============================================================================#
mutable struct Defaults
	rendersvg::Bool #Might want to dissalow SVG renderings for performance reasons
	mimeshowopt::ShowOptions

	#When no builder is specified:
	guibuilder::Optional{AbstractBuilder}
	filebuilder::Optional{AbstractBuilder} #Writing to file
	mimebuilder::Optional{AbstractBuilder} #Responds directly to show(::IO, ::MIME, ...)
end
Defaults() = Defaults(false, ShowOptions(), nothing, nothing, nothing)

#==Data
===============================================================================#
const global defaults = Defaults()

#==Helper functions
===============================================================================#
function readdefaults(::Type{AbstractPlotDimensions}, envstr::String, default::AbstractPlotDimensions)
	result = default
	optstr = "{(w, h), :default}"
	val = get(ENV, envstr, nothing)
	if isa(val, Tuple{Int,Int})
		return PlotDim(val...)
	elseif isnothing(val)
		#use default
	elseif s == :auto
		result = plotautosize
	else
		@warn("$envstr valid settings are: $optstr")
	end

	return result
end

function readdefaults(::Type{Bool}, envstr::String, default::Bool)
	bstr = ["FALSE", "TRUE"] #WANTCONST
	val = get(ENV, envstr, string(default))
	uval = uppercase(val)

	if !(in(uval, bstr))
		optstr = join(bstr, ", ")
		@warn("$envstr valid settings are: $optstr")
	end

	#Return default if not recognized as !default:
	return bstr[!default+1] != uval ? default : !default
end

#==Registration functions
===============================================================================#
function register(displayid::Symbol)
	#global defaults
	#TODO
end

#==Initialization
===============================================================================#
function _initialize(dflt::Defaults)
	dflt.rendersvg = readdefaults(Bool, "EASYPLOT_RENDERSVG", true)
	dim = readdefaults(AbstractPlotDimensions, "EASYPLOT_MIMESIZE", plotautosize)
	dflt.mimeshowopt = ShowOptions(dim)
	return
end

#Last line


text:
using Animations
using Javis
using NBodySimulator
using StaticArrays

# Define the bodies
body1 = MassBody(SVector(0.0, 1.0, 0.0), SVector(5.775e-6, 0.0, 0.0), 2.0)
body2 = MassBody(SVector(0.0, -1.0, 0.0), SVector(-5.775e-6, 0.0, 0.0), 2.0)

G = 6.673e-11

system = GravitationalSystem([body1, body2], G)
tspan = (0.0, 1111150.0)

# run the simulation
simulation = NBodySimulation(system, tspan)
sim_result = run_simulation(simulation, saveat = 10)
sim_data = sim_result.solution.u

# extract the data
planet1_pos = [(sim_data[x][1, 1], sim_data[x][2, 1]) for x in 1:length(sim_data)]
planet2_pos = [(sim_data[x][1, 2], sim_data[x][2, 2]) for x in 1:length(sim_data)]

#Convert to Points
planet1_point =
    Point.((200 .* getindex.(planet1_pos, 1)), ((200 .* getindex.(planet1_pos, 2) .- 100)))
planet2_point =
    Point.((200 .* getindex.(planet2_pos, 1)), ((200 .* getindex.(planet2_pos, 2) .+ 100)))

# Start the animation
function ground(args...)
    background("black")
    sethue("white")
end

function planet(p, r, color)
    sethue(color)
    return circle(p, r, :fill)
end

function make_animation()
    nframes = 200

    vid = Video(500, 500)
    Background(1:nframes, ground)

    planet1 = Object(1:nframes, (args...) -> planet(planet1_point[1], 10.0, "red"))
    planet2 = Object(1:nframes, (args...) -> planet(planet2_point[1], 10.0, "blue"))

    planet1_anim = Animation(
        collect(range(0, stop = 1, length = length(planet1_point))),
        planet1_point,
        [linear() for _ in 1:(length(planet1_point) - 1)],
    )

    planet2_anim = Animation(
        collect(range(0, stop = 1, length = length(planet2_point))),
        planet2_point,
        [linear() for _ in 1:(length(planet2_point) - 1)],
    )

    act!(planet1, Action(1:nframes, planet1_anim, translate()))
    act!(planet2, Action(1:nframes, planet2_anim, translate()))

    render(vid; pathname = "two_body.gif")
end

make_animation()

max_stars_repo_name: Ved-Mahajan/graphic-dev

max_issues_repo_name: Ved-Mahajan/graphic-dev

max_forks_repo_name: Ved-Mahajan/graphic-dev


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

mutable struct OSDiskImage <: SwaggerModel
    operatingSystem::Any # spec type: Union{ Nothing, String } # spec name: operatingSystem

    function OSDiskImage(;operatingSystem=nothing)
        o = new()
        validate_property(OSDiskImage, Symbol("operatingSystem"), operatingSystem)
        setfield!(o, Symbol("operatingSystem"), operatingSystem)
        o
    end
end # type OSDiskImage

const _property_map_OSDiskImage = Dict{Symbol,Symbol}(Symbol("operatingSystem")=>Symbol("operatingSystem"))
const _property_types_OSDiskImage = Dict{Symbol,String}(Symbol("operatingSystem")=>"String")
Base.propertynames(::Type{ OSDiskImage }) = collect(keys(_property_map_OSDiskImage))
Swagger.property_type(::Type{ OSDiskImage }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OSDiskImage[name]))}
Swagger.field_name(::Type{ OSDiskImage }, property_name::Symbol) =  _property_map_OSDiskImage[property_name]

const _allowed_OSDiskImage_operatingSystem = ["Windows", "Linux"]

function check_required(o::OSDiskImage)
    true
end

function validate_property(::Type{ OSDiskImage }, name::Symbol, val)
    if name === Symbol("operatingSystem")
        Swagger.validate_param(name, "OSDiskImage", :enum, val, _allowed_OSDiskImage_operatingSystem)
    end
end


text:
"""
    Checkpointer(folder)

Saves `learner.model` to `folder` after every [`AbstractTrainingPhase`].

Use `FluxTraining.`[`loadmodel`](#) to load a model.
"""
struct Checkpointer <: Callback
    folder
    function Checkpointer(folder)
        mkpath(folder)
        return new(folder)
    end
end

stateaccess(::Checkpointer) = (
    model = Read(),
    cbstate = (metricsepoch = Read(), history = Read())
)

function on(::EpochEnd, phase::AbstractTrainingPhase, checkpointer::Checkpointer, learner)
    loss = last(learner.cbstate.metricsepoch[phase], :Loss)[2]
    epoch = learner.cbstate.history.epochs
    filename = "checkpoint_epoch_$(lpad(string(epoch), 3, '0'))_loss_$loss.bson"
    savemodel(learner.model, joinpath(checkpointer.folder, filename))
end

# TODO: replace with JLD2?
function savemodel(model, path)
    @save path model = cpu(model)
end

"""
    loadmodel(path)

Loads a model that was saved to `path` using `FluxTraining.`[`savemodel`](#).
"""
function loadmodel(path)
    @load path model
    return model
end


text:
using CombUtils

function main()
    a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    for i = 2:1_000_000
        next_permutation!(a)
    end
    println(reduce(*, map((n) -> string(n), a)))
end

main()