text:
module Algos

include("./Environments.jl");
using .Environments
using Setfield
using Plots
using Distributions

export QLearning, Sarsa, train, evaluate

DEFAULT_EPSILON = 0.8
DEFAULT_DISCOUNT = 0.9
DEFAULT_LEARNING_RATE = 0.1
ACTIONS = ("left", "upper-left", "upper-right", "right", "lower-right", "lower-left")
ACTIONS_NUM = length(ACTIONS)

struct QLearning
    qTable::Array{Float64, 3}
    learningRate::Float64
    discountFactor::Float64
    epsilon::Float64
    env::HexEnv
end;

function QLearning(tiles::Matrix{Int64}; learningRate::Float64 = DEFAULT_LEARNING_RATE,
                    discountFactor::Float64 = DEFAULT_DISCOUNT, epsilon::Float64 = DEFAULT_EPSILON)

    qTableSize = (findmax(tiles[:, 1])[1], findmax(tiles[:, 2])[1]);
    qTable = rand(Uniform(-1., 0.), tuple(vcat(collect(qTableSize), [ACTIONS_NUM])...));
    # qTable = zeros(tuple(vcat(collect(qTableSize), [ACTIONS_NUM])...));
    env = HexEnv(tiles);
    return QLearning(qTable, learningRate, discountFactor, epsilon, env);
end;

function train(agent::QLearning, episodes::Int64; visualizeEps = nothing, gifFile = nothing, prob::Bool = true, samplingRate::Float64 = 0.001)

    print("\nStep cost: $(Environments.STEP_COST) \n")
    print("Out of bounds penalty: $(Environments.OUT_OF_BOUNDS_PENALTY) \n")

    EPSILON_DECAY_START_EP = 1
    EPSILON_DECAY_END_EP = episodes // 2;
    EPSILON_DECAY = agent.epsilon / EPSILON_DECAY_END_EP;
    SAMPLING_EP = trunc(Int, 1 / samplingRate);

    # accumulated reward
    y = [];
    accumulatedReward = 0;

    anim = Environments.Plots.Animation()

    if visualizeEps === nothing
        visualizeEps = episodes + 1
    end

    for episode in 1:episodes

        if (episode % visualizeEps == 0)
            print("Visualizing episode $episode/$episodes\r");
            flush(stdout);
            if (gifFile !== nothing)
                visualize(agent.env, show = false);
                Environments.Plots.frame(anim);
            else
                visualize(agent.env, lowDpi = true);
                sleep(0.5);
            end;
            visualizeEpsLabel(episode);
        end;

        terminal = false
        agent = @set agent.env = Environments.reset(agent.env)

        while (!terminal)

            nextAction = getNextAction(agent.env.currentState, agent.epsilon, agent.qTable, prob = prob)

            env, nextState, reward, terminal = Environments.step(agent.env, nextAction)
            accumulatedReward += reward;

            currentQValue = agent.qTable[env.currentState[1], env.currentState[2], nextAction]

            newQValue = (1 - agent.learningRate) * currentQValue + agent.learningRate *
                (reward + agent.discountFactor * findmax(agent.qTable[nextState[1], nextState[2], :])[1])

            agent.qTable[env.currentState[1], env.currentState[2], nextAction] = newQValue

            if (episode % visualizeEps == 0)
                if (gifFile !== nothing)
                    visualize(agent.env.currentState, show = false);
                    Environments.Plots.frame(anim);
                else
                    visualize(agent.env.currentState);
                    sleep(0.1);
                end
            end

            env = @set env.currentState = nextState

            agent = @set agent.env = env

        end

        if (episode <= EPSILON_DECAY_END_EP && episode >= EPSILON_DECAY_START_EP)
            agent = @set agent.epsilon = agent.epsilon - EPSILON_DECAY
        end

        if (episode % visualizeEps == 0)
            if (gifFile === nothing)
                sleep(1);
            end;
        end;

        if (visualizeEps > episodes && episode % SAMPLING_EP == 0)
            push!(y, accumulatedReward / episode);
            print("Trained $episode/$episodes episodes\r");
            flush(stdout);
        end

    end;

    print("\n");

    if (gifFile !== nothing)
        print("Writing GIF\n");
        Environments.Plots.gif(anim, gifFile, fps = 10);
    else
        if (visualizeEps <= episodes)
            print("Press any key to continue...\n");
            readline();
        end
    end;

    return y;
end;

struct Sarsa
    qTable::Array{Float64, 3}
    learningRate::Float64
    discountFactor::Float64
    epsilon::Float64
    env::HexEnv
end;

function Sarsa(tiles::Matrix{Int64}; learningRate::Float64 = DEFAULT_LEARNING_RATE,
                    discountFactor::Float64 = DEFAULT_DISCOUNT, epsilon::Float64 = DEFAULT_EPSILON)

    qTableSize = (findmax(tiles[:, 1])[1], findmax(tiles[:, 2])[1]);
    qTable = rand(Uniform(-1., 0.), tuple(vcat(collect(qTableSize), [ACTIONS_NUM])...));
    env = HexEnv(tiles);
    return Sarsa(qTable, learningRate, discountFactor, epsilon, env);
end;

function train(agent::Sarsa, episodes::Int64; visualizeEps = nothing, gifFile = nothing, prob::Bool = true, samplingRate::Float64 = 0.001)

    print("Step cost: $(Environments.STEP_COST) \n")
    print("Out of bounds penalty: $(Environments.OUT_OF_BOUNDS_PENALTY) \n")

    EPSILON_DECAY_START_EP = 1
    EPSILON_DECAY_END_EP = episodes // 2;
    EPSILON_DECAY = agent.epsilon / EPSILON_DECAY_END_EP;
    SAMPLING_EP = trunc(Int, 1 / samplingRate);

    # accumulated reward
    y = [];
    accumulatedReward = 0;

    anim = Environments.Plots.Animation()

    if visualizeEps === nothing
        visualizeEps = episodes + 1
    end

    for episode in 1:episodes

        if (episode % visualizeEps == 0)
            print("Visualizing episode $episode/$episodes\r");
            flush(stdout);
            if (gifFile !== nothing)
                visualize(agent.env, show = false);
                Environments.Plots.frame(anim);
            else
                visualize(agent.env, lowDpi = true);
                sleep(0.5);
            end;
            visualizeEpsLabel(episode);
        end;

        terminal = false
        agent = @set agent.env = Environments.reset(agent.env)
        nextAction = getNextAction(agent.env.currentState, agent.epsilon, agent.qTable, prob = prob)

        while (!terminal)

            if (episode % visualizeEps == 0)
                if (gifFile !== nothing)
                    visualize(agent.env.currentState, show = false);
                    Environments.Plots.frame(anim);
                else
                    visualize(agent.env.currentState);
                    sleep(0.1);
                end
            end

            env, nextState, reward, terminal = Environments.step(agent.env, nextAction)
            accumulatedReward += reward;

            currentQValue = agent.qTable[env.currentState[1], env.currentState[2], nextAction]

            nextNextAction = getNextAction(nextState, agent.epsilon, agent.qTable, prob = prob)

            newQValue = (1 - agent.learningRate) * currentQValue + agent.learningRate *
                (reward + agent.discountFactor * agent.qTable[nextState[1], nextState[2], nextNextAction])

            agent.qTable[env.currentState[1], env.currentState[2], nextAction] = newQValue

            env = @set env.currentState = nextState

            agent = @set agent.env = env

            nextAction = nextNextAction

        end

        if (episode <= EPSILON_DECAY_END_EP && episode >= EPSILON_DECAY_START_EP)
            agent = @set agent.epsilon = agent.epsilon - EPSILON_DECAY
        end

        if (episode % visualizeEps == 0)
            if (gifFile === nothing)
                sleep(1);
            end;
        end;

        if (visualizeEps > episodes && episode % SAMPLING_EP == 0)
            push!(y, accumulatedReward / episode);
            print("Trained $episode/$episodes episodes\r");
            flush(stdout);
        end

    end;

    print("\n");

    if (gifFile !== nothing)
        print("Writing GIF\n");
        Environments.Plots.gif(anim, gifFile, fps = 10);
    else
        if (visualizeEps <= episodes)
            print("Press any key to continue...\n");
            readline();
        end
    end;
    return y;
end;

function evaluate(agent, times = 100000; samplingRate = 0.001)
    result = [];
    SAMPLING_TIME = trunc(Int, 1 / samplingRate);
    accumulatedReward = 0;
    print("\n");
    for time in 1:times
        agent = @set agent.env = Environments.reset(agent.env);
        accumulatedReward += Environments.visualize(agent.env, agent.env.currentState, agent.qTable, show = false, render = false);
        if (time % SAMPLING_TIME == 0)
            push!(result, accumulatedReward / time);
            print("Evaluated $time/$times times\r");
            flush(stdout);
        end
    end
    return result;
end;

end;


text:
"""Generate documents with LDA and HLDA
"""

module LDA
using Distributions
using StatsBase
using DirichletProcess

vocabulary = ["see", "spot", "run"]
n_terms = length(vocabulary)
mean_document_length = 5
term_dirichlet_parameter = 1
topic_dirichlet_parameter = 1

mutable struct Document
    topics::Array{Vector{Float64}}
    words::Array{String}
    len::Int
    Document() = new([], [], 0)
end

function add_word!(doc::Document, word, topic)
    push!(doc.topics, topic)
    push!(doc.words, word)
    doc.len += 1
    return nothing
end

function print_topic_counts(documents)
    for (i, doc) in enumerate(documents)
        println("Doc: $(i)")
        for (topic, count) in countmap(doc.topics)
            println("      count: $(count) | topic: $(round.(topic, digits=2))")
        end
    end
end

function finite(n_documents)
    n_topics = 2

    term_dirichlet_vector = ones(n_terms) * term_dirichlet_parameter
    term_distributions = [rand(Dirichlet(term_dirichlet_vector)) for _ in 1:n_topics]
    base_distribution = () -> rand(term_distributions)

    documents = [Document() for _ in 1:n_documents]

    for doc in documents
        topic_distribution = DP(base_distribution, topic_dirichlet_parameter)
        document_length = rand(Poisson(mean_document_length))
        for _ in 1:document_length
            p = dp_sample!(topic_distribution)
            word = vocabulary[rand(Categorical(p))]
            add_word!(doc, word, p)
        end
        println(doc.words)
    end

    print_topic_counts(documents)

    return nothing
end

function infinite(n_documents)
    base_dp_parameter = 10
    base_distribution = () -> rand(Dirichlet(term_dirichlet_parameter * ones(n_terms)))
    base_dp = DP(base_distribution, base_dp_parameter)

    documents = [Document() for _ in 1:n_documents]

    for doc in documents
        topic_distribution = DP(() -> dp_sample!(base_dp), topic_dirichlet_parameter)
        document_length = rand(Poisson(mean_document_length))
        for _ in 1:document_length
            p = dp_sample!(topic_distribution)
            word = vocabulary[rand(Categorical(p))]
            add_word!(doc, word, p)
        end
        println(doc.words)
    end

    print_topic_counts(documents)

    return nothing
end

end

ext: jl


text:
module IrisSpectral

using IrisBase,
ArnoldiMethod,
ArnoldiMethodWrapper,
LinearAlgebra,
SparseArrays,
NonlinearEigenproblems,
Distributed,
LinearAlgebra,
Random,
ProgressMeter

export resonance_eigenproblem,
cf_eigenproblem,
resonance_nonlinear_eigenproblem,
eig_kl,
eig_cf,
eig_knl

include("linear.jl")
include("nonlinear.jl")
include("contour_beyn_progress_bar.jl")

"""
    eig_kl(sim, k[, ka=0, kb=0]; display=false, kwargs...) -> k,ψ
"""
function eig_kl(sim::Simulation, k::Number, ka::Number=0, kb::Number=0; display::Bool=false, kwargs...)
    A, B, σ = resonance_eigenproblem(sim, k, ka, kb)
    λ, ψ, history = partialeigen(A, B, σ; diag_inv_B=true, kwargs...)
    @assert history.converged history
    display ? println(history) : nothing

    k = sqrt.(λ)
    # Normalize wavefunctions according to (ψ₁,ψ₂)=δ₁₂, which requires transformed ε or F
    normalize!(sim,ψ,B)
    return k, ψ
end

"""
    eig_cf(sim, k[, ka=0, kb=0]; display=false, kwargs...) -> k,ψ
"""
function eig_cf(sim::Simulation, k::Number, ka::Number=0, kb::Number=0; η::Number=0, display::Bool=false, kwargs...)
    A, B, σ = cf_eigenproblem(sim, k, ka, kb)

    decomp, history = partialschur(A, B, η; diag_inv_B=true, kwargs...)
    @assert history.converged history
    display ? println(history) : nothing

    # Normalize wavefunctions according to (ψ₁,ψ₂)=δ₁₂, which requires transformed ε or F
    normalize!(sim,decomp.Q,B)
    return decomp.eigenvalues, Array(decomp.Q)
end

"""
    eig_knl(sim, k, ka=0, kb=0; method=contour_beyn, nk=3, display=false, quad_n=100, kwargs...) -> k,ψ
"""
function eig_knl(sim::Simulation, k::Number, ka::Number=0, kb::Number=0;
        quad_n::Int=100,
        display::Bool=false,
        method::Function=contour_beyn,
        nev::Int=3,
        quad_method=nprocs()>1 ? :ptrapz_parallel : :ptrapz,
        kwargs...
        )

    nep = resonance_nonlinear_eigenproblem(sim, k, ka, kb; check_consistency=false)
    displaylevel = display ? 1 : 0
    if display && method==contour_beyn
        k, ψ = contour_beyn(nep, true; N=quad_n, σ=k, quad_method=quad_method, neigs=nev, kwargs...)
    else
        k, ψ = method(nep; N=quad_n, σ=k, displaylevel=displaylevel, neigs=nev, quad_method=quad_method, kwargs...)
    end
    return k, ψ
end

function LinearAlgebra.normalize!(sim::Simulation,ψ,B)
    dx = sim.dis.dx
    for i ∈ 1:size(ψ,2)
        𝒩² = sum((ψ[:,i].^2).*diag(B))*(isinf(dx[1]) ? 1 : dx[1])*(isinf(dx[2]) ? 1 : dx[2])
        ψ[:,i] /= sqrt(𝒩²)*exp(complex(0,angle(ψ[end÷2-1,i])))
    end
    return nothing
end

end # module


text:
using JSON

"""
export_json(fd::FluxiData)
  returns String

creates the json representation of the fluxi data
"""
function export_json( fd::FluxiData )
    rxns = map( x::Flux -> export_json(x) , map( xs -> fd.net.fluxes_db[xs] ,fd.net.fluxes ) )
    species = map( x::Species -> export_json(x) , map( xs -> fd.net.species_db[xs] , fd.net.species ) )

    substrate_configs = Dict{String,Any}()
    for sci in collect(keys(fd.substrate_configs))
        substrate_configs[sci] = export_json(fd.substrate_configs[sci])
    end

    measurements = Dict{String,Any}()
    for msi in collect(keys(fd.measurements))
        measurements[msi] = export_json(fd.measurements[msi])
    end

    linconstraints = Dict{String,Any}()
    for lci in keys(fd.linear_constraints)
        linconstraints[lci] = export_json(fd.linear_constraints[lci])
    end

    json_buffer::IOBuffer = IOBuffer()

    fd = Dict{String,Any}()
    fd["reactions"] = rxns
    fd["metabolites"] = species
    fd["substrate_configs"] = substrate_configs
    fd["measurements"] = measurements
    fd["linear_constraints"] = linconstraints

    JSON.print(json_buffer, fd )

    return String(take!(json_buffer))
end

function export_json( f::Flux )
    fd = Dict{String,Any}()

    fd["name"] = f.id
    fd["id"]   = f.id

    if(!is_sx_flux(f))
        fd["lower_bound"] = 0.0
        fd["upper_bound"] = 1000.0

        fd["metabolites"] = add( f.P , mult( f.E , -1. ) )
        fd["metabolites_E"] = f.E
        fd["metabolites_P"] = f.P

        #fd["transition"] = fd.Transition
        # add atom transition data:
        fd["atomtransition"] = write_atom_transition(f)
    else
        fd["lower_bound"] = -1000.0
        fd["upper_bound"] =  1000.0

        if(haskey(f.E,"*"))
            fd["metabolites"] = f.P
            fd["metabolites_E"] = Dict{String,Any}()
            fd["metabolites_P"] = f.P
        else
            if(!haskey(f.P,"*"))
                error("Strange flux encountered: $(f)")
            end
            fd["metabolites"] = f.E
            fd["metabolites_E"] = Dict{String,Any}()
            fd["metabolites_P"] = f.E
        end
    end
    return fd
end

function export_json( s::Species )
    js = Dict{String,Any}()

    js["name"] = s.id
    js["id"]   = s.id
    js["notes"] = Dict{String,Any}( "CARBONS" => s.na )

    return js
end

"""
import_json_to_fluxi(str_json::String)
  returns FluxiData

creates the fluxi data object from the fluxi json representation
"""
function import_json_to_fluxi( str_json::String )
    imp = import_json(str_json)
    linear_constraints_sets = imp[7]
    linear_constraints_sets["rxn_bounds"] = imp[2]
    return FluxiData(imp[1],imp[5],imp[6],imp[7])
end

function import_json( str_json::String )
    json_d = JSON.parse(str_json)
    import_json(json_d)
end

function import_json( json_d::Dict{String,Any} )
    import_s = import_json_species(json_d["metabolites"])
    import_f = import_json_fluxes(json_d["reactions"])

    import_measurements      = import_json_measurements(json_d["measurements"])
    import_substrate_configs = import_json_substrconfig(json_d["substrate_configs"])
    import_linconstraints    = import_json_linconstraints(json_d["linear_constraints"])

    return ( Network(import_f[1],import_s[1], import_f[2],import_s[2]), import_f[4],import_s,import_f,import_measurements,import_substrate_configs,import_linconstraints)
end

function import_json_species( species_list::Vector{Any})
    map_id2speciesdata = Dict{String,Dict}()
    map_id2species     = Dict{String,Species}()
    list_species       = Vector{String}()

    for dsi in species_list
        id        = ""
        num_atoms = -1

        if( !isa(dsi,Dict) )
            error("unexpected species.. no dict?")
        end
        id = dsi["id"]
        map_id2speciesdata[id] = dsi

        if(haskey(dsi,"notes"))
            if( haskey(dsi["notes"],"CARBONS") )
                if(isa(dsi["notes"]["CARBONS"][1],String))
                    num_atoms = parse(Float64,dsi["notes"]["CARBONS"][1])
                end
                if(isa(dsi["notes"]["CARBONS"][1],Real))
                    num_atoms = Int(dsi["notes"]["CARBONS"][1])
                end
            end
        end
        #println(num_atoms)
        s1 = Species(id,Integer(round(num_atoms)))
        push!(list_species,id)
        map_id2species[id] = s1
    end

    return (list_species,map_id2species,map_id2speciesdata)
end

function import_json_fluxes( fluxes_list::Vector{Any})
    map_id2fluxdata   = Dict{String,Dict}()
    map_id2fluxes     = Dict{String,Flux}()
    list_fluxes       = Vector{String}()

    lin_constraints   = Vector{LinConstr}()

    for dfi in fluxes_list
        stoich::Dict{String,Float64} = Dict{String,Float64}()
        for ki in keys(dfi["metabolites"])
            stoich[ki] = dfi["metabolites"][ki]
        end
        if( dfi["lower_bound"] < 0 )
            #split into two..
            id = dfi["id"]
            id_fw = string(id,"_fw")
            id_bw = string(id,"_bw")

            fi_fw = Flux(id, stoich )
            fi_bw = Flux(id, mult(stoich,-1.) )

            map_id2fluxdata[id_fw] = dfi
            map_id2fluxdata[id_bw] = dfi
            push!(list_fluxes,id_fw)
            push!(list_fluxes,id_bw)
            map_id2fluxes[id_fw] = fi_fw
            map_id2fluxes[id_bw] = fi_bw

            push!(lin_constraints,parse_linconstr("$(id_fw) < $(dfi["upper_bound"])") )
            push!(lin_constraints,parse_linconstr("$(id_bw) < $(-dfi["lower_bound"])") )
        else
            #one flux
            id = dfi["id"]
            map_id2fluxdata[id] = dfi
            fi = Flux(id, stoich )
            push!(list_fluxes,id)
            map_id2fluxes[id] = fi

            push!(lin_constraints,parse_linconstr("$(id) < $(dfi["upper_bound"])") )
            if(dfi["lower_bound"]>0)
                push!(lin_constraints,parse_linconstr("$(id) > $(dfi["lower_bound"])") )
            end
        end
    end
    return (list_fluxes,map_id2fluxes,map_id2fluxdata,lin_constraints)
end

function export_json( substr_config:: Vector{ Dict{MetabolicFluxAnalysis.ConfiguredSEMU,Float64} } )
    #fd = Dict{String,Any}()
    substrates = Vector{Dict{String,Any}}()
    for conf_s in substr_config
        d_s = Dict{String,Any}()
        d_s["metabolite"] = collect( collect(keys(conf_s))[1].semu)[1].species
        d_s["pos"]        = get_positions_of_semu( collect(keys(conf_s))[1].semu )
        list_mix          = Vector{Dict{String,Any}}()
        for csr in keys( conf_s )
            d_mi = Dict{String,Any}()
            d_mi["conf"]  = csr.config
            d_mi["ratio"] = conf_s[csr]
            push!( list_mix , d_mi )
        end
        d_s["mix"] = list_mix
        push!( substrates , d_s)
    end
    return substrates
end

function array_to_intarray( va::Vector{Any})
    vi = Vector{Integer}(length(va))
    for zi in 1:length(va)
        vi[zi] = Integer(va[zi])
    end
    return vi
end

#function import_json_substrconfig( substr_conf::Dict{String,Vector{Dict{String,Any}}} )
function import_json_substrconfig( d_substr_conf::Dict{String,Any} )
    substrate_configs = Dict{String,Vector{Dict{MetabolicFluxAnalysis.ConfiguredSEMU,Float64}}}()
    for key in collect(keys(d_substr_conf))
        substr_i = d_substr_conf[key]
        substrate_config_i = Vector{Dict{MetabolicFluxAnalysis.ConfiguredSEMU,Float64}}()
        for si in substr_i
            semu_i = constructSpeciesEMU( si["metabolite"] , array_to_intarray( si["pos"] ) )
            configured_substrate_i = Dict{ConfiguredSEMU,Float64}()
            for dsc in si["mix"]
                csemu_i = ConfiguredSEMU(semu_i, array_to_intarray( dsc["conf"] ) )
                configured_substrate_i[csemu_i] = dsc["ratio"]
            end
            push!(substrate_config_i,configured_substrate_i)
        end
        substrate_configs[key] = substrate_config_i
    end
    return substrate_configs
end

function export_json( measurements_set::Vector{MetabolicFluxAnalysis.MassDistributionMeasurement} )
    measurements = Vector{Dict{String,Any}}()
    for mi in measurements_set
        dm = Dict{String,Any}()
        dm["metabolite"] = collect( mi.sass)[1].species
        dm["pos"]        = get_positions_of_semu( mi.sass )

        dm_data          = Vector{Dict{String,Any}}()
        for msi in keys(mi.data)
            d_data_i = Dict{String,Any}()
            d_data_i["mass"] = msi
            d_data_i["data"] = mi.data[msi]
            push!(dm_data,d_data_i)
        end
        dm["data"] = dm_data
        push!(measurements,dm)
    end
    return measurements
end

function import_json_measurements(d_measurements::Dict{String,Any})
    measurements = Dict{String,Vector{MetabolicFluxAnalysis.MassDistributionMeasurement}}()
    for key in collect(keys(d_measurements))
        meas_i = d_measurements[key]
        measurements_i = Vector{MetabolicFluxAnalysis.MassDistributionMeasurement}()
        for msi in meas_i
            semu_i = constructSpeciesEMU( msi["metabolite"] , array_to_intarray( msi["pos"] ) )
            mass_data_i = Dict{Int64,Vector{Float64}}()
            for mass_di in msi["data"]
                mass_data_i[Int64(mass_di["mass"])] = mass_di["data"]
            end
            mdi_i = MassDistributionMeasurement(semu_i,mass_data_i)
            push!(measurements_i,mdi_i)
        end
        measurements[key] = measurements_i
    end
    return measurements
end

function export_json( lcs::Vector{MetabolicFluxAnalysis.LinConstr} )
    vlcs = Vector{Dict{String,Any}}()
    for lci in lcs
        d_lci = Dict{String,Any}()
        d_lci["lhs"] = lci.lhs
        d_lci["rhs"] = lci.rhs
        d_lci["op"]  = lci.constr
        push!(vlcs,d_lci)
    end
    return vlcs
end

function export_json(lci::LC)
    return lci
end

function import_json_lc(d_lc::Dict{String,Any})
    lci = LC{String}()
    for ki in keys(d_lc)
        lci[ki] = Float64(d_lc[ki])
    end
    return lci
end

function import_json_linconstraints(d_lcs::Dict{String,Any})
    constraints = Dict{String,Vector{MetabolicFluxAnalysis.LinConstr}}()
    for key in collect(keys(d_lcs))
        lcs_i = d_lcs[key]
        linconstr_i = Vector{MetabolicFluxAnalysis.LinConstr}()

        for d_lci in lcs_i
            dlhs = import_json_lc(d_lci["lhs"])
            drhs = import_json_lc(d_lci["rhs"])
            d_op = 'x'
            if(isa(d_lci["op"],String)); d_op = d_lci["op"][1]; end
            if(isa(d_lci["op"],Char)); d_op   = d_lci["op"]; end
            lc_new = LinConstr( d_op , dlhs , drhs )
            push!(linconstr_i,lc_new)
        end
        constraints[key] = linconstr_i
    end
    return constraints
end


text:
function test_susceptibility(scenario)
    @unpack tx, bfield, species, ground = scenario

    M1 = LMP.susceptibility(70e3, tx.frequency, bfield, species)
    M2 = LMP.susceptibility(70e3, tx.frequency, bfield, species; params=LMPParams())
    M3 = LMP.susceptibility(70e3, tx.frequency, bfield, species; params=LMPParams(earthradius=6350e3))
    @test M1 == M2
    @test !(M2 ≈ M3)

    waveguide = HomogeneousWaveguide(bfield, species, ground)
    modeequation = PhysicalModeEquation(tx.frequency, waveguide)

    M4 = LMP.susceptibility(70e3, tx.frequency, waveguide)
    M5 = LMP.susceptibility(70e3, tx.frequency, waveguide; params=LMPParams())
    M6 = LMP.susceptibility(70e3, tx.frequency, waveguide; params=LMPParams(earthradius=6350e3))

    M7 = LMP.susceptibility(70e3, modeequation)
    M8 = LMP.susceptibility(70e3, modeequation; params=LMPParams())
    M9 = LMP.susceptibility(70e3, modeequation; params=LMPParams(earthradius=6350e3))

    @test M4 == M5 == M1
    @test M6 == M3

    @test M7 == M8 == M1
    @test M9 == M3

    @inferred LMP.susceptibility(70e3, tx.frequency, bfield, species)
end

function test_spline(scenario)
    @unpack tx, bfield, species = scenario

    itp = LMP.susceptibilityspline(tx.frequency, bfield, species)

    zs = LMP.BOTTOMHEIGHT:1:LMPParams().topheight
    Ms = LMP.susceptibility.(zs, (tx.frequency,), (bfield,), (species,))
    Mitp = itp.(zs)

    function matrix(M)
        mat = Matrix{eltype(M[1])}(undef, length(M), 9)
        for i in eachindex(M)
            mat[i,:] = M[i]
        end
        return mat
    end

    @test matrix(Mitp) ≈ matrix(Ms) rtol=1e-5

    # plot(matrix(real(Ms)), zs/1000, legend=false)
    # plot!(matrix(real(Mitp)), zs/1000, linestyle=:dash)

    # Mdiff = matrix(real(Ms)) .- matrix(real(Mitp))
    # scaleddiff = Mdiff./matrix(real(Ms))
    # plot(Mdiff, zs/1000, legend=false, ylims=(0,110))

    # Make sure `params.susceptibilitysplinestep` is being used
    itp2 = LMP.susceptibilityspline(tx.frequency, bfield, species; params=LMPParams(susceptibilitysplinestep=10.0))
    itp3 = LMP.susceptibilityspline(tx.frequency, bfield, species; params=LMPParams(susceptibilitysplinestep=20.0))
    @test itp == itp2
    @test itp != itp3
end

@testset "magnetoionic.jl" begin
    @info "Testing magnetoionic"

    for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario,
            multiplespecies_scenario)
        test_susceptibility(scn)
        test_spline(scn)
    end
end


text:
module SplitExplicit01

export OceanDGModel,
    OceanModel,
    Continuity3dModel,
    HorizontalModel,
    BarotropicModel,
    AbstractOceanProblem

#using Printf
using StaticArrays
using LinearAlgebra: I, dot, Diagonal

using ...VariableTemplates
using ...MPIStateArrays
using ...DGMethods: init_ode_state, basic_grid_info
using ...Mesh.Filters: CutoffFilter, apply!, ExponentialFilter
using ...Mesh.Grids:
    polynomialorders,
    dimensionality,
    dofs_per_element,
    VerticalDirection,
    HorizontalDirection,
    min_node_distance

using ...BalanceLaws
#import ...BalanceLaws: nodal_update_auxiliary_state!

using ...DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder,
    NumericalFluxGradient,
    NumericalFluxSecondOrder,
    RusanovNumericalFlux,
    CentralNumericalFluxFirstOrder,
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder

using ..Ocean: AbstractOceanProblem

import ...DGMethods.NumericalFluxes:
    update_penalty!, numerical_flux_second_order!, NumericalFluxFirstOrder

import ...DGMethods: LocalGeometry, DGModel, calculate_dt

import ...BalanceLaws:
    vars_state,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    parameter_set,
    boundary_conditions,
    boundary_state!,
    update_auxiliary_state!,
    update_auxiliary_state_gradient!,
    compute_gradient_argument!,
    init_state_auxiliary!,
    init_state_prognostic!,
    compute_gradient_flux!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

import ...SystemSolvers: BatchedGeneralizedMinimalResidual, linearsolve!

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
∘(a::SVector, b::SVector) = StaticArrays.dot(a, b)

abstract type AbstractOceanModel <: BalanceLaw end

function ocean_init_aux! end
function ocean_init_state! end
function ocean_model_boundary! end
function set_fast_for_stepping! end
function initialize_fast_state! end
function initialize_adjustment! end

include("SplitExplicitLSRK2nMethod.jl")
include("SplitExplicitLSRK3nMethod.jl")
include("OceanModel.jl")
include("Continuity3dModel.jl")
include("VerticalIntegralModel.jl")
include("BarotropicModel.jl")
include("IVDCModel.jl")
include("Communication.jl")
include("OceanBoundaryConditions.jl")

end


text:
__precompile__(true)

import Base: convert, sparse, size, diag, eltype, ndims, ==, *, .*, issymmetric, A_mul_B!, length, Diagonal
export convert,
	SparseMatrix,
	GraphMatrix,
	Adjacency,
	adjacency,
	Laplacian,
	CombinatorialAdjacency,
	CombinatorialLaplacian,
	NormalizedAdjacency,
	NormalizedLaplacian,
	StochasticAdjacency,
	StochasticLaplacian,
	AveragingAdjacency,
	AveragingLaplacian,
	PunchedAdjacency,
	Noop,
	diag,
	degrees,
	symmetrize,
	prescalefactor,
	postscalefactor

const SparseMatrix{T} = SparseMatrixCSC{T,Int64}

"""
	GraphMatrix{T}

An abstract type to allow opertions on any type of graph matrix
"""
abstract type GraphMatrix{T} end

"""
	Adjacency{T}

The core Adjacency matrix structure. Keeps the vertex degrees around.
Subtypes are used to represent the different normalizations of the adjacency matrix.
Laplacian and its subtypes are used for the different Laplacian matrices.

Adjacency(lapl::Laplacian) provides a generic function for getting the
adjacency matrix of a Laplacian matrix. If your subtype of Laplacian does not provide
an field A for the Adjacency instance, then attach another method to this function to provide
an Adjacency{T} representation of the Laplacian. The Adjacency matrix here
is the final subtype that corresponds to this type of Laplacian
"""
abstract type Adjacency{T} <: GraphMatrix{T} end
abstract type Laplacian{T} <: GraphMatrix{T} end

"""
	CombinatorialAdjacency{T,S,V}

The standard adjacency matrix.
"""
struct CombinatorialAdjacency{T,S,V} <: Adjacency{T}
	A::S
	D::V
end

function CombinatorialAdjacency(A::SparseMatrix{T}) where T
	D = vec(sum(A, 1))
	return CombinatorialAdjacency{T,SparseMatrix{T},typeof(D)}(A, D)
end

@doc_str """
	NormalizedAdjacency{T}

The normalized adjacency matrix is ``\\hat{A} = D^{-1/2} A D^{-1/2}``.
If A is symmetric, then the normalized adjacency is also symmetric
with real eigenvalues bounded by [-1, 1].
"""
struct NormalizedAdjacency{T} <: Adjacency{T}
	A::CombinatorialAdjacency{T}
	scalefactor::Vector{T}
end
function NormalizedAdjacency(adjmat::CombinatorialAdjacency)
	sf = adjmat.D.^(-1 / 2)
	return NormalizedAdjacency(adjmat, sf)
end

"""
	StochasticAdjacency{T}

A transition matrix for the random walk.
"""
struct StochasticAdjacency{T} <: Adjacency{T}
	A::CombinatorialAdjacency{T}
	scalefactor::Vector{T}

end
function StochasticAdjacency(adjmat::CombinatorialAdjacency)
	sf = adjmat.D.^(-1)
	return StochasticAdjacency(adjmat, sf)
end

"""
	AveragingAdjacency{T}

The matrix whose action is to average over each neighborhood.
"""
struct AveragingAdjacency{T} <: Adjacency{T}
	A::CombinatorialAdjacency{T}
	scalefactor::Vector{T}
end
function AveragingAdjacency(adjmat::CombinatorialAdjacency)
	sf = adjmat.D.^(-1)
	return AveragingAdjacency(adjmat, sf)
end

perron(adjmat::NormalizedAdjacency) = sqrt.(adjmat.A.D) / norm(sqrt.(adjmat.A.D))

struct PunchedAdjacency{T} <: Adjacency{T}
	A::NormalizedAdjacency{T}
	perron::Vector{T}
end
function PunchedAdjacency(adjmat::CombinatorialAdjacency)
            perron = sqrt.(adjmat.D) / norm(sqrt.(adjmat.D))
            return PunchedAdjacency(NormalizedAdjacency(adjmat), perron)
end

perron(m::PunchedAdjacency) = m.perron

"""
	Noop

A type that represents no action.

### Implementation Notes
- The purpose of `Noop` is to help write more general code for the
different scaled GraphMatrix types.
"""
struct Noop end

Base.broadcast(::typeof(*), ::Noop, x) = x

Diagonal(::Noop) = Noop()

==(g::GraphMatrix, h::GraphMatrix) = typeof(g) == typeof(h) && (g.A == h.A)

postscalefactor(::Adjacency) = Noop()

postscalefactor(adjmat::NormalizedAdjacency) = adjmat.scalefactor

postscalefactor(adjmat::AveragingAdjacency) = adjmat.scalefactor

prescalefactor(::Adjacency) = Noop()

prescalefactor(adjmat::NormalizedAdjacency) = adjmat.scalefactor

prescalefactor(adjmat::StochasticAdjacency) = adjmat.scalefactor

struct CombinatorialLaplacian{T} <: Laplacian{T}
	A::CombinatorialAdjacency{T}
end

@doc_str """
	NormalizedLaplacian{T}

The normalized Laplacian is ``\\hat{L} = I - D^{-1/2} A D^{-1/2}``.
If A is symmetric, then the normalized Laplacian is also symmetric
with positive eigenvalues bounded by 2.
"""
struct NormalizedLaplacian{T} <: Laplacian{T}
	A::NormalizedAdjacency{T}
end

"""
	StochasticLaplacian{T}

Laplacian version of the StochasticAdjacency matrix.
"""
struct StochasticLaplacian{T} <: Laplacian{T}
	A::StochasticAdjacency{T}
end

"""
	AveragingLaplacian{T}

Laplacian version of the AveragingAdjacency matrix.
"""
struct AveragingLaplacian{T} <: Laplacian{T}
	A::AveragingAdjacency{T}
end

arrayfunctions = (:eltype, :length, :ndims, :size, :strides, :issymmetric)
for f in arrayfunctions
	@eval $f(a::GraphMatrix) = $f(a.A)
end

size(a::GraphMatrix, i::Integer) = size(a.A, i)
issymmetric(::StochasticAdjacency) = false
issymmetric(::AveragingAdjacency) = false

"""
	degrees(adjmat)

Return the degrees of a graph represented by the [CombinatorialAdjacency](@ref) `adjmat`.
"""
degrees(adjmat::CombinatorialAdjacency) = adjmat.D

"""
	degrees(graphmx)

Return the degrees of a graph represented by the graph matrix `graphmx`.
"""
degrees(mat::GraphMatrix) = degrees(adjacency(mat))

adjacency(lapl::Laplacian) = lapl.A
adjacency(lapl::GraphMatrix) = lapl.A

convert(::Type{CombinatorialAdjacency}, adjmat::Adjacency) = adjmat.A
convert(::Type{CombinatorialAdjacency}, adjmat::CombinatorialAdjacency) = adjmat

function sparse(lapl::M) where M<:Laplacian
	adjmat = adjacency(lapl)
	A = sparse(adjmat)
	L = spdiagm(diag(lapl)) - A
	return L
end

function sparse(adjmat::Adjacency)
    A = sparse(adjmat.A)
    return Diagonal(prescalefactor(adjmat)) * (A * Diagonal(postscalefactor(adjmat)))
end

function convert(::Type{SparseMatrix{T}}, adjmat::Adjacency{T}) where T
    A = sparse(adjmat.A)
    return Diagonal(prescalefactor(adjmat)) * (A * Diagonal(postscalefactor(adjmat)))
end

function convert(::Type{SparseMatrix{T}}, lapl::Laplacian{T}) where T
	adjmat = adjacency(lapl)
	A = convert(SparseMatrix{T}, adjmat)
	L = spdiagm(diag(lapl)) - A
	return L
end

diag(lapl::CombinatorialLaplacian) = lapl.A.D
diag(lapl::Laplacian) = ones(size(lapl)[2])

*(x::AbstractArray, ::Noop) = x
*(::Noop, x) = x
*{T<:Number}(adjmat::Adjacency{T}, x::AbstractVector{T}) =
	postscalefactor(adjmat) .* (adjmat.A * (prescalefactor(adjmat) .* x))

*(adjmat::CombinatorialAdjacency{T}, x::AbstractVector{T}) where T<:Number =
	adjmat.A * x

*(lapl::Laplacian{T}, x::AbstractVector{T}) where T<:Number =
	(diag(lapl) .* x) - (adjacency(lapl) * x)

function *(adjmat::PunchedAdjacency{T}, x::AbstractVector{T}) where T<:Number
    y = adjmat.A * x
    return y - dot(adjmat.perron, y) * adjmat.perron
end

function A_mul_B!(Y, A::Adjacency, B)
    # we need to do 3 matrix products
    # Y and B can't overlap in any one call to A_mul_B!
    # The last call to A_mul_B! must be (Y, postscalefactor, tmp)
    # so we need to write to tmp in the second step  must be (tmp, A.A, Y)
    # and the first step (Y, prescalefactor, B)
    tmp1 = Diagonal(prescalefactor(A)) * B
    tmp = similar(Y)
    A_mul_B!(tmp, A.A, tmp1)
    return A_mul_B!(Y, Diagonal(postscalefactor(A)), tmp)
end

A_mul_B!(Y, A::CombinatorialAdjacency, B) = A_mul_B!(Y, A.A, B)

# You can compute the StochasticAdjacency product without allocating a similar of Y.
# This is true for all Adjacency where the postscalefactor is a Noop
# at time of writing this is just StochasticAdjacency and CombinatorialAdjacency
function A_mul_B!(Y, A::StochasticAdjacency, B)
    tmp = Diagonal(prescalefactor(A)) * B
    A_mul_B!(Y, A.A, tmp)
    return Y
end

function A_mul_B!(Y, adjmat::PunchedAdjacency, x)
    y = adjmat.A * x
    Y[:] = y - dot(adjmat.perron, y) * adjmat.perron
    return Y
end

function A_mul_B!(Y, lapl::Laplacian, B)
	A_mul_B!(Y, lapl.A, B)
	z = diag(lapl) .* B
	Y[:] = z - Y[:]
	return Y
end

"""
	symmetrize(A::SparseMatrix, which=:or)

Return a symmetric version of graph (represented by sparse matrix `A`) as a sparse matrix.
`which` may be one of `:triu`, `:tril`, `:sum`, or `:or`. Use `:sum` for weighted graphs.
"""
function symmetrize(A::SparseMatrix, which=:or)
	if which == :or
		M = A + A'
		M.nzval[M.nzval .== 2] = 1
		return M
    end
	T = A
	if which == :triu
		T = triu(A)
	elseif which == :tril
	    T = tril(A)
    elseif which == :sum
	    T = A
    else
        error("$which is not a supported method of symmetrizing a matrix")
    end
	M = T + T'
	return M
end

"""
	symmetrize(adjmat, which=:or)

Return a symmetric version of graph (represented by [`CombinatorialAdjacency`](@ref) `adjmat`)
as a [`CombinatorialAdjacency`](@ref). `which` may be one of `:triu`, `:tril`, `:sum`, or `:or`.
Use `:sum` for weighted graphs.

### Implementation Notes
Only works on [Adjacency](@ref) because the normalizations don't commute with symmetrization.
"""
symmetrize(adjmat::CombinatorialAdjacency, which=:or) =
	CombinatorialAdjacency(symmetrize(adjmat.A, which))

# per #564
@deprecate A_mul_B!(Y, A::Noop, B) None
@deprecate convert(::Type{Adjacency}, lapl::Laplacian) None
@deprecate convert(::Type{SparseMatrix}, adjmat::CombinatorialAdjacency) sparse(adjmat)

"""
	LinAlg

A package for using the type system to check types of graph matrices.
"""
LinAlg


text:
    function contractmodes(tensor1::Array{Float64},tensor2::Array{Float64},ind_contract::Matrix{Int64},ind_rem::Matrix{Int64})
        # function that contracts modes specified in ind_contract, combines indices specified in ind_comn and leaves indices free specified in ind_free
        # contraction
        matrix1  = unfold(tensor1,ind_contract[:,1],"right");
        matrix2  = unfold(tensor2,ind_contract[:,2],"left");
        matrix   = matrix1*matrix2;
        sizes_1  = deleteat!(collect(size(tensor1)),sort(ind_contract[:,1]));
        sizes_2  = deleteat!(collect(size(tensor2)),sort(ind_contract[:,2]));
        newsizes = vcat(sizes_1,sizes_2);
        tensor   = reshape(matrix,newsizes...);

        # combining specified modes and freeing specified modes
        dimrem_1 = ndims(tensor1) - size(ind_contract,1);
        dimrem_2 = ndims(tensor2) - size(ind_contract,1)
        newperm  = Vector{Int64}(undef,dimrem_1+dimrem_2);
        for i = 1:size(ind_contract,1)
            for j = 1:size(ind_rem,1)
                if ind_contract[i,1] < ind_rem[j,1]
                    ind_rem[j,1] = ind_rem[j,1]-1;
                end
                if ind_contract[i,2] < ind_rem[j,2]
                    ind_rem[j,2] = ind_rem[j,2]-1;
                end
            end
        end
        k = 1;
        for i = 1:size(ind_rem,1)
            if ind_rem[i,1] > 0
                newperm[k] = ind_rem[i,1];
                k = k+1;
            end
            if ind_rem[i,2] > 0
                newperm[k] = ind_rem[i,2]+dimrem_1;
                k = k+1;
            end
        end
        tensor   = permutedims(tensor,Tuple(newperm));

        newsizes = Int64[];
        for i = 1:size(ind_rem,1)
            if ind_rem[i,1] > 0
                size_1 = sizes_1[ind_rem[i,1]];
            else
                size_1 = 1;
            end
            if ind_rem[i,2] > 0
                size_2 = sizes_2[ind_rem[i,2]];
            else
                size_2 = 1;
            end
            push!(newsizes,size_1*size_2);
        end
        reshape(tensor,newsizes...);
    end

    function contractmodes(tensor1::Array{Float64},tensor2::Array{Float64},ind_contract::Matrix{Int64})
        # contract specified indices and free the remaining indices
        matrix1  = unfold(tensor1,ind_contract[:,1],"right")
        matrix2  = unfold(tensor2,ind_contract[:,2],"left")
        matrix   = matrix1*matrix2
        sizes_1  = deleteat!(collect(size(tensor1)),sort(ind_contract[:,1]));
        sizes_2  = deleteat!(collect(size(tensor2)),sort(ind_contract[:,2]));
        newsizes = vcat(sizes_1,sizes_2);
        reshape(matrix,newsizes...);
    end

    function contractcores(core1::Array{Float64,4},core2::Array{Float64,3})
        contractmodes(core1,core2,[3 2],[1 1; 2 0; 4 3]);
    end

    function contractcores(core1::Array{Float64,4},core2::Array{Float64,4})
        contractmodes(core1,core2,[3 2],[1 1; 2 0; 0 3; 4 4]);
    end

    function contractcores(core1::Array{Float64,3},core2::Array{Float64,3})
        contractmodes(core1,core2,[2 2],[1 1; 3 3]);
    end

    function contractcores(core1::Array{Float64,3},core2::Array{Float64,4})
        contractmodes(core1,core2,[2 2],[1 1; 0 3; 3 4]);
    end

    function nmodeproduct(matrix::Matrix{Float64},tensor::Array{Float64},n::Int64)
    #        ___         ___
    #    ___/ A \_______/ X \_________________________
    #     J \___/  I_n  \___/ I_1...I_n-1 I_n+1...I_N
    #
    #              ____
    #             /XtnA\
    #    =====>   \____/
    #             /||||\
    #           I_1  J  I_N
    # Kolda, Bader: Tensor Decomposition and Applications, 2009. Page 460

        tensor = contractmodes(matrix,tensor,[2 n]); # contract mode n and leave all other modes free
        permutedims(tensor,[collect(2:n)..., 1, collect(n+1:ndims(tensor))...])
    end

    function vectorbymatrix(A::MPT{4},b::MPT{3})
        d = order(A)
        vec = contractcores(b[1],A[1])[1,:,:]
        for i = 2:d-1
            vec = vec * contractcores(b[i],A[i])[:,1,:]
        end
        return vec = vec * contractcores(b[d],A[d])[:,:,1]
    end

    function matrixbyvector(A::MPT{4},b::MPT{3})
        d = order(A)
        vec = contractcores(b[1],A[1])[1,:,:]
        for i = 2:d-1
            vec = vec * contractcores(b[i],A[i])[:,1,:]
        end
        return vec = (vec * contractcores(b[d],A[d])[:,:,1])'
    end

    # computation of vector represented by an MPS
    function mps2vec(mps::MPT{3})
        tensor = unfold(mps[1],[3],"right");
        for i = 2:order(mps)
            tensor = tensor*unfold(mps[i],[1],"left");
            tensor = reshape(tensor, (Int(length(tensor)/size(mps[i],3)), size(mps[i],3)));
        end
        vector = tensor;
        return vector;
    end

    # computation of matrix represented by an MPO
    function mpo2mat(mpo::MPT{4})
        tensor = unfold(mpo[1],[4],"right");
        sizes = size(mpo);
        newsizes = [sizes[1]...];
        matsize1 = sizes[1][1];
        matsize2 = sizes[1][2];
        for i = 2:order(mpo)
            tensor = tensor*unfold(mpo[i],[1],"left");
            tensor = reshape(tensor, (Int(length(tensor)/size(mpo[i],4)), size(mpo[i],4)));
            newsizes = (newsizes..., sizes[i]...);
            matsize1 = matsize1*sizes[i][1];
            matsize2 = matsize2*sizes[i][2];
        end
        tensor = reshape(tensor,newsizes);
        tensor = permutedims(tensor,[collect(1:2:length(newsizes)-1)..., collect(2:2:length(newsizes))...]);
        matrix = reshape(tensor,(matsize1,matsize2));
        return matrix
    end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
function make_unique(x::AbstractVector{K}, y::AbstractVector{K}) where {K}
    if is_dynamic(x)
        return make_unique!(copy(x), y)
    else
        # TODO Better method of dealing with non dynamic vectors in make_unique
        return convert(typeof(x), make_unique!(Vector(x), y))
    end
end

"""
    make_unique!(x, y) -> x
"""
make_unique!(x, y) = make_unique!(promote(x, y)...)

function make_unique!(
    x::AbstractVector{T},
    y::AbstractVector{T},
    offset::T=unique_offset(x)
   ) where {T}
    for y_i in y
        if y_i in x
            push!(x, add_offset(x, y_i, offset))
        else
            push!(x, y_i)
        end
    end
    return x
end

function unique_offset(x::AbstractVector{T}) where {T<:Number}
    return T(round(maximum(names), RoundUp, sigdigits=1))
end
unique_offset(x::AbstractVector{Symbol}) = :_2
unique_offset(x::AbstractVector{AbstractString}) = "_2"

function add_offset(x, v::AbstractString, f::AbstractString)
    newval = v * f
    if newval in x
        return add_offset(x, newval, f)
    else
        return newval
    end
end
function add_offset(x, v::Number, f::Number)
    newval = v + f
    if newval in x
        return add_offset(x, newval, f)
    else
        return newval
    end
end
function add_offset(x, v::Symbol, f::Symbol)
    newval = Symbol(v, f)
    if newval in x
        return add_offset(x, newval, f)
    else
        return newval
    end
end


text:
using LevelSetClosestPoints
using Test

@testset "LevelSetClosestPoints.jl" begin
    # Write your tests here.
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using GLMakie, DynamicalSystems, InteractiveDynamics
using DataStructures: CircularBuffer

# Basic figure layouting
fig = Figure(resolution = (1600, 500), figure_padding = 50, dpi = 600)
display(fig)
axs = [fig[i, 2] = Axis(fig) for i in 1:3]
ax_orig = fig[:, 1] = Axis3(fig)
ax_embe = fig[:, 3] = Axis3(fig)
# Link x axis of middle plots
for i in 1:2
    linkxaxes!(axs[3], axs[i])
    hidexdecorations!(axs[i]; grid = false)
end
# Set up a nice projection for Lorenz system
ax_orig.azimuth = 5.855
ax_embe.azimuth = 5.855
ax_orig.elevation = 0.46
ax_embe.elevation = 0.46
ax_orig.zlabeloffset = 20

# Generate a very large timeseries. The animation will run through the timeseries
Δt = 0.01
tail = 800 # in units of Δt
N = 5 # how many tails we will have
τ = 20 # delay time
ds = Systems.lorenz()
traj_orig = trajectory(ds, N*tail*Δt; dt = Δt, Ttr = 1000)
x = traj_orig[:, 1]
t = 0:Δt:N*tail*Δt # time axis
traj_embed = embed(x, 3, τ)

# Set up observables for original
traj_orig_obs = Observable(CircularBuffer{Point3f0}(tail))
append!(traj_orig_obs[], traj_orig[1:tail])
c0 = to_color(COLORSCHEME[1]) # color for 3D attractors
traj_colors = [RGBAf0(c0.r, c0.g, c0.b, i/tail) for i in 1:tail]
orig_endpoint = Observable(Point3f0(traj_orig[tail]...))

# Setup observables for middle timeseries and for embedded plot
cs = [to_color(COLORS[i+1]) for i in 1:3]
xs = [@view(x[(1 + (i-1)*τ):(end - (3-i)*τ)]) for i in 1:3]
t = t[1:length(xs[1])]
xs_points = [Observable(Point2f0(0.0, 0.0)) for i in 1:3]
traj_embed_obs = Observable(CircularBuffer{Point3f0}(tail))
append!(traj_embed_obs[], traj_embed[1:tail])
embed_endpoint = Observable(Point3f0(0,0,0))

# Plot initial trajectory
proj_line = lines!(ax_orig, xprojection; color = cs[1], linewidth = 1)
proj_line.visible = false
lines!(ax_orig, traj_orig_obs; color = traj_colors, linewidth = 2, transparency=true)
scatter!(ax_orig, orig_endpoint; markersize = 4000, color = c0)
# PLot delayed timeseries (animation will simply change time axis range)
for i in 1:3
    lines!(axs[i], t, xs[i]; color = cs[i], linewidth = 2)
    scatter!(axs[i], xs_points[i]; color = cs[i])
end
xlims!(axs[3], (0, tail*Δt))
# Plot reconstructed attractor
lines!(ax_embe, traj_embed_obs; color = traj_colors, linewidth = 2, transparency=true)
scatter!(ax_embe, embed_endpoint; markersize = 4000, color = c0)

# Update/animation function
function update_to_time!(j)
    append!(traj_orig_obs[], traj_orig[j:(tail+j-1)])
    append!(traj_embed_obs[], traj_embed[j:(tail+j-1)])
    # these two lines are needed to actually trigger update.
    traj_orig_obs[] = traj_orig_obs[]
    traj_embed_obs[] = traj_embed_obs[]
    orig_endpoint[] = traj_orig_obs[][end]
    embed_endpoint[] = traj_embed_obs[][end]
    # update points of timeseries
    for i in 1:3
        offset = tail # - (i-1)*τ
        xs_points[i][] = Point2f0((j-1 + offset)*Δt, xs[i][j + offset])
    end
    xlims!(axs[3], ((j-1 + tail÷2)*Δt, (j-1+tail + tail÷2)*Δt))
end
update_to_time!(1)

# Final labelling
ax_embe.xlabel = "x(t)"
ax_embe.ylabel = "x(t + τ)"
ax_embe.zlabel = "x(t + 2τ)"
axs[1].ylabel = "x(t)"
axs[2].ylabel = "x(t + τ)"
axs[3].ylabel = "x(t + 2τ)"
ax_orig.title = "original set"
axs[1].title = "measured timeseries"
ax_embe.title = "reconstructed set"
for prop in ("tick", "label", "ticklabel")
    for (i, f) in enumerate('x':'z')
        @eval setproperty!(ax_embe, Symbol($f, $prop, :color), cs[$i])
    end
end

# %% Perform animation
# for i in 2:(N-2)*tail
#     update_to_time!(i)
#     yield()
#     # sleep(0.001)
# end

record(fig, projectdir("animations", "6", "embedding.mp4"), 1:(N-2)*tail;
    framerate = 60) do i
    update_to_time!(i)
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# Contains the physical specification of the model

# A physical specification of a model.
# Contains the geometry information, but no discretization parameters.
# The exact model used is defined by the list of terms.
struct Model{T <: Real, VT <: Real}
    # T is the default type to express data, VT the corresponding bare value type (i.e. not dual)

    # Human-readable name for the model (like LDA, PBE, ...)
    model_name::String

    # Lattice and reciprocal lattice vectors in columns
    lattice::Mat3{T}
    recip_lattice::Mat3{T}
    # Dimension of the system; 3 unless `lattice` has zero columns
    n_dim::Int
    # Useful for conversions between cartesian and reduced coordinates
    inv_lattice::Mat3{T}
    inv_recip_lattice::Mat3{T}
    # Volumes
    unit_cell_volume::T
    recip_cell_volume::T

    # Electrons, occupation and smearing function
    n_electrons::Int  # usually consistent with `atoms` field, but doesn't have to

    # spin_polarization values:
    #     :none       No spin polarization, αα and ββ density identical,
    #                 αβ and βα blocks zero, 1 spin component treated explicitly
    #     :collinear  Spin is polarized, but everywhere in the same direction.
    #                 αα ̸= ββ, αβ = βα = 0, 2 spin components treated
    #     :full       Generic magnetization, non-uniform direction.
    #                 αβ, βα, αα, ββ all nonzero, different (not supported)
    #     :spinless   No spin at all ("spinless fermions", "mathematicians' electrons").
    #                 The difference with :none is that the occupations are 1 instead of 2
    spin_polarization::Symbol
    n_spin_components::Int  # 2 if :collinear, 1 otherwise

    # If temperature==0, no fractional occupations are used.
    # If temperature is nonzero, the occupations are
    # `fn = max_occ*smearing((εn-εF) / temperature)`
    temperature::T
    smearing::Smearing.SmearingFunction # see Smearing.jl for choices

    # Particle types (elements) and particle positions and in fractional coordinates.
    # Possibly empty. It's up to the `term_types` to make use of this (or not).
    # `atom_groups` contains the groups of indices into atoms and positions, which
    # point to identical atoms. It is computed automatically on Model construction and may
    # be used to optimise the term instantiation.
    atoms::Vector{Element}
    positions::Vector{Vec3{T}}  # positions[i] is the location of atoms[i] in fract. coords
    atom_groups::Vector{Vector{Int}}  # atoms[i] == atoms[j] for all i, j in atom_group[α]

    # each element t must implement t(basis), which instantiates a
    # term in a given basis and gives back a term (<: Term)
    # see terms.jl for some default terms
    term_types::Vector

    # list of symmetries of the model
    symmetries::Vector{SymOp{VT}}
end

_is_well_conditioned(A; tol=1e5) = (cond(A) <= tol)

"""
    Model(lattice, atoms, positions; n_electrons, magnetic_moments, terms, temperature,
          smearing, spin_polarization, symmetries)

Creates the physical specification of a model (without any discretization information).

`n_electrons` is taken from `atoms` if not specified.

`spin_polarization` is :none by default (paired electrons)
unless any of the elements has a non-zero initial magnetic moment.
In this case the spin_polarization will be :collinear.

`magnetic_moments` is only used to determine the symmetry and the
`spin_polarization`; it is not stored inside the datastructure.

`smearing` is Fermi-Dirac if `temperature` is non-zero, none otherwise

The `symmetries` kwarg allows (a) to pass `true` / `false` to enable / disable
the automatic determination of lattice symmetries or (b) to pass an explicit list
of symmetry operations to use for lowering the computational effort.
The default behaviour is equal to `true`, namely that the code checks the
specified model in form of the Hamiltonian `terms`, `lattice`, `atoms` and `magnetic_moments`
parameters and from these automatically determines a set of symmetries it can safely use.
If you want to pass custom symmetry operations (e.g. a reduced or extended set) use the
`symmetry_operations` function. Notice that this may lead to wrong results if e.g. the
external potential breaks some of the passed symmetries. Use `false` to turn off
symmetries completely.
"""
function Model(lattice::AbstractMatrix{T},
               atoms::Vector{<:Element}=Element[],
               positions::Vector{<:AbstractVector}=Vec3{T}[];
               model_name="custom",
               n_electrons::Int=sum(n_elec_valence, atoms; init=0),
               magnetic_moments=T[],
               terms=[Kinetic()],
               temperature=zero(T),
               smearing=nothing,
               spin_polarization=default_spin_polarization(magnetic_moments),
               symmetries=default_symmetries(lattice, atoms, positions, magnetic_moments,
                                             spin_polarization, terms),
               ) where {T <: Real}
    lattice = Mat3{T}(lattice)
    temperature = T(austrip(temperature))

    # Atoms and electrons
    if length(atoms) != length(positions)
        error("Length of atoms and positions vectors need to agree.")
    end
    n_electrons < 0 && error("n_electrons should be non-negative. Ensure to provide a " *
                             "non-empty atoms list or an appropriate `n_electrons` kwarg.")
    isempty(terms) && error("Model without terms not supported.")
    atom_groups = [findall(Ref(pot) .== atoms) for pot in Set(atoms)]

    # Special handling of 1D and 2D systems, and sanity checks
    n_dim = count(!iszero, eachcol(lattice))
    n_dim > 0 || error("Check your lattice; we do not do 0D systems")
    for i = n_dim+1:3
        norm(lattice[:, i]) == norm(lattice[i, :]) == 0 || error(
            "For 1D and 2D systems, the non-empty dimensions must come first")
    end
    _is_well_conditioned(lattice[1:n_dim, 1:n_dim]) || @warn (
        "Your lattice is badly conditioned, the computation is likely to fail.")

    # Note: In the 1D or 2D case, the volume is the length/surface
    inv_lattice       = compute_inverse_lattice(lattice)
    recip_lattice     = compute_recip_lattice(lattice)
    inv_recip_lattice = compute_inverse_lattice(recip_lattice)
    unit_cell_volume  = compute_unit_cell_volume(lattice)
    recip_cell_volume = compute_unit_cell_volume(recip_lattice)

    # Spin polarization
    spin_polarization in (:none, :collinear, :full, :spinless) ||
        error("Only :none, :collinear, :full and :spinless allowed for spin_polarization")
    spin_polarization == :full && error("Full spin polarization not yet supported")
    !isempty(magnetic_moments) && !(spin_polarization in (:collinear, :full)) && @warn(
        "Non-empty magnetic_moments on a Model without spin polarization detected."
    )
    n_spin = length(spin_components(spin_polarization))

    if isnothing(smearing)
        @assert temperature >= 0
        # Default to Fermi-Dirac smearing when finite temperature
        smearing = temperature > 0.0 ? Smearing.FermiDirac() : Smearing.None()
    end

    if !allunique(string.(nameof.(typeof.(terms))))
        error("Having several terms of the same name is not supported.")
    end

    # Determine symmetry operations to use
    if symmetries === true
        symmetries = default_symmetries(lattice, atoms, positions, magnetic_moments,
                                        spin_polarization, terms)
    elseif symmetries === false
        symmetries = [one(SymOp)]
    end
    @assert !isempty(symmetries)  # Identity has to be always present.

    Model{T,value_type(T)}(model_name,
                           lattice, recip_lattice, n_dim, inv_lattice, inv_recip_lattice,
                           unit_cell_volume, recip_cell_volume,
                           n_electrons, spin_polarization, n_spin, T(temperature), smearing,
                           atoms, positions, atom_groups, terms, symmetries)
end
function Model(lattice::AbstractMatrix{<:Integer}, atoms::Vector{<:Element},
               positions::Vector{<:AbstractVector}; kwargs...)
    Model(Float64.(lattice), atoms, positions; kwargs...)
end
function Model(lattice::AbstractMatrix{<:Quantity}, atoms::Vector{<:Element},
               positions::Vector{<:AbstractVector}; kwargs...)
    Model(austrip.(lattice), atoms, positions; kwargs...)
end

normalize_magnetic_moment(::Nothing)::Vec3{Float64}          = (0, 0, 0)
normalize_magnetic_moment(mm::Number)::Vec3{Float64}         = (0, 0, mm)
normalize_magnetic_moment(mm::AbstractVector)::Vec3{Float64} = mm

"""
:none if no element has a magnetic moment, else :collinear or :full
"""
function default_spin_polarization(magnetic_moments)
    isempty(magnetic_moments) && return :none
    all_magmoms = normalize_magnetic_moment.(magnetic_moments)
    all(iszero, all_magmoms) && return :none
    all(iszero(magmom[1:2]) for magmom in all_magmoms) && return :collinear

    :full
end

"""
Default logic to determine the symmetry operations to be used in the model.
"""
function default_symmetries(lattice, atoms, positions, magnetic_moments, spin_polarization,
                            terms; tol_symmetry=SYMMETRY_TOLERANCE)
    dimension = count(!iszero, eachcol(lattice))
    if spin_polarization == :full || dimension != 3
        return [one(SymOp)]  # Symmetry not supported in spglib
    elseif spin_polarization == :collinear && isempty(magnetic_moments)
        # Spin-breaking due to initial magnetic moments cannot be determined
        return [one(SymOp)]
    elseif any(breaks_symmetries, terms)
        return [one(SymOp)]  # Terms break symmetry
    end

    # Standard case from here on:
    if length(positions) != length(atoms)
        error("Length of atoms and positions vectors need to agree.")
    end
    if !isempty(magnetic_moments) && length(magnetic_moments) != length(atoms)
        error("Length of atoms and magnetic_moments vectors need to agree.")
    end
    magnetic_moments = normalize_magnetic_moment.(magnetic_moments)
    symmetry_operations(lattice, atoms, positions, magnetic_moments; tol_symmetry)
end

"""
Maximal occupation of a state (2 for non-spin-polarized electrons, 1 otherwise).
"""
function filled_occupation(model)
    if model.spin_polarization in (:spinless, :collinear)
        return 1
    elseif model.spin_polarization == :none
        return 2
    else
        error("Not implemented $(model.spin_polarization)")
    end
end

"""
Explicit spin components of the KS orbitals and the density
"""
function spin_components(spin_polarization::Symbol)
    spin_polarization == :collinear && return (:up, :down  )
    spin_polarization == :none      && return (:both,      )
    spin_polarization == :spinless  && return (:spinless,  )
    spin_polarization == :full      && return (:undefined, )
end
spin_components(model::Model) = spin_components(model.spin_polarization)

# prevent broadcast
import Base.Broadcast.broadcastable
Base.Broadcast.broadcastable(model::Model) = Ref(model)

#=
There are two types of quantities, depending on how they transform under change of coordinates.

Positions transform with the lattice: r_cart = lattice * r_red. We term them vectors.

Linear forms on vectors (anything that appears in an expression f⋅r) transform
with the inverse lattice transpose: if f_cart ⋅ r_cart = f_red ⋅ r_red, then
f_cart = lattice' \ f_red. We term them covectors.
Examples of covectors are forces.

Reciprocal vectors are a special case: they are covectors, but conventionally have an
additional factor of 2π in their definition, so they transform rather with 2π times the
inverse lattice transpose: q_cart = 2π lattice' \ q_red = recip_lattice * q_red.
=#
vector_red_to_cart(model::Model, rred)        = model.lattice * rred
vector_cart_to_red(model::Model, rcart)       = model.inv_lattice * rcart
covector_red_to_cart(model::Model, fred)      = model.inv_lattice' * fred
covector_cart_to_red(model::Model, fcart)     = model.lattice' * fcart
recip_vector_red_to_cart(model::Model, qred)  = model.recip_lattice * qred
recip_vector_cart_to_red(model::Model, qcart) = model.inv_recip_lattice * qcart

#=
Transformations on vectors and covectors are matrices and comatrices.

Consider two covectors f and g related by a transformation matrix B. In reduced
coordinates g_red = B_red f_red and in cartesian coordinates we want g_cart = B_cart f_cart.
From g_cart = L⁻ᵀ g_red = L⁻ᵀ B_red f_red = L⁻ᵀ B_red Lᵀ f_cart, we see B_cart = L⁻ᵀ B_red Lᵀ.

Similarly for two vectors r and s with s_red = A_red r_red and s_cart = A_cart r_cart:
s_cart = L s_red = L A_red r_red = L A_red L⁻¹ r_cart, thus A_cart = L A_red L⁻¹.

Examples of matrices are the symmetries in real space (W)
Examples of comatrices are the symmetries in reciprocal space (S)
=#
matrix_red_to_cart(model::Model, Ared)    = model.lattice * Ared * model.inv_lattice
matrix_cart_to_red(model::Model, Acart)   = model.inv_lattice * Acart * model.lattice
comatrix_red_to_cart(model::Model, Bred)  = model.inv_lattice' * Bred * model.lattice'
comatrix_cart_to_red(model::Model, Bcart) = model.lattice' * Bcart * model.inv_lattice'


text:
using SymbolicUtils, Test
using SymbolicUtils: Term, symtype
@testset "polyform" begin
    @syms a b c d
    @test expand(a * (b + -1 * c) + -1 * (b * a + -1 * c * a)) == 0
    @eqtest expand(sin(a+b)+sin(c+d)) == sin(a+b) + sin(c+d)
    @eqtest simplify(expand(sin((a+b)^2)^2)) == simplify(sin(a^2+2*(b*a)+b^2)^2)
    @test simplify(expand(sin((a+b)^2)^2 + cos((a+b)^2)^2)) == 1
    @syms x1::Real f(::Real)::Real

    # issue 193
    @test isequal(expand(f(x1 + 2.0)), f(2.0 + x1))
    @test symtype(expand(f(x1 + 2.0))) == Real

    # cleanup rules
    @test expand(Term{Number}(identity, 0)) == 0
    @test expand(Term{Number}(one, 0)) == 1
    @test expand(Term{Number}(zero, 0)) == 0
    @test expand(identity(a * b) - b * a) == 0
    @test expand(a * b - b * a) == 0
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Parse SceneLang commands

"""
    parse_set_command(stream::InputStream, scene::Scene)

Push to the 'scene.variables' [`IdTable`](@ref) a constructed value and its identifier variable.
"""
function parse_set_command(stream::InputStream, scene::Scene)
    table = scene.variables
    expect_command(stream, SET)
    while true
        id = read_token(stream)
        isa(id.value, Identifier) || (unread_token(stream, id); break)
        id_name = id.value.value
        if (type = findfirst(d -> haskey(d, id_name), table)) |> !isnothing
            preexisting = table[type][id_name]
            iszero(preexisting.loc.line_num) ||
                throw(IdentifierRedefinition(id.loc, "Identifier '$(id_name)' has alredy been set at\n$(preexisting.loc)\nIf you want to redefine it first UNSET it.", id.length))
            parse_constructor(stream, scene)
            continue # if identifier was defined at the command line level throw no error and ignore the set instruction
        end
        value, id_type = parse_constructor(stream, scene)
        haskey(table, id_type) ?
            push!(table[id_type], id_name => ValueLoc(value, copy(id.loc))) :
            push!(table, id_type => Dict([id_name => ValueLoc(value, copy(id.loc))]))
    end
end

"""
    parse_unset_command(stream::InputStream, scene::Scene)

Pop identifier variable from the `scene.variables` [`IdTable`](@ref).
"""
function parse_unset_command(stream::InputStream, scene::Scene)
    table = scene.variables
    expect_command(stream, UNSET)
    while true
        id = read_token(stream)
        isa(id.value, Identifier) || (unread_token(stream, id); break)
        id_name = id.value.value
        type = findfirst(d -> haskey(d, id_name), table)
        isnothing(type) &&
            throw(UndefinedIdentifier(id.loc,"Undefined variable '$id_name'" ,id.length))
        pop!(table[type], id_name)
    end
    return
end

"""
    parse_dump_command(stream::InputStream, scene::Scene)

Show the contents of the [`Scene`](@ref). What is shown depends on the [`Keyword`](@ref) following the `DUMP` [`Command`](@ref).

The valid [`Keyword`](@ref)s are `.ALL` for showing the whole [`Scene`](@ref), or the name of one of its fields (in lowercase letters) to show that specific field.
"""
function parse_dump_command(stream::InputStream, scene::Scene)
    table = scene.variables
    expect_command(stream, DUMP)
    valid_keywords = (:ALL, :variables, :world, :lights, :image, :camera, :renderer)
    next_token = read_token(stream)
    if isa(next_token.value, Keyword)
        unread_token(stream, next_token)
        keyword = expect_keyword(stream, valid_keywords).value.value
        keyword == :ALL ?
            display(scene) :
            display(getproperty(scene, keyword))
    elseif isa(next_token.value, Identifier)
        unread_token(stream, next_token)
        id_name = expect_identifier(stream).value.value
        type = findfirst(d -> haskey(d, id_name), table)
        display(table[type][id_name])
    else
        throw(WrongTokenType(next_token.loc, "Expected either a keyword or a valid identifier instead of '$(typeof(next_token.value))'\n"*
                             "Valid keywords: \n\t$(join(valid_keywords, "\n\t"))",next_token.length))
    end
end

"""
    parse_spawn_command(stream::InputStream, scene::Scene)

Push the given `ShapeType` or `LightType` to the `scene.world` [`World`](@ref) or `scene.lights` [`Lights`](@ref) respectively.

See also: [`Scene`](@ref)
"""
function parse_spawn_command(stream::InputStream, scene::Scene)
    table = scene.variables
    expect_command(stream, SPAWN)
    next_token = read_token(stream)
    unread_token(stream, next_token)
    next_val = next_token.value
    (isa(next_val, Identifier) ||
     isa(next_val, LiteralType) ||
     next_val ∈ (UNITE, INTERSECT, DIFF, FUSE)) ||
        throw(WrongTokenType(next_token.loc,"Expected either a constructor or a valid identifier instead of '$(typeof(next_val))'", next_token.length))
    while true
        if isa(next_val, Identifier)
            id_name = expect_identifier(stream).value.value
            type = findfirst(d -> haskey(d, id_name), table)
            isnothing(type) &&
                throw(UndefinedIdentifier(next_token.loc,"Undefined variable '$id_name'",next_token.length))
            if type == ShapeType
                shape = table[type][id_name].value
                push!(scene.world, shape)
            elseif type == LightType
                light = table[type][id_name].value
                push!(scene.lights, light)
            else
                throw(WrongValueType(next_token.loc, "Identifier '$id_name' stores a non-spawnable '$type' object\n" *
                                    "Variable '$id_name' defined at $(table[type][id_name].loc)\n" *
                                     "Spawnable types are:\n\tShapeType\n\tLightType", next_token.length))
            end
        elseif isa(next_val, LiteralType)
            type_token = expect_type(stream, (ShapeType, LightType))
            unread_token(stream, type_token)
            type = type_token.value
            if type == ShapeType
                shape = parse_shape(stream, scene)
                push!(scene.world, shape)
            elseif type == LightType
                light = parse_light(stream, scene)
                push!(scene.lights, light)
            else
                @assert false "@ $(next_token.loc): expect_type returned a non-spawnable type '$type'"
            end
        elseif next_val ∈ (UNITE, INTERSECT, DIFF, FUSE)
            push!(scene.world, parse_shape_from_command(stream, scene))
        else
            break
        end
        next_token = read_token(stream)
        unread_token(stream, next_token)
        next_val = next_token.value
    end
    return
end

"""
    parse_using_command(stream::InputStream, scene::Scene)

Set the given `CameraType`, `ImageType` or `RendererType` to the `scene.camera`, `scene.image`, or `scene.renderer` respectively.

Can only be used once per type in a SceneLang script.

See also: [`Scene`](@ref)
"""
function parse_using_command(stream::InputStream, scene::Scene)
    table = scene.variables
    expect_command(stream, USING)
    next_token = read_token(stream)
    unread_token(stream, next_token)

    (isa(next_token.value, Identifier) ||
     isa(next_token.value, LiteralType) ||
     next_token.value ∈ (LOAD,)) ||
        throw(WrongTokenType(next_token.loc,"Expected either a constructor or a valid identifier instead of '$(typeof(next_token.value))'", next_token.length))

    function already_defined_exception(type::LiteralType)
        @assert type ∈ (CameraType, ImageType, RendererType)
        SettingRedefinition(next_token.loc, "Scene setting of type '$(type)' already in use.", next_token.length)
    end
    while true
        if isa(next_token.value, Identifier)
            id_name = expect_identifier(stream).value.value
            type = findfirst(d -> haskey(d, id_name), table)
            isnothing(type) &&
                throw(UndefinedIdentifier(next_token.loc,"Undefined variable '$id_name'",next_token.length))
            if type == CameraType
                camera = table[type][id_name].value
                isnothing(scene.camera) || throw(already_defined_exception(type))
                scene.camera = camera
            elseif type == ImageType
                image = table[type][id_name].value
                isnothing(scene.image) || throw(already_defined_exception(type))
                scene.image = image
            elseif type == RendererType
                renderer = table[type][id_name].value
                isnothing(scene.renderer) || throw(already_defined_exception(type))
                scene.renderer = renderer
            elseif type == TracerType
                tracer = table[type][id_name].value
                isnothing(scene.tracer) || throw(already_defined_exception(type))
                scene.tracer = tracer
            else
                throw(WrongValueType(next_token.loc, "Variable '$id_name' stores a non-usable '$type' object\n" *
                                        "Variable '$id_name' defined at $(table[type][id_name].loc)\n" *
                                        "Usable types are:\n\tCameraType\n\tImageType\n\tRendererType", next_token.length))
            end
        elseif isa(next_token.value, LiteralType)
            type_token = expect_type(stream, (CameraType, ImageType, RendererType, TracerType))
            unread_token(stream, type_token)
            type =type_token.value
            if type == CameraType
                camera = parse_camera(stream, scene)
                isnothing(scene.camera) || throw(already_defined_exception(type))
                scene.camera = camera
            elseif type == ImageType
                image = parse_explicit_image(stream, scene)
                isnothing(scene.image) || throw(already_defined_exception(type))
                scene.image = image
            elseif type == RendererType
                renderer = parse_renderer_settings(stream, scene)
                isnothing(scene.renderer) || throw(already_defined_exception(type))
                scene.renderer = renderer
            elseif type == TracerType
                tracer = parse_tracer_settings(stream, scene)
                isnothing(scene.tracer) || throw(already_defined_exception(type))
                scene.tracer = tracer
            else
                @assert false "@ $(next_token.loc): expect_type returned a non-spawnable type '$type'"
            end
        elseif isa(next_token.value, Command)
            next_token.value ∈ (SET, UNSET, SPAWN, USING, DUMP) && break
            image = parse_image_from_command(stream, scene)
            isnothing(scene.image) || throw(already_defined_exception(type))
            scene.image = image
        else
            break
        end
        next_token = read_token(stream)
        unread_token(stream, next_token)
    end
    return
end


text:
"""
$(TYPEDSIGNATURES)

Check whether in abstract vector `p` is a positive vector, i.e. that pᵢ≥0∀i.
"""
isposvec(p::AbstractVector{<:Real}) = all(p .≥ zero(eltype(p)))

"""
$(TYPEDSIGNATURES)

Check whether in abstract vector `p` sums to 1.
"""
isnormvec(p::AbstractVector{<:Real}) = isapprox(sum(p), one(eltype(p)))

"""
$(TYPEDSIGNATURES)

Check whether in abstract vector `p` is a probability vector, i.e. that ∑p=1 and pᵢ≥0∀i.
"""
isprobvec(p::AbstractVector{<:Real}) = isposvec(p) && isnormvec(p)

"""
$(TYPEDSIGNATURES)

Return non-zero elements of a vector.
"""
nonzero(p::AbstractVector{<:Real}) = p[p.!=zero(eltype(p))]

"""
$(TYPEDSIGNATURES)

Convert a vector to a vector of integer weights
"""
convert_prob_to_intweight(target::Vector{<:Real}) =
    round.(Int, target / minimum(nonzero(target)))

"""
$(TYPEDSIGNATURES)
"""
calculate_min_blocksize(target::Vector{<:Real}) = sum(convert_prob_to_intweight(target))

function random_sample(p::AbstractVector{<:Real})
    isprobvec(p) ? searchsortedlast([0.0; cumsum(p)], rand()) :
    error("p must be a probability vector")
end
function random_sample(p::AbstractVector{<:Real}, u::Real)
    if u ≤ 0 || u ≥ 1
        error("Must have 0 < u < 1")
    end
    isprobvec(p) ? searchsortedlast([0.0; cumsum(p)], u) :
    error("p must be a probability vector")
end


text:
using Pkg;
Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate()

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter, NearestNeighborDescent

makedocs(
    modules=[NearestNeighborDescent],
    format=Documenter.HTML(analytics="UA-154568964-1",
                           canonical="https://dillondaudert.github.io/NearestNeighborDescent.jl/stable/"),
    sitename="NearestNeighborDescent.jl",
    authors="Dillon Daudert",
    pages = [
        "Home" => "index.md",
        "Basic Usage" => "basic.md",
        "KNN Graphs" => [
            "graphs.md",
        ],
        "Reference" => [
            "Public" => "ref/public.md",
        ],
    ],
)

deploydocs(
    repo="github.com/dillondaudert/NearestNeighborDescent.jl.git",
)


text:
include("HashTest.jl")
include("RandTest.jl")
include("IterationTest.jl")
include("PaddingTest.jl")
include("BlockCiphersTest.jl")
include("CipherModesTest.jl")


text:
using TinyG
using Base.Test

TinyG.download_firmware()


text:
module AddIntegersF90
using CxxInterface

const libAddIntegersF90 = joinpath(pwd(), "libAddIntegersF90")

eval(f90setup())
eval(f90newfile("AddIntegersF90.f90", ""))

eval(f90function(FnName(:add_int, "add_int", libAddIntegersF90), FnResult(Cint, "integer", Int, expr -> :(convert(Int, $expr))),
                 [FnArg(:x, Ref{Cint}, "x", "integer", Integer, identity), FnArg(:y, Ref{Cint}, "y", "integer", Integer, identity)],
                 "add_int = x + y"))
end

################################################################################
# Extract Fortran code
using CxxInterface
using .AddIntegersF90
AddIntegersF90.f90_write_code!()

################################################################################
# Compile Fortran code
# (This fails if there is no Fortran compiler available)
# macOS does not have a Fortran compiler installed by default
if !Sys.isapple()
    using Libdl: dlext
    run(`gfortran -fPIC -c AddIntegersF90.f90`)
    run(`gfortran -shared -o libAddIntegersF90.$dlext AddIntegersF90.o`)

    # Please, DO NOT call a Fortran compiler manually in your own Julia
    # packages. This works only in very controlled environments such as on
    # CI infrastructure. If you do, your package will be fragile, and will
    # create lots of headaches for your users in the wild. Instead, use
    # [BinaryBuilder](https://binarybuilder.org) and store your build
    # recipes on [Yggdrasil](https://github.com/JuliaPackaging/Yggdrasil).

    ################################################################################
    # Call the wrapped function
    # (This fails if the Fortran compiler is not compatible with the
    # currently running Julia executable)

    # Only test cases where the Github CI environment supports this (:i686 is not supported)
    if Sys.ARCH ≡ :x86_64
        @test AddIntegersF90.add_int(2, 3) == 5
    end
end


text:
module CutPools

using JuMP

export Separator, Oracle, CuttingPlane, CutPool,
       separate, separator_callback, addfamily!

abstract Separator
abstract Oracle <: Separator

type CuttingPlane
    ineq::LinearConstraint
    added::Bool
    CuttingPlane(ineq::LinearConstraint) = new(ineq,false)
end

type CutFamily
    cuts::Vector{CuttingPlane}
    name::String
end

type CutPool <: Separator
    pool::Vector{CuttingPlane}
    clear_used::Bool
    groups::Vector{CutFamily}
end
CutPool() = CutPool(CuttingPlane[], true, CutFamily[])

function addfamily!(pool::CutPool, cp::Vector{LinearConstraint}, name::String)
    cuts = [CuttingPlane(c) for c in cp]
    append!(pool.pool, cuts)
    push!(pool.groups, CutFamily(cuts,name))
end

Base.push!(x::CutPool,y) = push!(x.pool,CuttingPlane(y))
Base.append!(x::CutPool,y) = append!(x.pool,[CuttingPlane(c) for c in y])

const ε = 1e-4
const ratio_pool_reshuffle = 0.1

function separate(cutpool::CutPool)
    separated = CuttingPlane[]
    num_sep = 0
    for cut in cutpool.pool
        if cut.added
            num_sep += 1
            continue
        end
        val = getValue(cut.ineq.terms)
        lb, ub = cut.ineq.lb, cut.ineq.ub
        if (lb - val) / abs(min(lb,val,ε)) > ε ||
           (val - ub) / abs(min(ub,val,ε)) > ε
            push!(separated, cut)
            cut.added = true
            num_sep += 1
        end
    end

    if cutpool.clear_used &&
    num_sep > ratio_pool_reshuffle*length(cutpool.pool)
        filter!(cutpool.pool) do cut
            !cut.added
        end
    end

    return separated
end

function _separator_callback(cb, cutpool::CutPool, maxnum::Int)
    tic()
    cuts_added = 0
    cuts = separate(cutpool)
    for cut in cuts
        cuts_added += 1
        addUserCut(cb, cut.ineq)
        cuts_added ≥ maxnum && break
    end
    el = toq()
    # println("Adding $(min(length(cuts),maxnum)) cuts in $el seconds")
end

separator_callback(cutpool::CutPool; maxnum=100) =
    (cb -> _separator_callback(cb, cutpool, maxnum))

end


text:
"""
    Jaro()

Creates the Jaro metric

The Jaro distance is defined as

``1 - (m / |s1| + m / |s2| + (m - t) / m) / 3``

where ``m`` is the number of matching characters and
``t`` is half the number of transpositions.
"""
struct Jaro <: SemiMetric end
isnormalized(::Jaro) = true

## http://alias-i.com/lingpipe/docs/api/com/aliasi/spell/JaroWinklerDistance.html
## accepts any iterator, including AbstractString
function evaluate(dist::Jaro, s1, s2, max_dist = nothing)
    (ismissing(s1) | ismissing(s2)) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    # if both are empty, m = 0 so should be 1.0 according to wikipedia.
    # Add this line so that not the case
    len2 == 0 && return 0.0
    maxdist = max(0, div(len2, 2) - 1)
    flag = fill(false, len2)
    prevstate1 = firstindex(s1)
    i1_match = fill(prevstate1, len1)
    #  m counts number matching characters
    m = 0
    i1 = 1
    i2 = 1
    x1 = iterate(s1)
    x2 = iterate(s2)
    while x1 !== nothing
        ch1, state1 = x1
        if i2 <= i1 - maxdist - 1
            ch2, state2 = x2
            i2 += 1
            x2 = iterate(s2, state2)
        end
        i2curr = i2
        x2curr = x2
        while x2curr !== nothing
            i2curr > i1 + maxdist && break
            ch2, state2 = x2curr
            if (ch1 == ch2) && !flag[i2curr]
                m += 1
                flag[i2curr] = true
                i1_match[m] = prevstate1
                break
            end
            x2curr = iterate(s2, state2)
            i2curr += 1
        end
        x1 = iterate(s1, state1)
        i1 += 1
        prevstate1 = state1
    end
    m == 0 && return 1.0
    # t counts number of transpositions
    t = 0
    i1 = 0
    i2 = 0
    for ch2 in s2
        i2 += 1
        if flag[i2]
            i1 += 1
            t += ch2 != iterate(s1, i1_match[i1])[1]
        end
    end
    return 1.0 - (m / len1 + m / len2 + (m - t/2) / m) / 3.0
end

"""
    Levenshtein()

Creates the Levenshtein metric

The Levenshtein distance is the minimum number of operations (consisting of insertions, deletions,
substitutions of a single character) required to change one string into the other.
"""
struct Levenshtein <: Metric end

## Source: http://blog.softwx.net/2014/12/optimizing-levenshtein-algorithm-in-c.html
# Return max_dist +1 if distance higher than max_dist
# This makes it possible to differentiate distance equalt to max_dist vs strictly higher
# This is important for find_all
## accepts any iterator, including AbstractString
function evaluate(dist::Levenshtein, s1, s2, max_dist = nothing)
    (ismissing(s1) | ismissing(s2)) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    max_dist !== nothing && len2 - len1 > max_dist && return max_dist + 1
    # prefix common to both strings can be ignored
    k, x1, x2start = common_prefix(s1, s2)
    x1 == nothing && return len2 - k
    # distance initialized to first row of matrix
    # => distance between "" and s2[1:i}
    v = collect(1:(len2-k))
    current = 0
    i1 = 1
    while x1 !== nothing
        ch1, state1 = x1
        left = i1 - 1
        current = i1 - 1
        min_dist = i1 - 2
        i2 = 1
        x2 = x2start
        while x2 !== nothing
            ch2, state2 = x2
            #  update
            above, current, left = current, left, v[i2]
            if ch1 != ch2
                current = min(current + 1, above + 1, left + 1)
            end
            min_dist = min(min_dist, left)
            v[i2] = current
            x2 = iterate(s2, state2)
            i2 += 1
        end
        max_dist !== nothing && min_dist > max_dist && return max_dist + 1
        x1 = iterate(s1, state1)
        i1 += 1
    end
    max_dist !== nothing && current > max_dist && return max_dist + 1
    return current
end

"""
    DamerauLevenshtein()

Creates the DamerauLevenshtein metric

The DamerauLevenshtein distance is the minimum number of operations (consisting of insertions,
deletions or substitutions of a single character, or transposition of two adjacent characters)
required to change one string into the other.
"""
struct DamerauLevenshtein <: SemiMetric end

## http://blog.softwx.net/2015/01/optimizing-damerau-levenshtein_15.html
## accepts any iterator, including AbstractString
function evaluate(dist::DamerauLevenshtein, s1, s2, max_dist = nothing)
    (ismissing(s1) | ismissing(s2)) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    max_dist !== nothing && len2 - len1 > max_dist && return max_dist + 1
    # prefix common to both strings can be ignored
    k, x1, x2start = common_prefix(s1, s2)
    (x1 == nothing) && return len2 - k
    v = collect(1:(len2-k))
    w = similar(v)
    if max_dist !== nothing
        offset = 1 + max_dist - (len2 - len1)
        i2_start = 1
        i2_end = max_dist
    end
    i1 = 1
    current = i1
    prevch1, = x1
    while x1 !== nothing
        ch1, state1 = x1
        left = (i1 - 1)
        current = i1
        nextTransCost = 0
        prevch2, = x2start
        if max_dist !== nothing
            i2_start += (i1 > offset) ? 1 : 0
            i2_end = min(i2_end + 1, len2)
        end
        x2 = x2start
        i2 = 1
        while x2 !== nothing
            ch2, state2 = x2
            if max_dist == nothing || (i2_start <= i2 <= i2_end)
                above = current
                thisTransCost = nextTransCost
                nextTransCost = w[i2]
                # cost of diagonal (substitution)
                w[i2] = current = left
                # left now equals current cost (which will be diagonal at next iteration)
                left = v[i2]
                if ch1 != ch2
                    # insertion
                    if left < current
                        current = left
                    end
                    # deletion
                    if above < current
                        current = above
                    end
                    current += 1
                    if (i1 != 1) & (i2 != 1) & (ch1 == prevch2) & (prevch1 == ch2)
                        thisTransCost += 1
                        if thisTransCost < current
                            current = thisTransCost
                        end
                    end
                end
                v[i2] = current
            end
            x2 = iterate(s2, state2)
            i2 += 1
            prevch2 = ch2
        end
        max_dist !== nothing && v[i1 + len2 - len1] > max_dist && return max_dist + 1
        x1 = iterate(s1, state1)
        i1 += 1
        prevch1 = ch1
    end
    max_dist !== nothing && current > max_dist && return max_dist + 1
    return current
end

"""
    RatcliffObershelp()

Creates the RatcliffObershelp metric

The distance between two strings is defined as one minus  the number of matching characters
divided by the total number of characters in the two strings. Matching characters are those
in the longest common subsequence plus, recursively, matching characters in the unmatched
region on either side of the longest common subsequence.
"""
struct RatcliffObershelp <: SemiMetric end

isnormalized(::RatcliffObershelp) = true

function evaluate(dist::RatcliffObershelp, s1, s2, max_dist = nothing)
    (ismissing(s1) | ismissing(s2)) && return missing
    n_matched = sum(last.(matching_blocks(s1, s2)))
    len1, len2 = length(s1), length(s2)
    len1 + len2 == 0 ? 0. : 1.0 - 2 *  n_matched / (len1 + len2)
end

function matching_blocks(s1, s2)
    matching_blocks!(Set{Tuple{Int, Int, Int}}(), s1, s2, length(s1), length(s2), 1, 1)
end

function matching_blocks!(x::Set{Tuple{Int, Int, Int}}, s1, s2,
    len1::Integer, len2::Integer, start1::Integer, start2::Integer)
    a = longest_common_pattern(s1, s2, len1 , len2)
    # exit if there is no common substring
    a[3] == 0 && return x
    # add the info of the common to the existing set
    push!(x, (a[1] + start1 - 1, a[2] + start2 - 1, a[3]))
    # add the longest common substring that happens before
    s1before = _take(s1, a[1] - 1)
    s2before = _take(s2, a[2] - 1)
    matching_blocks!(x, s1before, s2before, a[1] - 1, a[2] - 1, start1, start2)
    # add the longest common substring that happens after
    s1after = _drop(s1, a[1] + a[3] - 1)
    s2after = _drop(s2, a[2] + a[3] - 1)
    matching_blocks!(x, s1after, s2after, len1 - (a[1] + a[3]) + 1,
        len2 - (a[2] + a[3]) + 1, start1 + a[1] + a[3] - 1, start2 + a[2] + a[3] - 1)
    return x
end

function longest_common_pattern(s1, s2, len1::Integer, len2::Integer)
    if len1 > len2
        start2, start1, len = longest_common_pattern(s2, s1, len2, len1)
    else
        start1, start2, len = 0, 0, 0
        p = zeros(Int, len2)
        i1 = 0
        for ch1 in s1
            i1 += 1
            i2 = 0
            oldp = 0
            for ch2 in s2
                i2 += 1
                newp = 0
                if ch1 == ch2
                    newp = oldp > 0 ? oldp : i2
                    currentlength = (i2 - newp + 1)
                    if currentlength > len
                        start1, start2, len = i1 - currentlength + 1, newp, currentlength
                    end
                end
                p[i2], oldp = newp, p[i2]
            end
        end
    end
    return start1, start2, len
end


text:
``` this is a julia source file, and random stuff
```

using InstantiateFromURL
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.4.0")

x = 2

using Plots
function snapabove(g, a)
  function f(x)
    if x > a     # "a" is captured in the closure f
      return g(x)
    else
      return g(a)
    end
  end
  return f  # closure with the embedded a
end

f(x) = x^2
h = snapabove(f, 2.0)

gr(fmt=:png);
plot(h, 0.0:0.1:5.0)

for i in 1:3
    print(i)
end
x_values = 1:5

for x = x_values
  println(x*x)
end

[ i + j for i in 1:3, j in 4:6 ]
function foo(x)
    if x > 0
        return "positive"
    end
    return "nonpositive"
end

x_vec = [0.0 2.0 3.0 4.0]
y_vec = sin.(x_vec)

using LinearAlgebra, Statistics

function chisq(k)
  @assert k > 0
  z = randn(k)
  return sum(z -> z^2, z)
end

y = chisq.([2 4 6])

f(x, y) = [1, 2, 3] ⋅ x + y   # "⋅" can be typed by \cdot<tab>
f([3, 4, 5], 2)   # uses vector as first parameter
f.(Ref([3, 4, 5]), [2, 3])

twice(f,x) = f(f(x))
f(x) = x^2
@show twice(f, 2.0)

twice(x -> x^2, 2.0)

a = 5
g(x) = a * x
twice(g,2.0)

using Expectations, Distributions
@show d = Exponential(2.0)
f(x) = x^2
@show expectation(f, d);  # E(f(x))

function multiplyit(a, g)
    return x -> a * g(x)  # function with `g` used in the closure
end
f(x) = x^2
h = multiplyit(2.0, f)    # use our quadratic, returns a new function which doubles the result
h(2)     #

a = zeros(2,2)
println(a)
x = Array{Float64}(undef, 2, 2)
println(x)
x[1,1] = 1;
println(x)

a = [1 2; 3 4]
a[:, 2] = [4 5]

a[:,2] = [3;3]
println(a)

a = [1 2; 3 4]
b = a' # transpose
c = Matrix(b)  # convert to matrix
d = collect(b)
a = [0,1,-1]
b = sort(a, rev=false)
println(b)
println(a)

a = zeros(1,3,4);
b = dropdims(a,dims = 1)

using Parameters
function f(parameters)
    @unpack α, β = parameters # poor style, error prone if adding parameters
    return α + β
end
parameters = (α = 0.1, γ = "test", β = 0.2)
f(parameters)

function f(x)
  if x > 0
    return x
  else
    return NaN
  end
end

using Base: show_supertypes
show_supertypes(Int32)
@show subtypes(Int32)

struct Foo
    a::Float64
    b::Int64
    c::Vector{Float64}
end
foo_nt = FooNotTyped(2.0, 3, [1.0, 2.0, 3.0])
typeof(foo_nt)

f(x) = x
x_iv = [1.0]
#f(x) = 1.1 * x # fixed-point blows

# BAD!
xstar = fixedpoint(f, x_iv).zero # assumes convergence
xsol = nlsolve(f, x_iv).zero # assumes convergence

# GOOD!
result = fixedpoint(f, x_iv)
converged(result) || error("Failed to converge in $(result.iterations) iterations")
xstar = result.zero

f(x) = 1.1 * x

result = nlsolve(f, x_iv)
converged(result) || error("Failed to converge in $(result.iterations) iterations")
xsol = result.zero

using LinearAlgebra, Statistics
using QuantEcon, QuadGK, FastGaussQuadrature, Distributions, Expectations
using Interpolations, Plots, LaTeXStrings, ProgressMeter

li = LinearInterpolation(x, y)
li_spline = CubicSplineInterpolation(x, y)

@show li(0.3) # evaluate at a single point

scatter(x, y, label = "sampled data", markersize = 4)
plot!(xf, li.(xf), label = "linear")
plot!(xf, li_spline.(xf), label = "spline")

f(x) = x > 0.0 ? x : nothing
@code_warntype f(1)

function f(a, b)
    y = (a + 8b)^2
    return 7y
end

@code_native f(1.0,2.0)

max_stars_repo_name: robsunch/quanteconuser

max_issues_repo_name: robsunch/quanteconuser

max_forks_repo_name: robsunch/quanteconuser


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

"""
    instantiate(p)

Create a suitable object according to the type specified in the ProtoBuf object `p`.
Return value will be a subtype of `AbstractType`.
"""
function instantiate(t::Proto._Type)
    label = hasproperty(t, :scope) ? Label(t.scope, t.label) : Label(t.label)
    if t.encoding == Proto.Type_Encoding.ROLE_TYPE
        return RoleType(label, t.root)
    elseif t.encoding == Proto.Type_Encoding.RELATION_TYPE
        return RelationType(label, t.root)
    elseif t.encoding == Proto.Type_Encoding.ATTRIBUTE_TYPE
        return AttributeType(label, t.root, t.value_type)
    elseif t.encoding == Proto.Type_Encoding.ENTITY_TYPE
        return EntityType(label, t.root)
    elseif t.encoding == Proto.Type_Encoding.THING_TYPE
        return ThingType(label, t.root)
    else
        throw(TypeDBClientException(CONCEPT_BAD_ENCODING, t.encoding))
    end
end

"""
    encoding(t)

Return the ProtoBuf encoding value for the given concept type `t`.
"""
function encoding end

encoding(t::AbstractEntityType) = Proto.Type_Encoding.ENTITY_TYPE
encoding(t::AbstractRelationType) = Proto.Type_Encoding.RELATION_TYPE
encoding(t::AbstractAttributeType) = Proto.Type_Encoding.ATTRIBUTE_TYPE
encoding(t::AbstractRoleType) = Proto.Type_Encoding.ROLE_TYPE
encoding(t::AbstractThingType) = Proto.Type_Encoding.THING_TYPE

# TODO The UNRECOGNIZED enum is missing; need to regenerate from proto?
# encoding(t::AbstractType) = Proto.Type_Encoding.UNRECOGNIZED

# I don't think this is unnecessary since the user can broadcast themselves
# proto(types::AbstractVector{T}) where {T <: AbstractType} = proto.(types)