text:
using Test
using JSON3
using Pinecone

APIKEY = ENV["PINECONE_API_KEY"]
CLOUDENV="us-west1-gcp"

context = Pinecone.init(APIKEY, CLOUDENV)

INDEX = "filter-example"
NAMESPACE = "mynamespace"
index = Pinecone.Index(INDEX);

@testset verbose = true "Create/Delete" begin
    testindexname = "unittestindex"
    @testset "Create" begin
       #delete in case already present from previous failure
       result = Pinecone.delete_index(context, Pinecone.Index(testindexname))
       #sleep to wait for delete to go thru, backend takes a bit
       sleep(10)
       indexconfig = Dict{String, Any}("k_bits"=>512, "hybrid"=>true)
       result = Pinecone.create_index(context, testindexname, 10, metric="euclidean", indextype="approximated",pods=1,replicas=2, shards=1, indexconfig=indexconfig)
       println("CREATE(): ", result)
       @test result == true

       #max dims is 10000, check you cannot create 10001
       @test_throws ArgumentError Pinecone.create_index(context, testindexname, 10001)
       #broken pod type
       @test_throws ArgumentError Pinecone.create_index(context, "shouldfail", 10, metric="euclidean", indextype="approximated",pods=1,replicas=2, shards=1, podtype="fail", indexconfig=indexconfig)
       #bad pods
       @test_throws ArgumentError Pinecone.create_index(context, "shouldfail", 10, metric="euclidean", indextype="approximated",pods=0,replicas=2, shards=1, podtype="p1", indexconfig=indexconfig)
       #bad shards
       @test_throws ArgumentError Pinecone.create_index(context, "shouldfail", 10, metric="euclidean", indextype="approximated",pods=1,replicas=2, shards=0, podtype="p1", indexconfig=indexconfig)
       #bad pods and shards
       @test_throws ArgumentError Pinecone.create_index(context, "shouldfail", 10, metric="euclidean", indextype="approximated",pods=0,replicas=2, shards=0, podtype="p1", indexconfig=indexconfig)
       #bad metric
       @test_throws ArgumentError Pinecone.create_index(context, "shouldfail", 10, metric="badmetric", indextype="approximated",pods=1,replicas=2, shards=1, podtype="p1", indexconfig=indexconfig)

      end
    @testset "Delete" begin
       result = Pinecone.delete_index(context, Pinecone.Index(testindexname))
       println("DELETE(): ", result)
       @test result == true
    end
 end


text:
@testset "Keywords" begin

using SExpressions.Keywords

@test Keyword("test") == Keyword("test")
@test Keyword("test") ≠ Keyword("Test")

x = Dict(Keyword("test") => "hello",
         Keyword("test2") => "world",
         Keyword("test3") => "!")

@test x[Keyword("test")] == "hello"
x[Keyword("test")] = "modified"
@test x[Keyword("test")] == "modified"

end

ext: jl


text:
#saha function
function fHI(T,ne)
     χ   = 13.60
    temp = 15.6826 - log10(ne) + log10(2) + 1.5 * log10(T) - χ * (5040 / T)

    return temp
end

function fHeI(T,ne)
     χ   = 24.58
    temp = 15.6826 - log10(ne) + log10(2) + 1.5 * log10(T) - χ * (5040 / T)

    return temp
end

function fHeII(T,ne)
     χ   = 54.41
    temp = 15.6826 - log10(ne) + log10(2) + 1.5 * log10(T) - χ * (5040 / T)

    return temp
end

function saha(T,N;ionization = false)
    Ne  = N * 0.1
    X   = 0.9
    Y   = 0.1
    NH  = X * N
    NHe = Y * N
    δ   = 1e-5

    X1  = fHI(  T, Ne)
    X2  = fHeI( T, Ne)
    X3  = fHeII(T, Ne)

    Y1  = 10^(X1) #H0
    Y2  = 10^(X2) #He0
    Y3  = 10^(X3) #He+

    X4  = X * N  / (1 + Y1)
    X5  = Y * N  / (1 + Y2 * (1 + Y3))
    Ne_last = 0

    AHI    = 0
    AHII   = 0
    AHeI   = 0
    AHeII  = 0
    AHeIII = 0
    ANe    = 0

    while true
        HI    = X4
        HII   = X * N - HI
        HeI   = X5
        HeII  = Y2 * HeI
        HeIII = Y3 * HeII
        Ne    = HII + HeII + 2 * HeIII

        if (1 - Ne / Ne_last) < δ

            AHI    = HI
            AHII   = HII
            AHeI   = HeI
            AHeII  = HeIII
            AHeIII = HeIII
            ANe    = Ne

            break
        else
            Ne_last = (Ne_last + Ne) /2
        end
    end

    if ionization == true
        return [AHI,AHII,AHeI,AHeII,AHeIII,ANe]
    else
        return ANe, AHI, AHII, AHeI, AHeII, AHeIII
    end
end

export fHI, fHeI, fHeII, saha


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "GeographicLibWrapper"
version = v"0.1.0"

julia_version = v"1.6"

sources = [DirectorySource("./src")]

script = raw"""
mkdir -p build
cd build
cmake \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DJulia_PREFIX=${prefix} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
cmake --build . --config Release --target install -- -j${nproc}
install_license /usr/share/licenses/MIT
"""

include("../../L/libjulia/common.jl")
platforms = expand_cxxstring_abis(libjulia_platforms(julia_version))

products = [LibraryProduct("libGeographicLibWrapper", :libGeographicLibWrapper)]

dependencies = [
    Dependency("GeographicLib_jll"),
    Dependency("libcxxwrap_julia_jll"),
    BuildDependency(PackageSpec(; name="libjulia_jll", version=julia_version)),
]

build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    preferred_gcc_version=v"7",
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
const _SUPPORTED_FORMATS = ("csv",)

mutable struct SimulationResultsExport
    models::Dict{Symbol, ProblemResultsExport}
    start_time::Dates.DateTime
    end_time::Dates.DateTime
    path::Union{Nothing, String}
    format::String
end

function SimulationResultsExport(
    models::Vector{ProblemResultsExport},
    params::SimulationStoreParams;
    start_time=nothing,
    end_time=nothing,
    path=nothing,
    format="csv",
)
    # This end time is outside the bounds of the simulation.
    sim_end_time = params.initial_time + params.step_resolution * params.num_steps

    if start_time === nothing
        start_time = params.initial_time
    elseif start_time < params.initial_time || start_time >= sim_end_time
        throw(IS.InvalidValue("invalid start_time: $start_time"))
    end

    if end_time === nothing
        # Reduce the end_time to be within the simulation.
        end_time = sim_end_time - Dates.Second(1)
    elseif end_time < params.initial_time || end_time >= sim_end_time
        throw(IS.InvalidValue("invalid end_time: $end_time"))
    end

    if !(format in list_supported_formats(SimulationResultsExport))
        throw(IS.InvalidValue("format = $format is not supported"))
    end

    return SimulationResultsExport(
        Dict(x.name => x for x in models),
        start_time,
        end_time,
        path,
        format,
    )
end

function SimulationResultsExport(filename::AbstractString, params::SimulationStoreParams)
    if splitext(filename)[2] != ".json"
        throw(IS.InvalidValue("only JSON files are supported: $filename"))
    end

    return SimulationResultsExport(read_json(filename), params)
end

function SimulationResultsExport(data::AbstractDict, params::SimulationStoreParams)
    models = Vector{ProblemResultsExport}()
    for model in get(data, "models", [])
        if !haskey(model, "name")
            throw(IS.InvalidValue("model data does not define 'name'"))
        end

        problem_params = params.decision_models_params[Symbol(model["name"])]
        duals = Set(
            deserialize_key(problem_params, x) for
            x in get(model, "duals", Set{ConstraintKey}())
        )
        parameters = Set(
            deserialize_key(problem_params, x) for
            x in get(model, "parameters", Set{ParameterKey}())
        )
        variables = Set(
            deserialize_key(problem_params, x) for
            x in get(model, "variables", Set{VariableKey}())
        )
        aux_variables = Set(
            deserialize_key(problem_params, x) for
            x in get(model, "variables", Set{AuxVarKey}())
        )
        problem_export = ProblemResultsExport(
            model["name"],
            duals=duals,
            parameters=parameters,
            variables=variables,
            optimizer_stats=get(model, "optimizer_stats", false),
            store_all_duals=get(model, "store_all_duals", false),
            store_all_parameters=get(model, "store_all_parameters", false),
            store_all_variables=get(model, "store_all_variables", false),
            store_all_aux_variables=get(model, "store_all_aux_variables", false),
        )
        push!(models, problem_export)
    end

    start_time = get(data, "start_time", nothing)
    if start_time isa AbstractString
        start_time = Dates.DateTime(start_time)
    end

    end_time = get(data, "end_time", nothing)
    if end_time isa AbstractString
        end_time = Dates.DateTime(end_time)
    end

    return SimulationResultsExport(
        models,
        params;
        start_time=start_time,
        end_time=end_time,
        path=get(data, "path", nothing),
        format=get(data, "format", "csv"),
    )
end

function get_problem_exports(x::SimulationResultsExport, model_name)
    name = Symbol(model_name)
    if !haskey(x.models, name)
        throw(IS.InvalidValue("model $name is not stored. keys = $(keys(x.models))"))
    end

    return x.models[name]
end

function get_export_file_type(exports::SimulationResultsExport)
    if exports.format == "csv"
        return CSV.File
    end

    throw(IS.InvalidValue("format not supported: $(exports.format)"))
end

list_supported_formats(::Type{SimulationResultsExport}) = ("csv",)

function should_export(exports::SimulationResultsExport, tstamp::Dates.DateTime)
    return tstamp >= exports.start_time && tstamp <= exports.end_time
end

function should_export_dual(exports::SimulationResultsExport, tstamp, model, name)
    return _should_export(exports, tstamp, model, STORE_CONTAINER_DUALS, name)
end

function should_export_parameter(exports::SimulationResultsExport, tstamp, model, name)
    return _should_export(exports, tstamp, model, STORE_CONTAINER_PARAMETERS, name)
end

function should_export_variable(exports::SimulationResultsExport, tstamp, model, name)
    return _should_export(exports, tstamp, model, STORE_CONTAINER_VARIABLES, name)
end

function should_export_expression(exports::SimulationResultsExport, tstamp, model, name)
    return _should_export(exports, tstamp, model, STORE_CONTAINER_EXPRESSIONS, name)
end

function should_export_aux_variable(exports::SimulationResultsExport, tstamp, model, name)
    return _should_export(exports, tstamp, model, STORE_CONTAINER_AUX_VARIABLES, name)
end

function _should_export(exports::SimulationResultsExport, tstamp, model, field_name, name)
    if tstamp < exports.start_time || tstamp >= exports.end_time
        return false
    end

    problem_exports = get_problem_exports(exports, model)
    return _should_export(problem_exports, field_name, name)
end


text:
module BayesianLinearRegressors

using AbstractGPs
using LinearAlgebra
using Random
using Statistics
using PDMats

using AbstractGPs: AbstractGP, _cholesky, FiniteGP

include("bayesian_linear_regression.jl")
include("basis_function_regression.jl")
include("sampling_functions.jl")

export logpdf, rand, mean, std, cov, BayesianLinearRegressor, marginals, posterior
export BasisFunctionRegressor

end # module


text:
__precompile__(false)

module TimeManager

export programStartTime
export currentTime
export programTime
export programTimeStr
export OnTime
export GetTimer
export SetTimer
export UpdateTimers
export GetTimePassed

using Dates
using Printf

TIMERS = Dict{Any,Base.RefValue{Float64}}()
time_start = time()

currentTime(startTime::Real) = (time() - startTime)
programStartTime() = time_start
programTime() = currentTime(time_start)
programTimeStr() = @sprintf("%.3f", programTime())

"""
TODO
"""
GetTimer(key) = Base.getindex(TIMERS[key])

"""
TODO
"""
SetTimer(key, time::Number) = (TIMERS[key] = Ref{Float64}(time))

"""
TODO
"""
function UpdateTimers()
  global TIMERS
  time = GetTimer("FRAME_TIMER")

  LOCAL_TIMERS = filter(tuple -> (f=first(tuple); isa(f,AbstractFloat) || isa(f,Integer)), collect(TIMERS))
  for (milisec,prevTime) in LOCAL_TIMERS
    r=(time - Base.getindex(prevTime)) >= milisec
    if r Base.setindex!(prevTime, time) end
  end
  SetTimer("FRAME_TIMER", Dates.time())
end

timePassed = 0
GetTimePassed() = timePassed
SetTimePassed(t) = global timePassed=t

"""
TODO
"""
function OnTime(milisec::Number)
  global TIMERS

  if !haskey(TIMERS,milisec)
    prevTime = 0.0
    SetTimer(milisec, 0.0)
  else
    prevTime = GetTimer(milisec)
  end

  SetTimePassed(GetTimer("FRAME_TIMER") - prevTime) >= milisec
end

"""
TODO
"""
function OnTime(milisec::Number, prevTime::Ref{Float64}; time=time()) #RefValue
  r = (time - Base.getindex(prevTime)) >= milisec
  if r Base.setindex!(prevTime,time) end
  r
end

SetTimer("FRAME_TIMER", time())

end #TimeManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module TikzPictures

export TikzPicture, PDF, TEX, SVG, save, tikzDeleteIntermediate, tikzCommand, TikzDocument, push!
import Base: push!
import LaTeXStrings: LaTeXString, @L_str, @L_mstr
export LaTeXString, @L_str, @L_mstr

using Compat

_tikzDeleteIntermediate = true
_tikzCommand = "lualatex"

# standalone workaround:
# see http://tex.stackexchange.com/questions/315025/lualatex-texlive-2016-standalone-undefined-control-sequence
_standaloneWorkaround = false

function standaloneWorkaround()
    global _standaloneWorkaround
    _standaloneWorkaround
end

function standaloneWorkaround(value::Bool)
    global _standaloneWorkaround
    _standaloneWorkaround = value
    nothing
end

function tikzDeleteIntermediate(value::Bool)
    global _tikzDeleteIntermediate
    _tikzDeleteIntermediate = value
    nothing
end

function tikzDeleteIntermediate()
    global _tikzDeleteIntermediate
    _tikzDeleteIntermediate
end

function tikzCommand(value::AbstractString)
    global _tikzCommand
    _tikzCommand = value
    nothing
end

function tikzCommand()
    global _tikzCommand
    _tikzCommand
end

type TikzPicture
    data::AbstractString
    options::AbstractString
    preamble::AbstractString
    usePDF2SVG::Bool
    enableWrite18::Bool
    TikzPicture(data::AbstractString; options="", preamble="", usePDF2SVG=true, enableWrite18=true) = new(data, options, preamble, usePDF2SVG, enableWrite18)
end

type TikzDocument
    pictures::Vector{TikzPicture}
    captions::Vector{AbstractString}
end

TikzDocument() = TikzDocument(TikzPicture[], ASCIIString[])

function push!(td::TikzDocument, tp::TikzPicture; caption="")
    push!(td.pictures, tp)
    push!(td.captions, caption)
end

function removeExtension(filename::AbstractString, extension::AbstractString)
    if endswith(filename, extension) || endswith(filename, uppercase(extension))
        return filename[1:(end - length(extension))]
    else
        return filename
    end
end

type PDF
    filename::AbstractString
    PDF(filename::AbstractString) = new(removeExtension(filename, ".pdf"))
end

type TEX
    filename::AbstractString
    include_preamble::Bool
    TEX(filename::AbstractString; include_preamble::Bool=true) = new(removeExtension(filename, ".tex"), include_preamble)
end

type SVG
    filename::AbstractString
    SVG(filename::AbstractString) = new(removeExtension(filename, ".svg"))
end

Base.mimewritable(::MIME"image/svg+xml", tp::TikzPicture) = true

function save(f::TEX, tp::TikzPicture)
    filename = f.filename
    tex = open("$(filename).tex", "w")
    if f.include_preamble
        if standaloneWorkaround()
            println(tex, "\\RequirePackage{luatex85}")
        end
        println(tex, "\\documentclass[tikz]{standalone}")
        println(tex, tp.preamble)
        println(tex, "\\begin{document}")
    end
    print(tex, "\\begin{tikzpicture}[")
    print(tex, tp.options)
    println(tex, "]")
    println(tex, tp.data)
    println(tex, "\\end{tikzpicture}")
    if f.include_preamble
        println(tex, "\\end{document}")
    end
    close(tex)
end

function save(f::TEX, td::TikzDocument)
    if isempty(td.pictures)
        error("TikzDocument does not contain pictures")
    end
    filename = f.filename
    tex = open("$(filename).tex", "w")
    if f.include_preamble
        println(tex, "\\documentclass{article}")
        println(tex, "\\usepackage{caption}")
        println(tex, "\\usepackage{tikz}")
        println(tex, td.pictures[1].preamble)
        println(tex, "\\begin{document}")
    end
    println(tex, "\\centering")
    @assert length(td.pictures) == length(td.captions)
    i = 1
    for tp in td.pictures
        println(tex, "\\centering")
        print(tex, "\\begin{tikzpicture}[")
        print(tex, tp.options)
        println(tex, "]")
        println(tex, tp.data)
        println(tex, "\\end{tikzpicture}")
        print(tex, "\\captionof{figure}{")
        print(tex, td.captions[i])
        println(tex, "}")
        println(tex, "\\vspace{5ex}")
        println(tex)
        i += 1
    end
    if f.include_preamble
        println(tex, "\\end{document}")
    end
    close(tex)
end

function latexerrormsg(s)
    beginError = false
    for l in split(s, '\n')
        if beginError
            if !isempty(l) && l[1] == '?'
                return
            else
                println(l)
            end
        else
            if !isempty(l) && l[1] == '!'
                println(l)
                beginError = true
            end
        end
    end
end

function save(f::PDF, tp::TikzPicture)
    # Generate the .tex file and make pass along any possible errors
    foldername = dirname(f.filename)
    if isempty(foldername)
        foldername = "."
    end
    save(TEX(f.filename * ".tex"), tp)        # Save the tex file in the directory that was given
    #   This will throw an error if the directory doesn't exist

    # From the .tex file, generate a pdf within the specified folder
    latexCommand = ``
    if tp.enableWrite18
        latexCommand = `$(tikzCommand()) --enable-write18 --output-directory=$(foldername) $(f.filename)`
    else
        latexCommand = `$(tikzCommand()) --output-directory=$(foldername) $(f.filename)`
    end
    latexSuccess = success(latexCommand)

    if !latexSuccess
        s = readstring("$(f.filename).log")
        if !standaloneWorkaround() && contains(s, "\\sa@placebox ->\\newpage \\global \\pdfpagewidth")
            standaloneWorkaround(true)
            save(f, tp)
            return
        end
        latexerrormsg(s)
        error("LaTeX error")
    end

    try
        # Shouldn't need to be try-catched anymore, but best to be safe
        # This failing is NOT critical either, so just make it a warning
        if tikzDeleteIntermediate()
            rm("$(f.filename).tex")
            rm("$(f.filename).aux")
            rm("$(f.filename).log")
        end
    catch
        println("WARNING! Your intermediate files are not being deleted.")
    end
end

function save(f::PDF, td::TikzDocument)
    foldername = dirname(f.filename)
    if isempty(foldername)
        foldername = "."
    end
    if isempty(td.pictures)
        error("TikzDocument does not contain pictures")
    end
    try
        save(TEX(f.filename * ".tex"), td)
        if td.pictures[1].enableWrite18
            success(`$(tikzCommand()) --enable-write18 --output-directory=$(foldername) $(f.filename)`)
        else
            success(`$(tikzCommand()) --output-directory=$(foldername) $(f.filename)`)
        end

        if tikzDeleteIntermediate()
            rm("$(f.filename).tex")
            rm("$(f.filename).aux")
            rm("$(f.filename).log")
        end
    catch
        println("Error saving as PDF.")
        rethrow()
    end
end

function save(f::SVG, tp::TikzPicture)
    try
        filename = f.filename
        if tp.usePDF2SVG
            save(PDF(filename), tp)
            success(`pdf2svg $filename.pdf $filename.svg`) || error("pdf2svg failure")
            if tikzDeleteIntermediate()
                rm("$filename.pdf")
            end
        else
            save(TEX("$filename.tex"), tp)
            if tp.enableWrite18
                success(`$(tikzCommand()) --enable-write18 --output-format=dvi $filename`)
            else
                success(`$(tikzCommand()) --output-format=dvi $filename`)
            end
            success(`dvisvgm --no-fonts $filename`)
            if tikzDeleteIntermediate()
                rm("$filename.tex")
                rm("$filename.aux")
                rm("$filename.dvi")
                rm("$filename.log")
            end
        end
    catch
        println("Error saving as SVG")
        rethrow()
    end
end

# this is needed to work with multiple images in ijulia (kind of a hack)
global _tikzid = round(UInt64, time() * 1e6)

@compat function Base.show(f::IO, ::MIME"image/svg+xml", tp::TikzPicture)
    global _tikzid
    filename = "tikzpicture"
    save(SVG(filename), tp)
    s = readstring("$filename.svg")
    s = replace(s, "glyph", "glyph-$(_tikzid)-")
    s = replace(s, "\"clip", "\"clip-$(_tikzid)-")
    s = replace(s, "#clip", "#clip-$(_tikzid)-")
    s = replace(s, "\"image", "\"image-$(_tikzid)-")
    s = replace(s, "#image", "#image-$(_tikzid)-")
    s = replace(s, "linearGradient id=\"linear", "linearGradient id=\"linear-$(_tikzid)-")
    s = replace(s, "#linear", "#linear-$(_tikzid)-")
    s = replace(s, "image id=\"", "image style=\"image-rendering: pixelated;\" id=\"")
    _tikzid += 1
    println(f, s)
    if tikzDeleteIntermediate()
        rm("$filename.svg")
    end
end

end # module


text:
using LinearAlgebra
using Yao
using Yao.Blocks
using QuAlgorithmZoo
using Statistics

#=
"""wavefunction through generator"""
function psiGen(qg::QuGAN)
    regGen = copy(qg.reg0)
    regGenStore = zero_state(0)
    for i=1:2qg.dG ##each layer consists of 1 layer of rotation gates and 1 layer of cnot gates.
        regGen |> qg.generator[i]
        if i % 2qg.dGB == 0
            regGenStore = join(regGenStore, copy(focus!(regGen, 1:qg.RBit)))
            relax!(regGen, 1:qg.RBit)
            measure_reset!(regGen, 1, val = 0)
        end
    end
    regGenStore = join(regGenStore, copy(focus!(regGen, (qg.RBit+1):qg.nBitGB))
    regGenStore
end

"""Target qubit through discriminator"""
function psiTarDis(qg::QuGAN)
    regGenStore = psiGen(qg)
    regDis = join(copy(focus!(copy(regGenStore),(qg.RBit-qg.VBit-1+1):qg.nBitGB)),zero_state(1))
    regDisStore = zero_state(0)
    for i=1:(2qg.dD+1)
        regDis |> qg.discriminator[i]
        if i % 2qg.dDB == 0
            regDisStore = join(regDisStore, copy(focus!(regDis, 1:qg.RBit)))
            RBitDis = focus!(copy(regGenStore),(qg.RBit-qg.VBit-1+1))
            measure_reset!(regBit, 1, val = measure!(RBitDis))
        end
    end
    regDisStore = join(regDisStore, copy(focus!(regDis, (qg.RBit+1):(qg.nBitGB+1)))
    regDisStore
end
=#

function expectCB(reg::DefaultRegister)
    dim = 1<<nqubits(reg)
    CpOp = diagm(0=>[0:1:(dim-1);])
    expect(matrixgate(CpOp), reg) |>abs
end

function distance(qg::QuMGAN)
    regGenAll = psiGen(qg, qg.reg0)
    GenP1 = regGenAll[1]
    if (qg.N-qg.nBitGB) != (qg.RBit*(qg.Nblock-1))
        println("Dim mismatch!!")
    end
    GenP2 = product_state(ComplexF64, (qg.N-qg.nBitGB), (regGenAll[2] |> round |> Int64) )
    Gen = join(GenP2, GenP1)
    tracedist(qg.target, Gen)[]
end


text:
import Base.Test: @testset, @test, @test_throws

using EffectiveWaves

@testset "Summary" begin

include("types_constructors.jl")

include("large_vol_low_freq_effective.jl")
include("integrated_reflection.jl")

end


text:
using EngEconomics

# Given
aInit = 1700
bInit = 2100
cInit = 3750
aAnnualBenefit = 1000
bAnnualBenefit = 1000
cAnnualBenefit = 1000
aLife = 2
bLife = 3
cLife = 6
i = 0.08

# Find: Present Worth Value of each Item then pick the best option
# Start by determining the number of periods for each product
years = 6
aPeriods = years / aLife
bPeriods = years / bLife
cPeriods = years / cLife

# Next, determine the present worth of product A
aInitPW = -aInit * (1 + presentWorthFactor(i, aLife) + presentWorthFactor(i, 2 * aLife))
aAnnualBenefitPW = (aAnnualBenefit * seriesPresentAmountFactor(i, aLife)) * (1 + presentWorthFactor(i, aLife) + presentWorthFactor(i, 2 * aLife))
aPW = aInitPW + aAnnualBenefitPW

# Determine PW of Product B
bInitPW = -bInit * (1 + presentWorthFactor(i, bLife))
bAnnualBenefitPW = (bAnnualBenefit * seriesPresentAmountFactor(i, bLife)) * (1 + presentWorthFactor(i, bLife))
bPW = bInitPW + bAnnualBenefitPW

# Determine PW of Product C
cInitPW = -cInit * (1)
cAnnualBenefitPW = (cAnnualBenefit * seriesPresentAmountFactor(i, cLife)) * (1)
cPW = cInitPW + cAnnualBenefitPW

# Compare Alternatives, pick the greatest one
maxPW = max(aPW, bPW, cPW)

if maxPW == aPW
	println("Pick A")
elseif maxPW == bPW
	println("Pick B")
else
	println("Pick C")
end


text:
"""
Run gibbs sampler.
"""
function gibbs_sample!(
    model::DistributedSeqModel,
    spikes::Vector{Spike},
    initial_assignments::Vector{Int64},
    num_samples::Int64,
    extra_split_merge_moves::Int64,
    split_merge_window::Float64,
    save_every::Int64;
    verbose=false
)

    if extra_split_merge_moves > 0
        @warn "Split merge not implemented for distributed model."
    end

    num_partitions = model.num_partitions
    max_time = model.primary_model.max_time

    # Partition spikes.
    split_points = Tuple(p * max_time / num_partitions for p in 0:num_partitions)
    spk_partition = Tuple(Spike[] for m in model.submodels)
    assgn_partition = [Int64[] for m in model.submodels]  # TODO -- Tuple?
    partition_ids = [Int64[] for m in model.submodels]    # TODO -- Tuple?

    for s = 1:length(spikes)
        for p = 1:model.num_partitions
            if split_points[p] <= spikes[s].timestamp <= split_points[p + 1]
                push!(spk_partition[p], spikes[s])
                push!(assgn_partition[p], initial_assignments[s])
                push!(partition_ids[p], s)
                break
            end
        end
    end

    # Dense rank assignments within each partition.
    for p = 1:model.num_partitions
        idx = assgn_partition[p] .> 0 # ignore background spikes.
        assgn_partition[p][idx] .= denserank(assgn_partition[p][idx])
    end

    # Pass assignments to submodels
    for p in 1:model.num_partitions
        recompute!(
            model.submodels[p],
            spk_partition[p],
            assgn_partition[p],
        )
    end

    # Save spike assignments over samples.
    assignments = initial_assignments
    collect_assignments!(model, assignments, assgn_partition, partition_ids) # updates assignments.

    # Order to iterate over spikes.
    spike_order_partition = Tuple(
        collect(1:length(subspikes)) for subspikes in spk_partition
    )

    # ======== THINGS TO SAVE ======== #

    # number of samples to save.
    n_saved_samples = Int(round(num_samples / save_every))

    # log likelihoods.
    log_p_hist = zeros(n_saved_samples)

    # parent assignments.
    assignment_hist = zeros(
        Int64, length(spikes), n_saved_samples
    )

    # latent event assignment ids, times, types, amplitudes.
    latent_event_hist = Vector{EventSummaryInfo}[]

    # global variables (offsets, response amplitudes, etc.)
    globals_hist = SeqGlobals[]

     # ======== MAIN LOOP ======== #

    # Draw samples.
    for s = 1:num_samples

        # Update spike assignments in random order.
        Threads.@threads for p in 1:model.num_partitions
            Random.shuffle!(spike_order_partition[p])
            for i in spike_order_partition[p]
                remove_datapoint!(
                    model.submodels[p],
                    spk_partition[p][i],
                    assgn_partition[p][i]
                )
                assgn_partition[p][i] = gibbs_add_datapoint!(
                    model.submodels[p],
                    spk_partition[p][i]
                )
            end
        end

        # Pass spike assignment updates and latent event updates to primary model.
        # Latent events.
        collect_assignments!(model, assignments, assgn_partition, partition_ids)

        # Update latent events.
        gibbs_update_latents!(model.primary_model)

        # Update globals
        gibbs_update_globals!(model.primary_model, spikes, assignments)

        # No need to pass updates to submodels, since the globals
        # and events point to the same objects.

        # Store results
        if (s % save_every) == 0
            j = Int(s / save_every)

            # Save log likelihood
            log_p_hist[j] = log_like(model, spikes)

            # Save assignments.
            assignment_hist[:, j] .= assignments

            # Save latent event information.
            push!(
                latent_event_hist,
                event_list_summary(model.primary_model)
            )

            # Save global variables.
            push!(globals_hist, deepcopy(model.primary_model.globals))

            verbose && print(s, "-")
        end
    end
    verbose && println("Done")

    return (
        assignments,
        assignment_hist,
        log_p_hist,
        latent_event_hist,
        globals_hist
    )
end

"""
Collects assignments from submodels into `assignments` vector
and populates `model.primary_model.sequence_events` with
submodel's events.
"""
function collect_assignments!(
        model::DistributedSeqModel,
        assignments::Vector{Int64},
        assgn_partition::Vector{Vector{Int64}},
        partition_ids::Vector{Vector{Int64}}
    )

    empty!(model.primary_model.sequence_events.events)
    empty!(model.primary_model.sequence_events.indices)
    event_id = 1
    event_map = Dict{Tuple{Int64,Int64},Int64}()

    # Create flat assignment index for each event, to be used
    # in `primary_model`.
    for p = 1:model.num_partitions

        # Iterate over event indices in submodel-p.
        for ind in model.submodels[p].sequence_events.indices
            event_map[(p, ind)] = event_id
            push!(
                model.primary_model.sequence_events.events,
                model.submodels[p].sequence_events[ind]
            )
            push!(
                model.primary_model.sequence_events.indices,
                event_id
            )
            event_id += 1
        end

        # Assign each spike in `primary_model` using the
        # event ids generated above.
        for (s, assgn) in enumerate(assgn_partition[p])
            if assgn == -1
                assignments[partition_ids[p][s]] = -1
            else
                assignments[partition_ids[p][s]] = event_map[(p, assgn)]
            end
        end

    end
end


text:
import Pkg; Pkg.activate(".")

using Distributions
using Statistics
using StatsBase

# z-score of the oceanic network

# map values of posterior distribution
betab_posterior = CSV.read(joinpath("data", "posterior_distributions", "beta_binomial_posterior.csv"))
mu_map = median(betab_posterior[:mu])
phi_map = exp(median(betab_posterior[:phi]))

# Empirical number of species and links
Sx = 11367
Lx = 7062647

means = (Sx^2 - Sx + 1) * mu_map + Sx - 1
sds = sqrt((Sx^2 - Sx + 1) * mu_map * (1 - mu_map) * (1 + Sx* (Sx - 1) / (1 + phi_map)))

z_score = (Lx - means) / sds


text:
@testset "Storage data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type GenericBattery, consider changing the device models"
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5 = build_system("c_sys5")
    c_sys14 = build_system("c_sys14")
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Storage, model)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Storage, model)
end

@testset "Storage Basic Storage With DC - PF" begin
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = build_system("c_sys5_bat")
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 72, 0, 72, 72, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage Basic Storage With AC - PF" begin
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = build_system("c_sys5_bat")
    op_problem = OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 96, 0, 96, 96, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation DC - PF" begin
    model = DeviceModel(GenericBattery, BookKeepingwReservation)
    c_sys5_bat = build_system("c_sys5_bat")
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 96, 0, 72, 72, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation With AC - PF" begin
    model = DeviceModel(GenericBattery, BookKeepingwReservation)
    c_sys5_bat = build_system("c_sys5_bat")
    op_problem = OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 120, 0, 96, 96, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end


text:
# This file was generated, do not modify it. # hide
using MLJ, RDatasets, Random
MLJ.color_off() # hide

data = dataset("datasets", "USArrests")
names(data)

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
#Ellipse methods

Ellipse(x::Int, y::Int, args...; kwargs...) = Ellipse(Point(x,y), args...; kwargs...)
Ellipse(p::CartesianIndex{2}, args...; kwargs...) = Ellipse(Point(p), args...; kwargs...)
Ellipse(circle::CirclePointRadius) = Ellipse(circle.center, circle.ρ, circle.ρ; thickness = circle.thickness, fill = circle.fill)

function draw!(img::AbstractArray{T, 2}, ellipse::Ellipse, color::T) where T<:Colorant
	ys = Int[]
	xs = Int[]
	break_point = 0
	if ellipse.fill == false
		break_point = ((ellipse.ρy - ellipse.thickness) / ellipse.ρy) ^ 2 + ((ellipse.ρx - ellipse.thickness) / ellipse.ρx) ^ 2
	end
	for i in ellipse.center.y : ellipse.center.y + ellipse.ρy
		for j in ellipse.center.x : ellipse.center.x + ellipse.ρx
			val = ((i - ellipse.center.y) / ellipse.ρy) ^ 2 + ((j - ellipse.center.x) / ellipse.ρx) ^ 2
			if val < 1 && val >= break_point
				push!(ys, i)
				push!(xs, j)
			end
		end
	end
	for (yi, xi) in zip(ys, xs)
		drawifinbounds!(img, yi, xi, color)
		drawifinbounds!(img, 2 * ellipse.center.y - yi, xi, color)
		drawifinbounds!(img, yi, 2 * ellipse.center.x - xi, color)
		drawifinbounds!(img, 2 * ellipse.center.y - yi, 2 * ellipse.center.x - xi, color)
	end
	img
end


text:
# wrappers of the low-level CUBLAS functionality

function cufftGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  cufftGetProperty(value_ref)
  value_ref[]
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export cyto_cytoscape

"""
    cyto_cytoscape(;kwargs...)

A Cytoscape component.
A Component Library for Dash aimed at facilitating network visualization in
Python, wrapped around [Cytoscape.js](http://js.cytoscape.org/).
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `autoRefreshLayout` (Bool; optional): Whether the layout should be refreshed when elements are added or removed.
- `autolock` (Bool; optional): Whether nodes should be locked (not draggable at all) by default
(if true, overrides individual node state).
- `autoungrabify` (Bool; optional): Whether nodes should be ungrabified (not grabbable by user) by
default (if true, overrides individual node state).
- `autounselectify` (Bool; optional): Whether nodes should be unselectified (immutable selection state) by
default (if true, overrides individual element state).
- `boxSelectionEnabled` (Bool; optional): Whether box selection (i.e. drag a box overlay around, and release it
to select) is enabled. If enabled, the user must taphold to pan the graph.
- `className` (String; optional): Sets the class name of the element (the value of an element's html
class attribute).
- `contextmenu` (optional): Displays a context menu on right click. Requires extra layouts loaded.
Accepts a list of dictionaries, each of which describes a context
menu option. Options are rendered in the order presented.. contextmenu has the following type: Array of lists containing elements 'id', 'selector', 'content', 'tooltipText', 'disabled'.
Those elements have the following types:
  - `id` (String; optional): ID associated with option.
  - `selector` (String; optional): Determines which Cytoscape elements the option is attached to. Takes in a Cytoscape selector
(see Cytoscape documentation for more information). Examples of valid selectors include node,
edge, and core.
  - `content` (String; optional): Label assigned to option.
  - `tooltipText` (String; optional): Hover tooltip text assigned to option.
  - `disabled` (Bool; optional): Toggles option disabled (greyed out).s
- `contextmenuData` (optional): Dictionary returned when a context menu option is selected. Read-only.. contextmenuData has the following type: lists containing elements 'id', 'position', 'timestamp', 'target', 'coordinates'.
Those elements have the following types:
  - `id` (String; optional): ID associated with option selected.
  - `position` (optional): Position associated with option selected.. position has the following type: lists containing elements 'x', 'y'.
Those elements have the following types:
  - `x` (Real; optional)
  - `y` (Real; optional)
  - `timestamp` (Real; optional): Time the option was selected.
  - `target` (Dict; optional): Dictionary containing information about the selected item. Information provided varies depending
on the type of the selected item (node, edge, core, etc.).
  - `coordinates` (Array of Reals; optional): Array containing latitude and longitude where context menu was opened if leaflet is enabled.
- `elements` (optional): A list of dictionaries representing the elements of the networks. Each dictionary describes an element, and
specifies its purpose. The [official Cytoscape.js documentation](https://js.cytoscape.org/#notation/elements-json)
offers an extensive overview and examples of element declaration.
Alternatively, a dictionary with the format { 'nodes': [], 'edges': [] } is allowed at initialization,
but arrays remain the recommended format.. elements has the following type: Array of lists containing elements 'group', 'data', 'position', 'selected', 'selectable', 'locked', 'grabbable', 'classes'.
Those elements have the following types:
  - `group` (String; optional): Either 'nodes' or 'edges'. If not given, it's automatically inferred.
  - `data` (optional): Element specific data.. data has the following type: lists containing elements 'id', 'label', 'parent', 'source', 'target'.
Those elements have the following types:
  - `id` (String; optional): Reference to the element, useful for selectors and edges. Randomly assigned if not given.
  - `label` (String; optional): Optional name for the element, useful when `data(label)` is given to a style's `content`
or `label`. It is only a convention.
  - `parent` (String; optional): Only for nodes. Optional reference to another node. Needed to create compound nodes.
  - `source` (String; optional): Only for edges. The id of the source node, which is where the edge starts.
  - `target` (String; optional): Only for edges. The id of the target node, where the edge ends.
  - `position` (optional): Only for nodes. The position of the node.. position has the following type: lists containing elements 'x', 'y'.
Those elements have the following types:
  - `x` (Real; optional): The x-coordinate of the node.
  - `y` (Real; optional): The y-coordinate of the node.
  - `selected` (Bool; optional): If the element is selected upon initialisation.
  - `selectable` (Bool; optional): If the element can be selected.
  - `locked` (Bool; optional): Only for nodes. If the position is immutable.
  - `grabbable` (Bool; optional): Only for nodes. If the node can be grabbed and moved by the user.
  - `classes` (String; optional): Space separated string of class names of the element. Those classes can be selected
by a style selector.s | lists containing elements 'nodes', 'edges'.
Those elements have the following types:
  - `nodes` (Array; optional)
  - `edges` (Array; optional)
- `generateImage` (optional): Dictionary specifying options to generate an image of the current cytoscape graph.
Value is cleared after data is received and image is generated. This property will
be ignored on the initial creation of the cytoscape object and must be invoked through
a callback after it has been rendered.

If the app does not need the image data server side and/or it will only be used to download
the image, it may be prudent to invoke `'download'` for `action` instead of
`'store'` to improve performance by preventing transfer of data to the server.. generateImage has the following type: lists containing elements 'type', 'options', 'action', 'filename'.
Those elements have the following types:
  - `type` (a value equal to: 'svg', 'png', 'jpg', 'jpeg'; optional): File type to output
  - `options` (Dict; optional): Dictionary of options to cy.png() / cy.jpg() or cy.svg() for image generation.
See https://js.cytoscape.org/#core/export for details. For `'output'`, only 'base64'
and 'base64uri' are supported. Default: `{'output': 'base64uri'}`.
  - `action` (a value equal to: 'store', 'download', 'both'; optional): `'store'`: Stores the image data (only jpg and png are supported)
in `imageData` and invokes server-side Dash callbacks. `'download'`: Downloads the image
as a file with all data handling done client-side. No `imageData` callbacks are fired.
`'both'`: Stores image data and downloads image as file. The default is `'store'`
  - `filename` (String; optional): Name for the file to be downloaded. Default: 'cyto'.
- `imageData` (String; optional): String representation of the image requested with generateImage. Null if no
image was requested yet or the previous request failed. Read-only.
- `layout` (optional): A dictionary specifying how to set the position of the elements in your
graph. The `'name'` key is required, and indicates which layout (algorithm) to
use. The keys accepted by `layout` vary depending on the algorithm, but these
keys are accepted by all layouts: `fit`,  `padding`, `animate`, `animationDuration`,
`boundingBox`.

 The complete list of layouts and their accepted options are available on the
 [Cytoscape.js docs](https://js.cytoscape.org/#layouts) . For the external layouts,
the options are listed in the "API" section of the  README.
 Note that certain keys are not supported in Dash since the value is a JavaScript
 function or a callback. Please visit this
[issue](https://github.com/plotly/dash-cytoscape/issues/25) for more information.. layout has the following type: lists containing elements 'name', 'fit', 'padding', 'animate', 'animationDuration', 'boundingBox'.
Those elements have the following types:
  - `name` (a value equal to: 'random', 'preset', 'circle', 'concentric', 'grid', 'breadthfirst', 'cose', 'cose-bilkent', 'fcose', 'cola', 'euler', 'spread', 'dagre', 'klay'; required): The layouts available by default are:
  `random`: Randomly assigns positions.
  `preset`: Assigns position based on the `position` key in element dictionaries.
  `circle`: Single-level circle, with optional radius.
  `concentric`: Multi-level circle, with optional radius.
  `grid`: Square grid, optionally with numbers of `rows` and `cols`.
  `breadthfirst`: Tree structure built using BFS, with optional `roots`.
  `cose`: Force-directed physics simulation.

Some external layouts are also included. To use them, run
  `dash_cytoscape.load_extra_layouts()` before creating your Dash app. Be careful about
  using the extra layouts when not necessary, since they require supplementary bandwidth
  for loading, which impacts the startup time of the app.
  The external layouts are:
  [cose-bilkent](https://github.com/cytoscape/cytoscape.js-cose-bilkent),
  [fcose](https://github.com/iVis-at-Bilkent/cytoscape.js-fcose),
  [cola](https://github.com/cytoscape/cytoscape.js-cola),
  [euler](https://github.com/cytoscape/cytoscape.js-dagre),
  [spread](https://github.com/cytoscape/cytoscape.js-spread),
  [dagre](https://github.com/cytoscape/cytoscape.js-dagre),
  [klay](https://github.com/cytoscape/cytoscape.js-klay),
  - `fit` (Bool; optional): Whether to render the nodes in order to fit the canvas.
  - `padding` (Real; optional): Padding around the sides of the canvas, if fit is enabled.
  - `animate` (Bool; optional): Whether to animate change in position when the layout changes.
  - `animationDuration` (Real; optional): Duration of animation in milliseconds, if enabled.
  - `boundingBox` (Dict; optional): How to constrain the layout in a specific area. Keys accepted are either
`x1, y1, x2, y2` or `x1, y1, w, h`, all of which receive a pixel value.
- `leaflet` (optional): Dictionary specifying configuration options to overlay a leaflet map on top of Cytoscape. All
configuration options are optional; provide an empty dictionary to use default options. Requires
latitude and longitude properties to be included in the node data for positional information.
Requires preset layout to be used. Requires extra layouts to be loaded.. leaflet has the following type: lists containing elements 'provider', 'tileUrl', 'attribution', 'maxZoom', 'latitudeId', 'longitudeId', 'zoomOffset', 'tileSize', 'view'.
Those elements have the following types:
  - `provider` (String; optional): Specify a tile preset from Leaflet providers instead of a manual tileUrl.  See http://leaflet-extras.github.io/leaflet-providers/preview/  See also https://github.com/leaflet-extras/leaflet-providers
  - `tileUrl` (String; optional): Endpoint used by leaflet to fetch map tiles (if not using provider).
  - `attribution` (String; optional): Attribution text displayed on the bottom right corner of the map (if not using provider).
  - `maxZoom` (Real; optional): Sets the max zoom allowed by leaflet (if not using provider). See leaflet documentation for more information about zoom.
  - `latitudeId` (String; optional): Specifies the name of the node property containing the latitude of the node. Default: 'lat'.
  - `longitudeId` (String; optional): Specifies the name of the node property containing the longitude of the node. Default: 'lon'.
  - `zoomOffset` (Real; optional): Sets the offset from the zoom number used in tile URLs
  - `tileSize` (Real; optional): Specifies the size of each tile image retrieved by leaflet
  - `view` (Array of Reals; optional): Specifies a pair of coordinates (latitude, longitude) to centre the viewport over and an optional zoom level for the viewport using an array
- `maxZoom` (Real; optional): A maximum bound on the zoom level of the graph. The viewport can not be
scaled larger than this zoom level.
- `minZoom` (Real; optional): A minimum bound on the zoom level of the graph. The viewport can not be
scaled smaller than this zoom level.
- `mouseoverEdgeData` (Dict; optional): The data dictionary of an edge returned when you hover over it. Read-only.
- `mouseoverNodeData` (Dict; optional): The data dictionary of a node returned when you hover over it. Read-only.
- `pan` (optional): Dictionary indicating the initial panning position of the graph. The
following keys are accepted:. pan has the following type: lists containing elements 'x', 'y'.
Those elements have the following types:
  - `x` (Real; optional): The x-coordinate of the node
  - `y` (Real; optional): The y-coordinate of the node
- `panningEnabled` (Bool; optional): Whether panning the graph is enabled (i.e., the position of the graph is
mutable overall).
- `responsive` (Bool; optional): Toggles intelligent responsive resize of Cytoscape graph with viewport size change
- `selectedEdgeData` (Array; optional): The list of data dictionaries of all selected edges (e.g. using
Shift+Click to select multiple nodes, or Shift+Drag to use box selection). Read-only.
- `selectedNodeData` (Array; optional): The list of data dictionaries of all selected nodes (e.g. using
Shift+Click to select multiple nodes, or Shift+Drag to use box selection). Read-only.
- `style` (Dict; optional): Add inline styles to the root element.
- `stylesheet` (optional): A list of dictionaries representing the styles of the elements.
Each dictionary requires the following keys: `selector` and `style`.

Both the [selector](https://js.cytoscape.org/#selectors) and
the [style](https://js.cytoscape.org/#style/node-body) are
exhaustively documented in the Cytoscape.js docs. Although methods such
as `cy.elements(...)` and `cy.filter(...)` are not available, the selector
string syntax stays the same.. stylesheet has the following type: Array of lists containing elements 'selector', 'style'.
Those elements have the following types:
  - `selector` (String; required): Which elements you are styling. Generally, you select a group of elements (node, edges, both),
a class (that you declare in the element dictionary), or an element by ID.
  - `style` (Dict; required): What aspects of the elements you want to modify. This could be the size or
color of a node, the shape of an edge arrow, or many more.s
- `tapEdge` (optional): The complete edge dictionary returned when you tap or click it. Read-only.. tapEdge has the following type: lists containing elements 'isLoop', 'isSimple', 'midpoint', 'sourceData', 'sourceEndpoint', 'targetData', 'targetEndpoint', 'timeStamp', 'classes', 'data', 'grabbable', 'group', 'locked', 'selectable', 'selected', 'style'.
Those elements have the following types:
  - `isLoop` (Bool; optional): Edge-specific item
  - `isSimple` (Bool; optional): Edge-specific item
  - `midpoint` (Dict; optional): Edge-specific item
  - `sourceData` (Dict; optional): Edge-specific item
  - `sourceEndpoint` (Dict; optional): Edge-specific item
  - `targetData` (Dict; optional): Edge-specific item
  - `targetEndpoint` (Dict; optional): Edge-specific item
  - `timeStamp` (Real; optional): Edge-specific item
  - `classes` (String; optional): General item (for all elements)
  - `data` (Dict; optional): General item (for all elements)
  - `grabbable` (Bool; optional): General item (for all elements)
  - `group` (String; optional): General item (for all elements)
  - `locked` (Bool; optional): General item (for all elements)
  - `selectable` (Bool; optional): General item (for all elements)
  - `selected` (Bool; optional): General item (for all elements)
  - `style` (Dict; optional): General item (for all elements)
- `tapEdgeData` (Dict; optional): The data dictionary of an edge returned when you tap or click it. Read-only.
- `tapNode` (optional): The complete node dictionary returned when you tap or click it. Read-only.. tapNode has the following type: lists containing elements 'edgesData', 'renderedPosition', 'timeStamp', 'classes', 'data', 'grabbable', 'group', 'locked', 'position', 'selectable', 'selected', 'style', 'ancestorsData', 'childrenData', 'descendantsData', 'parentData', 'siblingsData', 'isParent', 'isChildless', 'isChild', 'isOrphan', 'relativePosition'.
Those elements have the following types:
  - `edgesData` (Array; optional): node specific item
  - `renderedPosition` (Dict; optional): node specific item
  - `timeStamp` (Real; optional): node specific item
  - `classes` (String; optional): General item (for all elements)
  - `data` (Dict; optional): General item (for all elements)
  - `grabbable` (Bool; optional): General item (for all elements)
  - `group` (String; optional): General item (for all elements)
  - `locked` (Bool; optional): General item (for all elements)
  - `position` (Dict; optional): General item (for all elements)
  - `selectable` (Bool; optional): General item (for all elements)
  - `selected` (Bool; optional): General item (for all elements)
  - `style` (Dict; optional): General item (for all elements)
  - `ancestorsData` (Dict | Array; optional): Item for compound nodes
  - `childrenData` (Dict | Array; optional): Item for compound nodes
  - `descendantsData` (Dict | Array; optional): Item for compound nodes
  - `parentData` (Dict | Array; optional): Item for compound nodes
  - `siblingsData` (Dict | Array; optional): Item for compound nodes
  - `isParent` (Bool; optional): Item for compound nodes
  - `isChildless` (Bool; optional): Item for compound nodes
  - `isChild` (Bool; optional): Item for compound nodes
  - `isOrphan` (Bool; optional): Item for compound nodes
  - `relativePosition` (Dict; optional): Item for compound nodes
- `tapNodeData` (Dict; optional): The data dictionary of a node returned when you tap or click it. Read-only.
- `userPanningEnabled` (Bool; optional): Whether user events (e.g. dragging the graph background) are allowed to
pan the graph.
- `userZoomingEnabled` (Bool; optional): Whether user events (e.g. dragging the graph background) are allowed
to pan the graph.
- `zoom` (Real; optional): The initial zoom level of the graph. You can set `minZoom` and
`maxZoom` to set restrictions on the zoom level.
- `zoomingEnabled` (Bool; optional): Whether zooming the graph is enabled (i.e., the zoom level of the graph
is mutable overall).
"""
function cyto_cytoscape(; kwargs...)
        available_props = Symbol[:id, :autoRefreshLayout, :autolock, :autoungrabify, :autounselectify, :boxSelectionEnabled, :className, :contextmenu, :contextmenuData, :elements, :generateImage, :imageData, :layout, :leaflet, :maxZoom, :minZoom, :mouseoverEdgeData, :mouseoverNodeData, :pan, :panningEnabled, :responsive, :selectedEdgeData, :selectedNodeData, :style, :stylesheet, :tapEdge, :tapEdgeData, :tapNode, :tapNodeData, :userPanningEnabled, :userZoomingEnabled, :zoom, :zoomingEnabled]
        wild_props = Symbol[]
        return Component("cyto_cytoscape", "Cytoscape", "dash_cytoscape", available_props, wild_props; kwargs...)
end

max_stars_repo_name: zyphr-solutions/dash-cytoscape

max_issues_repo_name: zyphr-solutions/dash-cytoscape

max_forks_repo_name: zyphr-solutions/dash-cytoscape


text:
using ShaleDrillingLikelihood
using SparseArrays
using Random
using Distributions
using Dates
using StatsBase
using LinearAlgebra
using JLD2

using Base.Iterators: product, OneTo

using ShaleDrillingLikelihood: ObservationDrill,
    SimulationDraw,
    idx_ρ,
    actionspace,
    ichars_sample,
    j1chars,
    tptr,
    _x, _y,
    zchars,
    ztrange,
    InitialDrilling,
    _i,
    j1_sample,
    j1_range,
    jtstart,
    exploratory_terminal,
    ssprime,
    total_leases,
    zero_out_small_probs,
    lrvar,
    lrstd,
    lrmean,
    simulate,
    theta_royalty_ρ,
    revenue,
    cost,
    extend,
    theta_ρ,
    vw_revenue,
    vw_cost,
    vw_extend,
    DataDynamicDrill,
    max_states,
    _D,
    simloglik!,
    simloglik_drill_data!,
    SimulationDraws,
    DCDPTmpVars,
    _model,
    theta_drill_ρ,
    update!,
    value_function,
    EV, dEV,
    _x, _y,
    end_ex0, end_ex1, end_lrn, end_inf,
    state_idx,
    j2ptr, tptr, j1ptr, ichars, zchars

Random.seed!(7)

discount = ((0x1.006b55c832502p+0)^12 / 1.125) ^ (1/4)  # real discount rate

# set up coefs
θρ = 0.0
αψ = 0.33
αg = 0.56

# set up coefs
#            θρ  ψ        g    x_r     κ₁   κ₂
θ_royalty = [θρ, 1.0,    0.1, -0.3,  -0.6, 0.6]
#                αψ, αg, γx   σ2η, σ2u   # η is iwt, u is iw
θ_produce = vcat(αψ, αg, 0.2, 0.3, 0.4)
#            drill  ext    α0  αg  αψ  θρ
θ_drill_u = [-5.5, -2.0, -2.8, αg, αψ, θρ]
# θ_drill_u = [-5.8,       -2.8, αg, αψ, θρ]
θ_drill_c = vcat(θ_drill_u[1:2], θρ)

# model
rwrd_c = DrillReward(DrillingRevenue(Constrained(;log_ogip=αg, α_ψ=αψ),NoTrend(),NoTaxes()), DrillingCost_constant(), ExtensionCost_Constant())
rwrd_u = UnconstrainedProblem(rwrd_c)
rwrd = rwrd_u

# parms
num_i = 350

# grid sizes
nψ =  13
nz = 15

# simulations
M = 250

# observations
num_zt = 150          # drilling
obs_per_well = 10:20  # pdxn

# roylaty
royalty_rates = [1/8, 3/16, 1/4,]
num_royalty_rates = length(royalty_rates)
nk_royalty = length(θ_royalty) - (length(royalty_rates)-1)-2

# geology
ogip_dist = Normal(4.68, 0.31)

# price parameters
zrho = 0.8                # rho
zmean = 1.33*(1-zrho)     # mu
zvar = 0.265^2*(1-zrho^2) # sigsq

# state space, psi-space
wp = LeasedProblem(8, 8, 12, 12, 8)
ψs = range(-4.5; stop=4.5, length=nψ)

# simulate price process
zprocess = AR1process(zmean, zrho, zvar)
zvec = simulate(zprocess, num_zt)
ztrng = range(Date(2003,10); step=Quarter(1), length=num_zt)
_zchars = ExogTimeVars(tuple.(zvec), ztrng)

# price grid2
zmin = min(minimum(zvec), lrmean(zprocess)-3*lrstd(zprocess))
zmax = max(maximum(zvec), lrmean(zprocess)+3*lrstd(zprocess))
zrng = range(zmin, zmax; length=nz)
zs = tuple(zrng)
ztrans = sparse(zero_out_small_probs(tauchen_1d(zprocess, zrng), 1e-4))
@show zrng

# random vars
log_ogip = rand(ogip_dist, num_i)
Xroyalty = vcat(log_ogip', randn(nk_royalty-1, num_i))
u,v = randn(num_i), randn(num_i)

# model
ddm_with_t1ev = DynamicDrillModel(rwrd_u, discount, wp, zs, ztrans, ψs, true)

# construct royalty data
data_roy = DataRoyalty(u, v, Xroyalty, θ_royalty, num_royalty_rates)
royrates_endog = royalty_rates[_y(data_roy)]
royrates_exog = sample(royalty_rates, num_i)

# construct ichars for drilling
_ichars = [gr for gr in zip(log_ogip, royrates_exog)]

ddm_opts = (minmaxleases=0:0, nper_initial=0:0, nper_development=60:60, tstart=1:75, xdomain=0:0)
data = DataDrill(u, v, _zchars, _ichars, ddm_with_t1ev, θ_drill_u; ddm_opts...)

modulepath = dirname(pathof(ShaleDrillingLikelihood))
filepath = joinpath(modulepath, "..", "scratch", "simulated-data.jld2")

jldopen(filepath, "w") do file
    file["j1ptr"]   = j1ptr(data)
    file["j2ptr"]   = j2ptr(data)
    file["tptr"]    = tptr(data)
    file["jtstart"] = jtstart(data)
    file["j1chars"] = j1chars(data)
    file["ichars"]  = ichars(data)
    file["y"]       = _y(data)
    file["x"]       = _x(data)
    file["zchars"]  = zchars(data).timevars
    file["zrng"] = zrng
    file["psis"] = ψs
    file["ztrans"] = ztrans
    file["theta"] = θ_drill_u
    file["discount"] = discount
end


text:
"""
    Exo1Quaternion{T <: Real} = Parabolic{Quaternion{T}}

An 1-exo-quaternion is a parabolic Cayley-Dickson construct with quaternions.
"""
const Exo1Quaternion{T <: Real} = Parabolic{Quaternion{T}}

function Exo1Quaternion(z::Quaternion{T}) where T <: Real
    Exo1Quaternion{T}(z, zero(Quaternion{T}))
end

function Exo1Quaternion(z::Binion{T}) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(z, zero(Binion{T})),
        zero(Quaternion{T})
    )
end

function Exo1Quaternion(a::T, b::T, c::T, d::T, f::T, g::T, h::T, j::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Quaternion{T}(
            Binion{T}(f, g),
            Binion{T}(h, j)
        )
    )
end

function Exo1Quaternion(a::T, b::T, c::T, d::T, f::T, g::T, h::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Quaternion{T}(
            Binion{T}(f, g),
            Binion{T}(h, zero(T))
        )
    )
end

function Exo1Quaternion(a::T, b::T, c::T, d::T, f::T, g::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Quaternion{T}(
            Binion{T}(f, g),
            zero(Binion{T})
        )
    )
end

function Exo1Quaternion(a::T, b::T, c::T, d::T, f::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Quaternion{T}(
            Binion{T}(f, zero(T)),
            zero(Binion{T})
        )
    )
end

function Exo1Quaternion(a::T, b::T, c::T, d::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        zero(Quaternion{T})
    )
end

function Exo1Quaternion(a::T, b::T, c::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, b),
            Binion{T}(c, zero(T))
        ),
        zero(Quaternion{T})
    )
end

function Exo1Quaternion(a::T, b::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, b),
            zero(Binion{T})
        ),
        zero(Quaternion{T})
    )
end

function Exo1Quaternion(a::T) where T <: Real
    Exo1Quaternion{T}(
        Quaternion{T}(
            Binion{T}(a, zero(T)),
            zero(Binion{T})
        ),
        zero(Quaternion{T})
    )
end

function Exo1Quaternion(a::Real, b::Real)
    Exo1Quaternion(promote(a, b)...)
end

function Exo1Quaternion(a::Real, b::Real, c::Real)
    Exo1Quaternion(promote(a, b, c)...)
end

function Exo1Quaternion(a::Real, b::Real, c::Real, d::Real)
    Exo1Quaternion(promote(a, b, c, d)...)
end

function Exo1Quaternion(a::Real, b::Real, c::Real, d::Real, f::Real)
    Exo1Quaternion(promote(a, b, c, d, f)...)
end

function Exo1Quaternion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real)
    Exo1Quaternion(promote(a, b, c, d, f, g)...)
end

function Exo1Quaternion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real)
    Exo1Quaternion(promote(a, b, c, d, f, g, h)...)
end

function Exo1Quaternion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real, j::Real)
    Exo1Quaternion(promote(a, b, c, d, f, g, h, j)...)
end

function Exo1Quaternion(z::Exo1Real{T}) where T <: Real
    Exo1Quaternion(z.l, zero(T), zero(T), zero(T), z.r, zero(T), zero(T), zero(T))
end

"""

    unreal(z::Exo1Quaternion)

The unreal part of an exo-quaternion is a 7-vector.
"""
function unreal(z::Exo1Quaternion)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::Exo1Quaternion)
    print(io, "[1: ")
    print(io, z.l.l.l)
    print(io, ", i: ")
    print(io, z.l.l.r)
    print(io, ", j: ")
    print(io, z.l.r.l)
    print(io, ", ij: ")
    print(io, z.l.r.r)
    print(io, ", a: ")
    print(io, z.r.l.l)
    print(io, ", ia: ")
    print(io, z.r.l.r)
    print(io, ", ja: ")
    print(io, z.r.r.l)
    print(io, ", (ij)a: ")
    print(io, z.r.r.r)
    print(io, "]")
end

function random(::Type{Exo1Quaternion{T}}) where T <: Real
    Exo1Quaternion{T}(random(Quaternion{T}), random(Quaternion{T}))
end


text:
include("consts.jl")
include("model.jl")
include("equations.jl")

idealmodel(model::GERG2008) = model
molecular_weight(model::GERG2008,z=SA[1.0]) = comp_molecular_weight(mw(model),z)
split_model(model::GERG2008) = simple_split_model(model)
mw(model::GERG2008) = model.Mw
function Base.show(io::IO, mime::MIME"text/plain",sp::GERG2008)
    ln = length(sp.components)
    println(io,"GERG008 model with ",ln," component",ifelse(isone(ln),"","s"),":")
    tick = "\""
    for i in 1:ln-1
        println(io,tick,sp.components[i],tick)
    end
    print(io,tick,sp.components[ln],tick)
end

function Base.show(io::IO,model::GERG2008)
    return eosshow(io,model)
end

#uses extended corresponding states and the propane ancilliary eqs
#naively uses Tc as scaling factor
#ideally we would use a eos-agnostic methodology, the noro-frenkel law could give us some insight
#also,, if an eos agnostic shape factor is implemented, we could use GERG2008 as a provider of shape factors

function x0_sat_pure(model::GERG2008,T)
    Ts = T_scale(model)
    vs = _v_scale(model)*0.001 #remember, vc constants in L/mol
    h = vs*5000.0
    T0 = 369.89*T/Ts
    vl = (1.0/_propaneref_rholsat(T0))*h
    vv = (1.0/_propaneref_rhovsat(T0))*h
    return [log10(vl),log10(vv)]
end

function x0_volume_liquid(model::GERG2008,T,z)
    return 1.01*lb_volume(model,z)
end

export GERG2008


text:
# -*- coding: utf-8 -*-
# # Augmenting Model Workflows
#
# We can apply our model augmentation framework to models that are compositions of component models.

# We are going to use the model augmentation presented in `examples/agentgraft.jl` as a baseline simulation and build a workflow to compose that model with the example in `examples/polynomial_regression.jl`. It is strongly recommended that you understand those examples before following this notebook.

# This example combines an agent based model of SIR diseases with a statistical model of polynomial regression to quantify the response of the agent based model with respect to one of its parameters. The input models have to be composed carefully in order to make the software work.
#
# As taught by the scientific computing education group [Software Carpentry](https://swcarpentry.github.io/), the best practice for composing scientific models is to have each component write files to disk and then use a workflow tool such as [Make](https://swcarpentry.github.io/make-novice/) to orchestrate the execution of the modeling scripts.
#
# An alternative approach is to design modeling frameworks for representing the models. The problem with this avenue becomes apparent when models are composed. The frameworks must be interoperable in order to make combined models. ModelTools avoids this problem by representing the models as code and manipulating the codes. The interoperation of two models is defined by user supplied functions in a fully featured programming language.

# Let $m_1,m_2$ be models, and $t_1,t_2$ be tranformations and define $M_i = t_i(m_i)$. If we denote the creation of pipelines with the function composition symbol $g\circ f$ then we want to implement everything such that the following diagram commutes.
#

#
# This example shows how you can use a pipeline to represent the combination of models and then apply combinations of transformations to that pipeline. Transforming models and composing them into pipelines are two operations that commute, you can transform then compose or compose and then transform.

using SemanticModels.Parsers
using SemanticModels.ModelTools
using SemanticModels.ModelTools.ExpStateModels
import Base: push!
import SemanticModels.ModelTools: model, isexpr
using Random

samples = 100
nsteps = 25
finalcounts = Any[]

println("Running Agent Based Simulation Augmentation Demo")
println("================================================")
println("demo parameters:\n\tsamples=$samples\n\tnsteps=$nsteps")

# ## Baseline SIRS model
#
# Here is the baseline model, which is read in from a text file. You could instead of using `parsefile` use a `quote/end` block to code up the baseline model in this script.
#
#

expr = parsefile("../examples/agentbased.jl")
m = model(ExpStateModel, expr)

function returns(block::Vector{Any})
    filter(x->(head(x)==:return), block)
end
returntuples = (bodyblock(filter(x->isa(x, Expr), findfunc(m.expr, :main))[end])
    |> returns
    .|> x-> x.args[1].args )
push!(returntuples[1], :((ρ=ρ, μ=μ, n=n)))

magents = m
println("\nRunning basic model")
AgentModels = eval(m.expr)
@show AgentModels
for i in 1:samples
    println(("======= . Simulation $i  ========"))
    newsam, counts, params = AgentModels.main(nsteps)
    push!(finalcounts, (model=:basic, counts=counts, params=params))
end

@show length(finalcounts)
finalcounts

ModelTools.invoke(magents, 10)[2:end]

# ## Statistical Regression Model
# The following expression defines a univariate polynomial regression model of degree 0, which just computes the average of target variable. This model can be augmented to an polynomial regression model using transformations
# $T_1,T_x$ which will be defined later.

# +

expr = quote
    module Regression
    using Random
    using LsqFit
    using LinearAlgebra

    function f(x, β)
        # This .+ node is added so that we have something to grab onto
        # in the metaprogramming. It is the ∀a .+(a) == a.
        return .+(β[1].* x.^0)
    end

    function sample(g::Function, n)
        x = randn(Float64, n)
        target = g(x) .+ randn(Float64, n[1])./1600
        return x, target
    end

    function describe(fit)
        if !fit.converged
            error("Did not converge")
        end
        return (β = fit.param, r=norm(fit.resid,2), n=length(fit.resid))
    end
    #setup
    a₀ = [1.0]
    function main(X, target)
        #solving
        fit = curve_fit(f, X, target, a₀)#; autodiff=:forwarddiff)
        result = describe(fit)
        return fit, result
    end
end
end
Regression = eval(expr.args[2])
# -

function connector(finalcounts)
    n = length(finalcounts)
    Data = zeros(n,length(finalcounts[1].counts))
    @show size(Data)
    for i in 1:n
        c = finalcounts[i].counts
        Data[i, : ] = map(last, c)
    end
    # multivariate regression not yet supported
    # X = Data[:, 1:end-2]
    # Y = Data[:, end]
    # @assert(size(X) == (n,size(Data,2)-2))
    # @assert(size(Y) == (n,))
    X = Data[:, 1]
    Y = Data[:, 2]
    @assert size(X,1) == size(Y,1)
    return X,Y
end

X,Y = connector(finalcounts)
collect(zip(X,Y))

fit, result = Regression.main(connector(finalcounts)...)
result

if sum(fit.resid)/length(fit.resid) > 1e-8
    @warn "Mean residual is large, regression might have failed"
end

# ## Implementation Details
# The following code is the implementation details for representing the models as an `AbstractProblem` and representing the transformations as `Product{Tuple{Pow{Int}, Pow{Int}}}` and applying the transformations onto the models.
#
# See the `examples/polynomial_regression.jl` example for details of what this code does.

using LinearAlgebra
using SemanticModels
using SemanticModels.ModelTools
using SemanticModels.ModelTools.Transformations
import SemanticModels.ModelTools: model, AbstractModel, isexpr
import SemanticModels.Parsers: findfunc, findassign
import Base: show

# +
"""    Lsq

A program that solves min_β || f(X,β) - y ||_2

Example:

`f(X, β) = β[1].*X.^p .+ β[2].*X.^q`

See also [`(t::Pow)(m::MultivariateLsq)`](@ref)
"""
struct Lsq <: AbstractModel
    expr
    f
    coefficient
    p₀
end

function show(io::IO, m::Lsq)
    write(io, "Lsq(\n  f=$(repr(m.f)),\n  coefficient=$(repr(m.coefficient)),\n  p₀=$(repr(m.p₀))\n)")
end

function model(::Type{Lsq}, expr::Expr)
    if expr.head == :block
        return model(Lsq, expr.args[2])
    end
    objective = :l2norm
    f = callsites(expr, :curve_fit)[end].args[2]
    coeff = callsites(expr, f)[1].args[end]
    p₀ = callsites(expr, :curve_fit)[end].args[end]
    return Lsq(expr, f, coeff, p₀)
end

"""    poly(m::Lsq)::Expr

find the part of the model that implements the polynomial model for regression.
"""
function poly(m::Lsq)
    func = findfunc(m.expr, m.f)[1]
    poly = func.args[2].args[end].args[1]
    return poly
end

"""    (t::Pow)(m::Lsq)

Example:

If `m` is a program implementing `f(X, β) = β[1]*X^p + β[2]*X^q`

a) and `t = Pow(2)` then `t(m)` is the model implementing
`f(X, β) = β[1]*X^p+2 + β[2]*X^q+q`.

"""
function (t::Pow)(m::Lsq)
    p = poly(m)
    for i in 2:length(p.args)
        slot = p.args[i]
        pow = callsites(slot, :(.^))
        pow[end].args[3] += t.inc
    end
    return m
end

# +
struct AddConst <: Transformations.Transformation end

"""    (c::AddConst)(m::MultivariateLsq)

Example:

If `m` is a program implementing `f(X, β) = β[1]*X^p + β[2]*X^q`

a) and `c = AddConst()` then `c(m)` is the model implementing
`f(X, β) = β[1]*X^p + β[2]*X^0`.

"""
function (c::AddConst)(m::Lsq)
    p = poly(m)
    ix = map(t->t.args[2].args[2], p.args[2:end])
    i = maximum(ix)+1
    @show p
    push!(p.args, :(β[$i].*x.^0))
    assigns = findassign(m.expr, m.p₀)
    @show assigns
    b = assigns[end].args[2].args
    push!(b, 1)
    return m
end
# -

# ## Model Augmentation yields polynomial regression
# Given our transformations `f(x) -> xf(x)` and `f(x) -> f(x) + beta` we are able to generate all possible polynomial regression using composition of these transformations.

# Let's build an instance of the model object from the code snippet expr
m = model(Lsq, deepcopy(expr))
mstats = deepcopy(m)
@show m
poly(m)

# Some *generator elements* will come in handy for building elements of the transformation group.
# $T_x,T_1$ are *generators* for our group of transformations $T = \langle T_x, T_1 \rangle$. $T_1$ adds a constant to our polynomial and $T_x$ increments all the powers of the terms by 1. Any polynomial can be generated by these two operations. The proof of Horner's rule for evaluating $p(x)$ gives a construction for how to create $f(x,\beta) = p(x)$ from these two operations.

@show Tₓ = Pow(1)
@show T₁ = AddConst()
Tₓ, T₁

m′ = deepcopy(m)
Tₓ(m′)
T₁(m′)
@show poly(m)
Regression = eval(m′.expr)

@show m′.expr
fit′, result′ = Regression.main(connector(finalcounts)...)
result′

result′.r

# # Model workflows
#
# Models can be chained together into workflows, the most basic type is a pipeline where the outputs from model $m_i$ are passed to model $m_{i+1}$. One area where traditional modeling frameworks get in trouble is the fact that the connections between the models can be arbitrarily complex. Thus any modeling framework that supports worflows, must embed a programming language for describing the connectors between the steps of the workflow.
#
# Since we are already embedded in Julia, we will use regular Julia functions for the connectors.
#
# Mathematically, a pipeline is defined as $r_n = P(m_1,\dots,m_n, c_1,\dots,c_n)$ based on the recurrence,
#
# $r_0 = m_1(c)$ where $c$ is a constant value, and
#
# $r_i = m_i(c_i(r_{i-1}))$
#
#
# We store the values of $r_i$ in the field results so that they can be accessed later by visualization and  analysis programs.

# +
module Pipelines
using SemanticModels.ModelTools
using Random
struct Pipeline <: AbstractModel
    steps
    connectors
    results
end

function run!(p::Pipeline)
    stages = length(p.steps)
    connectors = p.connectors
    for s in 1:stages
        data = p.results[end]
        r = connectors[s](p.steps[s], data...)
        push!(p.results, r)
    end
end

function reset!(p::Pipeline)
    while length(p.results) > 1
        pop!(p.results)
    end
    return p
end
end
# -

# ## Running the baseline workflow
#
# This workflow connects the two models so that we simulate the agent based model and then perform a regression on the outputs.

P = Pipelines.Pipeline(deepcopy.([magents, mstats]),
    [(m, args...) -> begin
            Random.seed!(42)
            results = Any[]
            Mod = eval(m.expr)
            for i in 1:samples
                r = Base.invokelatest(Mod.main, args...)
                push!(results, (model=:basic, counts=r[2], params=r[3]))
                #push!(results, r)
            end
            return [results]
                end,
        (m, results...) -> begin
            data = connector(results...)
            Mod = eval(m.expr)
            Base.invokelatest(Mod.main, data...) end
        ],
        Any[(10)]
        )

# Warning: Pipelines can only be run once. Recreate the pipeline and run it again if necessary.

Pipelines.run!(P)

for r in P.results[end]
    println(r)
    println("==========")
end

Tₓ(P.steps[2])
T₁(P.steps[2])
Tₓ(P.steps[2])
T₁(P.steps[2])

Pipelines.reset!(P)
Pipelines.run!(P)
P.results[end][2]

# 1. TODO implement the composition monoid on ExpStateModels
# 2. TODO use Product acting on Pipeline to represent Tuple.coordinates == Pipeline.steps

# Product = ModelTools.Transformations.Product
function (t::Product)(m::Pipelines.Pipeline)
    for (i, s) in enumerate(m.steps)
        t.dims[i](s)
    end
end

# ## Adding the Dead State
#
#
#
# We are going to add an additional state to the model to represent the infectious disease fatalities. The user must specify what that concept means in terms of the name for the new state and the behavior of that state. `D` is a terminal state for a finite automata.

function addstate!(m::ExpStateModel)
    println("\nThe system states are $(m.states.args)")
    println("\nAdding un estado de los muertos")

    put!(m, ExpStateTransition(:D, :((x...)->:D)))

    println("\nThe system states are $(m.states.args)")
    # once you are dead, you are dead forever
    println("\nThere is no resurrection in this model")
    println("\nInfected individuals recover or die in one step")

    # replace!(m, ExpStateTransition(:I, :((x...)->rand(Bool) ? :D : :I)))
    m[:I] = :((x...)->begin
            roll = rand()
            if roll < ρ
                return :R
            elseif rand(Bool)
                return :D
            else
                return :I
            end
        end
    )
    @show m[:I]
    return m
end

# Some utilities for manipulating functions at a higher level than expressions.

# +

struct Func end

function push!(::Func, func::Expr, ex::Expr)
    push!(bodyblock(func), ex)
end
# -

# ## Population Growth
#
# Another change we can make to our model is the introduction of population growth. Our model for population is that on each timestep, one new suceptible person will be added to the list of agents. We use the `tick!` function as an anchor point for this transformation.
#
#

function addgrowth!(m::ExpStateModel)
    println("\nAdding population growth to this model")
    stepr = filter(x->isa(x,Expr), findfunc(m.expr, :tick!))[1]
    @show stepr
    push!(Func(), stepr, :(push!(sm.agents, :S)))
    println("------------------------")
    @show stepr;
    return m
end

# ## Transforming the Pipeline
#
# We can apply the tranformations that we defined onto the pipeline. Remember that $T_1\times T_2$ acts on a pipeline by creating $P(T_1(m_1),T_2(m_2), c_1,c_2)$. The current implementation does not support transforming the connectors, but that would be straightforward to add.

function connector(finalcounts, i, j)
    n = length(finalcounts)
    Data = zeros(n,length(finalcounts[1].counts))
    params = zeros(n,length(finalcounts[1].params))
    @show size(Data)
    for i in 1:n
        c = finalcounts[i].counts
        Data[i, : ] = map(last, c)
        params[i,:] = collect(map(float, finalcounts[i].params))
    end
    # multivariate regression not yet supported
    # X = Data[:, 1:end-2]
    # Y = Data[:, end]
    # @assert(size(X) == (n,size(Data,2)-2))
    # @assert(size(Y) == (n,))
    X = params[:, i]
    Y = Data[:, j]
    # normalization
    #Y = Y ./ [sum(Data[i, :]) for i in 1:n]
    @assert size(X,1) == size(Y,1)
    return X,Y
end
P = Pipelines.Pipeline(deepcopy.([magents, mstats]),
    [(m, args...) -> begin
            Random.seed!(42)
            results = Any[]
            Mod = eval(m.expr)
            @show Mod
            for i in 1:samples
                r = Base.invokelatest(Mod.main, args...)
                push!(results, (model=:basic, counts=r[2],params=r[3]))
                #push!(results, r)
            end
            return [results]
                end,
        (m, results...) -> begin
            data = connector(results..., 1, 4)
            Mod = eval(m.expr)
            Base.invokelatest(Mod.main, data...)
        end
        ],
        Any[(10)]
        )
println("\nInitial Pipeline")
println("----------------")
@show P.steps[1].states
@show poly(P.steps[2])
println("\n\nApplying the first pair of transformations")
println(     "------------------------------------------")
Product((addstate!, one(Pow)))(P)
@show P.steps[1].states
@show poly(P.steps[2])
println("\n\nApplying the second pair of transformations")
println(    "-------------------------------------------")
Product((x->x, AddConst()))(P)
Product((addgrowth!, one(Pow)))(P)
Product((x->x, AddConst()))(P)
Product((x->x, one(Pow)))(P)
Product((x->x, AddConst()))(P)
println("\n\nThe final model state")
println(     "---------------------")
println(filter(isexpr, findfunc(P.steps[1].expr, :tick!))[end])
@show poly(P.steps[2]);

# ## Running the new pipeline
#
# now that we have transformed the pipeline we run it to build the results

Pipelines.run!(P)
P.results[end][2]

# Here is the data we observed when running the first stage of the pipeline, stage two fits a polynomial to these observations

table = map(x->(round(x.params.ρ, digits=4), last(x.counts[end])), P.results[2][1]) |> sort
try
    using Plots
catch
    @warn "Plotting is not available, make a table"
    for t in table
        println(join(t, "\t"))
    end
end

# ## Results
#
# The regression model that we have trained based on the simulated data from the SIRD model with population growth can be presented as a polynomial sampled over the domain. We construct this table to show the nonlinear dependence of the model on the recovery parameter $\rho$. The best fitting polynomial is shown below.

using Printf
eval(:(f(x,β) = $(poly(P.steps[2]))))
xdomain = (0.0:0.05:1.0)
println("ρ\tf(ρ,β)\n==============")
xŷ = zip(xdomain, f(xdomain, P.results[end][2].β))
z = collect(map(x->(@sprintf("%0.2f", x[1]),
                    @sprintf("%7.3f", x[2])),
        xŷ))
for t in z
    println(join(t, "\t"))
end

@info "Loading Plots, this may take a while"
using Plots

@info "Making plots, this may take a while"
p = scatter(first.(table), last.(table), label="obs")
plot!(first.(xŷ), last.(xŷ), label="fit")
xlabel!(p, "Probability of Recovery")
ylabel!(p, "Deaths")
println("β: ", P.results[end][2].β, "\n", string(poly(P.steps[2])))
p

# ## Conclusions
#
# This example shows that the SemanticModels approach to post-hoc modeling frameworks can enable metamodeling which is the combinations of model composed using different technologies into a coherent modeling workflow. Our ModelTools provides the basic building blocks for representing models and transformations in such a way that they transformations can be composed and models can be combined. Composition of transformations respects the combination of models. In this case the Product of transformation respects the Pipeline of models. Such that you can transform the models and then pipeline them, or pipeline them and then transform them.
#
# This example combined an agent based model of SIR diseases with a statistical model of polynomial regression to quantify the response of the agent based model with respect to one of its parameters. The input models have to be composed carefully in order to make the software work.
#
# As taught by the scientific computing education group [Software Carpentry](https://swcarpentry.github.io/), the best practice for composing scientific models is to have each component write files to disk and then use a workflow tool such as [Make](https://swcarpentry.github.io/make-novice/) to orchestrate the execution of the modeling scripts.
#
# An alternative approach is to design modeling frameworks for representing the models. The problem with this avenue becomes apparent when models are composed. The frameworks must be interoperable in order to make combined models. ModelTools avoids this problem by representing the models as code and manipulating the codes. The interoperation of two models is defined by user supplied functions in a fully featured programming language.
#
# SemanticModels.jl also provides transformations on these models that are grounded in category theory and abstract algebra. The concepts of category theory such as Functors and Product Categories allow us to build a general framework fit for any modeling task. In the language of category theory, the Pipelining functor on models commutes with the Product functor on transformations.
#
# This examples shows that metamodeling is feasible with SemanticModels and that the algebras of model transformations can be preserved when acting on metamodel workflows.


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/constraints.jl
# Defines constraints supported Optimizer and how to store them.
#############################################################################

##### Supports function and add_constraint for scalar affine functions
MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{S}) where {S <: INEQ_SETS} = true

macro define_addconstraint_linear(function_type, set_type, array_name, count_name)
    quote
        function MOI.add_constraint(m::Optimizer, func::$function_type, set::$set_type)
            check_inbounds!(m, func)
            push!(m._input_problem.$(array_name), (func, set))
            m._input_problem._last_constraint_index += 1
            m._input_problem.$(count_name) += 1
            indx = CI{$function_type, $set_type}(m._input_problem._last_constraint_index)
            return indx
        end
    end
end

@define_addconstraint_linear SAF LT _linear_leq_constraints _linear_leq_count
@define_addconstraint_linear SAF GT _linear_geq_constraints _linear_geq_count
@define_addconstraint_linear SAF ET _linear_eq_constraints _linear_eq_count

##### Supports function and add_constraint for scalar quadratic functions
MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{S}) where {S <: INEQ_SETS} = true

macro define_addconstraint_quadratic(function_type, set_type, array_name, count_name)
    quote
        function MOI.add_constraint(m::Optimizer, func::$function_type, set::$set_type)
            check_inbounds!(m, func)
            push!(m._input_problem.$(array_name), (func, set))
            m._input_problem._last_constraint_index += 1
            m._input_problem.$(count_name) += 1
            indx = CI{$function_type, $set_type}(m._input_problem._last_constraint_index)
            return indx
        end
    end
end

@define_addconstraint_quadratic SQF LT _quadratic_leq_constraints _quadratic_leq_count
@define_addconstraint_quadratic SQF GT _quadratic_geq_constraints _quadratic_geq_count
@define_addconstraint_quadratic SQF ET _quadratic_eq_constraints _quadratic_eq_count

##### Supports function and add_constraint for conic functions
#=
const CONE_SETS = Union{SOC}
MOI.supports_constraint(::Optimizer, ::Type{VECOFVAR}, ::Type{S}) where {S <: CONE_SETS} = true

function MOI.add_constraint(m::Optimizer, func::VECOFVAR, set::SOC)

    if length(func.variables) !== set.dimension
        error("Dimension of $(s) does not match number of terms in $(f)")
    end

    check_inbounds!(m, func)
    push!(m._input_problem._conic_second_order, (func, set))
    m._input_problem._last_constraint_index += 1
    m._input_problem._conic_second_order_count += 1

    return CI{VECOFVAR, SOC}(m._input_problem._last_constraint_index)
end
=#


text:
using MAT
using StatsBase
std = StatsBase.std
trials = 10
include("../include/HyperLocal/HyperLocal.jl")
include("../include/HyperLocal/SparseCardHyperLocal.jl")
include("../src/hypergraph-cond-utils.jl")
## Read in the matrix, if it isn't read in already
include("load_stackoverflow.jl")

## Mean runtime plots
for type = ["sqrt", "x9", "clique"]
typename = type
if type == "x9"
    typename = "x^{0.9}"
end
run_mean = zeros(45,3)
run_std = zeros(45,3)
f1_mean = zeros(45,3)
f1_std = zeros(45,3)
for jj = 1:3
    Druns = zeros(45,trials)
    Df1s = zeros(45,trials)
    for trial = 1:trials
        M = matread("Output/SCHL_$(type)_45_$trial.mat")
        stats = M["sparse$(jj)_stats"]
        drun = stats[:,6]
        df1 = stats[:,3]
        Druns[:,trial] = drun
        Df1s[:,trial] = df1
    end
    index = jj
    run_mean[:,index] = mean(Druns,dims = 2)
    run_std[:,index] = std(Druns,dims = 2)
    f1_mean[:,index] = mean(Df1s,dims = 2)
    f1_std[:,index] = std(Df1s,dims = 2)
end

function plotunif(pl, y, label, color,ms)
    scatter!(pl,y,markersize = ms, markerstrokecolor = color, color = color, markershape = :circle, label = label)
end
stepy = 2
s1 = 500
s2 = 350
p = sortperm(run_mean[:,1])
ln = names[p]
pl = plot(grid = false)
plotunif(pl,run_mean[p,3],"ε = 0.01",:black,5)
plotunif(pl,run_mean[p,2],"ε = 0.1",:grey,5)
plotunif(pl,run_mean[p,1],"ε = 1.0",:green,5)
plot!(pl,yscale = :log10,legend = :topleft,title = type*" function, mean runtime",
    ylabel = "Runtime (seconds)",size = (s1,s2))
mean(abs.(f1_mean[:,1]-f1_mean[:,2]))
plot!(pl,xticks = (1:stepy:length(ln), ln[1:stepy:length(ln)]),xrotation = 40,
xtickfont=font(8),
ytickfont=font(11),
guidefont=font(12),
titlefont=font(12),
legendfont=font(9))

savefig("Figures/$(type)_runtime_plot.pdf")
jump = mean(run_mean[:,2]./run_mean[:,1])
jump2 = mean(run_mean[:,3]./run_mean[:,2])
println("")
println("Runtime jumps = $jump, $jump2")
end

## Meanwhile, F1 doesn't change
for type = ["sqrt", "x9", "clique"]
run_mean = zeros(45,3)
run_std = zeros(45,3)
f1_mean = zeros(45,3)
f1_std = zeros(45,3)
for jj = 1:3
    Druns = zeros(45,trials)
    Df1s = zeros(45,trials)
    for trial = 1:trials
        M = matread("Output/SCHL_$(type)_45_$trial.mat")
        stats = M["sparse$(jj)_stats"]
        drun = stats[:,6]
        df1 = stats[:,3]
        Druns[:,trial] = drun
        Df1s[:,trial] = df1
    end
    index = jj
    run_mean[:,index] = mean(Druns,dims = 2)
    run_std[:,index] = std(Druns,dims = 2)
    f1_mean[:,index] = mean(Df1s,dims = 2)
    f1_std[:,index] = std(Df1s,dims = 2)
end

function plotunif(pl, y, label, color,ms)
    scatter!(pl,y,markersize = ms, markerstrokecolor = color, color = color, markershape = :circle, label = label)
end
stepy = 2
s1 = 500
s2 = 350
p = sortperm(run_mean[:,1])
ln = names[p]
pl = plot(grid = false)
plotunif(pl,f1_mean[p,3],"ε = 0.01",:black,5)
plotunif(pl,f1_mean[p,2],"ε = 0.1",:grey,5)
plotunif(pl,f1_mean[p,1],"ε = 1.0",:green,5)
typename = type
if type == "x9"
    typename = "x^{0.9}"
end
plot!(pl,yscale = :identity,legend = :topleft,title = typename*" function, mean F1",
    ylabel = "F1 Score",size = (s1,s2))
mean(abs.(f1_mean[:,1]-f1_mean[:,2]))
plot!(pl,xticks = (1:stepy:length(ln), ln[1:stepy:length(ln)]),xrotation = 40,
xtickfont=font(8),
ytickfont=font(11),
guidefont=font(12),
titlefont=font(12),
legendfont=font(9))

savefig("Figures/$(type)_f1_plot.pdf")
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
@testset "Running PELT tests..." begin

    n = 1000;        # Number of samples
    λ = 100;         # Frequencey of changes

    ########################
    # Normal mean segments #
    ########################
    μ, σ = Normal(0.0, 10.0), 1.0
    sample, cps = @changepoint_sampler n λ Normal(μ, σ)
    seg_costs = NormalMeanSegment(sample);
    @test_nowarn PELT(seg_costs, n);
    @test_nowarn @PELT sample Normal(:?, σ)

    #######################
    # Normal var segments #
    #######################
    μ, σ = 1.0, Uniform(2.0, 15.0)
    sample, cps = @changepoint_sampler n λ Normal(μ, σ)
    seg_costs = NormalVarSegment(sample, μ)
    @test_nowarn PELT(seg_costs, n)
    @test_nowarn @PELT sample Normal(μ, :?)

    ############################
    # Exponential changepoints #
    ############################
    μ = Uniform(0.0, 10.0)
    sample, cps = @changepoint_sampler n λ Exponential(μ)
    seg_costs = ExponentialSegment(sample)
    @test_nowarn PELT(seg_costs, n)
    @test_nowarn @PELT sample Exponential(:?)

    ########################
    # Poisson changepoints #
    ########################
    μ = Uniform(0.0, 10.0)
    sample, cps = @changepoint_sampler n λ Poisson(μ)
    seg_costs = PoissonSegment(sample)
    @test_nowarn PELT(seg_costs, n)
    @test_nowarn @PELT sample Poisson(:?)

end