text:
add_brackets(arr::AbstractArray, vars) = [add_brackets(element, vars) for element in arr]

add_brackets(s::Any, vars) = return s

function add_brackets(ex::Expr, vars)
    ex = postwalk(x -> x in vars ? "\\left[ $(convertSubscript(x)) \\right]" : x, ex)
    return ex
end

"""
    render(::LaTeXString; debug=false, name=tempname(), command="\\Large")

Display a standalone PDF with the given input.
"""
function render(s::LaTeXString; debug=false, name=tempname(), command="\\Large")
    doc = """
    \\documentclass[varwidth=100cm]{standalone}
    \\usepackage{amssymb}
    \\usepackage{amsmath}
    $(occursin("\\ce{", s) ? "\\usepackage{mhchem}" : "")
    \\begin{document}
    {
        $command
        $s
    }
    \\end{document}
    """
    doc = replace(doc, "\\begin{align}"=>"\\[\n\\begin{aligned}")
    doc = replace(doc, "\\end{align}"=>"\\end{aligned}\n\\]")
    doc = replace(doc, "\\require{mhchem}\n"=>"")
    open("$(name).tex", "w") do f
        write(f, doc)
    end
    cd(dirname(name)) do
        cmd = `lualatex --interaction=batchmode $(name).tex`
        debug || (cmd = pipeline(cmd, devnull))
        run(cmd)
    end
    if Sys.iswindows()
        run(`cmd /c "start $(name).pdf"`, wait=false)
    elseif Sys.islinux()
        run(`xdg-open $(name).pdf`, wait=false)
    elseif Sys.isapple()
        run(`open $(name).pdf`, wait=false)
    elseif Sys.isbsd()
        run(`xdg-open $(name).pdf`, wait=false)
    end
    return nothing
end

function expr_to_array(ex)
    ex.head == :typed_vcat && (ex = Expr(:vcat, ex.args[2:end]...))
    ex.head == :typed_hcat && (ex = Expr(:hcat, ex.args[2:end]...))
    ex.head == :ref && (ex = Expr(:vect, ex.args[2:end]...))
    ## If it is a matrix
    if ex.args[1] isa Expr && ex.args[1].head == :row
        return eval(ex.head)(map(y -> permutedims(y.args), ex.args)...)
    else
        if ex.head == :hcat
            return hcat( ex.args...)
        elseif ex.head in [:vcat, :vect]
            return vcat( ex.args...)
        end
    end
end


text:
"""
    function vertex3(para, extK = [DiagTree.getK(para.totalLoopNum, 1), DiagTree.getK(para.totalLoopNum, 2)],
        subdiagram = false; name = :Γ3, chan = [PHr, PHEr, PPr, Alli], resetuid = false)

    Generate 3-vertex diagrams using Parquet Algorithm.
    With imaginary-time variables, all vertex3 generated has the same bosonic Tidx ``extT[1]=para.firstTauIdx`` and the incoming fermionic Tidx ``extT[2]=para.firstTauIdx+1``.

#Arguments
- `para`            : parameters. It should provide internalLoopNum, interactionTauNum, firstTauIdx
- `extK`            : basis of external loops as a vector [bosonic leg, fermionic in, fermionic out].
- `subdiagram`      : a sub-vertex or not
- `name`            : name of the vertex
- `chan`            : vector of channels of the current 4-vertex.
- `resetuid`        : restart uid count from 1

# Output
- A DataFrame with fields :response, :extT, :diagram, :hash.
"""
function vertex3(para, extK=[DiagTree.getK(para.totalLoopNum, 1), DiagTree.getK(para.totalLoopNum, 2)],
    subdiagram=false; name=:Γ3, chan=[PHr, PHEr, PPr, Alli], resetuid=false)

    resetuid && uidreset()
    @assert para.diagType == Ver3Diag
    @assert para.innerLoopNum >= 1 "Only generates vertex corrections with more than one internal loops."
    for k in extK
        @assert length(k) >= para.totalLoopNum "expect dim of extK>=$(para.totalLoopNum), got $(length(k))"
    end
    extK = [k[1:para.totalLoopNum] for k in extK]

    q, Kin = extK[1], extK[2]
    Kout = length(extK) == 3 ? extK[3] : Kin .- q
    @assert ((q ≈ Kin) == false) && ((q ≈ Kout) == false) "The bosonic q cann't be same as the fermionic k. Ohterwise the proper diagram check will fail!"
    extK = [q, Kin, Kout]

    if Proper in para.filter
        para = reconstruct(para, transferLoop=q)
    end

    t0 = para.firstTauIdx
    vertex3 = DataFrame()

    # if para.innerLoopNum == 0
    #     push!(vertex3, (response = UpUp, extT = (t0, t0, t0), diagram = ver3diag))
    # end

    if (para.extra isa ParquetBlocks) == false
        para = reconstruct(para, extra=ParquetBlocks())
    end

    K = zero(q)
    LoopIdx = para.firstLoopIdx
    K[LoopIdx] = 1.0
    # extT = (t0, t0 + 1)
    legK = [Kin, Kout, K, K .+ q]

    ######################## Π0 = GG #########################################
    for (oVer4, oGin, oGout) in orderedPartition(para.innerLoopNum - 1, 3, 0)
        # ! Vertex4 must be in the first place, because we want to make sure that the TinL of the vertex4 start with t0+1

        idx, maxLoop = findFirstLoopIdx([oVer4, oGin, oGout], LoopIdx + 1)
        @assert maxLoop <= para.totalLoopNum "maxLoop = $maxLoop > $(para.totalLoopNum)"
        Ver4Kidx, GinKidx, GoutKidx = idx

        ver4t0 = para.hasTau ? para.firstTauIdx + 1 : para.firstTauIdx
        idx, maxTau = findFirstTauIdx([oVer4, oGin, oGout], [Ver4Diag, GreenDiag, GreenDiag], ver4t0, para.interactionTauNum)
        @assert maxTau <= para.totalTauNum "maxTau = $maxTau > $(para.totalTauNum)"
        Ver4Tidx, GinTidx, GoutTidx = idx

        if isValidG(para.filter, oGin) && isValidG(para.filter, oGout)
            paraGin = reconstruct(para, diagType=GreenDiag, innerLoopNum=oGin,
                firstLoopIdx=GinKidx, firstTauIdx=GinTidx)
            paraGout = reconstruct(para, diagType=GreenDiag, innerLoopNum=oGout,
                firstLoopIdx=GoutKidx, firstTauIdx=GoutTidx)
            paraVer4 = reconstruct(para, diagType=Ver4Diag, innerLoopNum=oVer4,
                firstLoopIdx=Ver4Kidx, firstTauIdx=Ver4Tidx)
            ver4 = vertex4(paraVer4, legK, chan, true)
            if isnothing(ver4) || isempty(ver4)
                continue
            end

            if para.hasTau
                @assert all(x -> x[INL] == ver4t0, ver4[:, :extT]) "The TinL of the inner Γ4 must be firstTauIdx+1"
            end

            #transform extT coloum into extT for Vertex4 and the extT for Gin and Gout
            df = transform(ver4, :extT => ByRow(x -> [(t0, x[INL], x[OUTL]), (t0, x[INR]), (x[OUTR], t0)]) => [:extT, :GinT, :GoutT])

            groups = mergeby(df, [:response, :GinT, :GoutT, :extT], operator=Sum())

            for v4 in eachrow(groups)
                response = v4[:response]
                @assert response == UpUp || response == UpDown
                #type: Instant or Dynamic
                ver3id = Ver3Id(para, response, k=extK, t=v4[:extT])
                gin = green(paraGin, K, v4[:GinT], true, name=:Gin)
                gout = green(paraGout, K .+ q, v4[:GoutT], true, name=:Gout)
                @assert gin isa Diagram && gout isa Diagram

                ver3diag = Diagram(ver3id, Prod(), [gin, gout, v4[:diagram]], name=name)
                push!(vertex3, (response=response, extT=v4[:extT], diagram=ver3diag))
            end
        end
    end

    if isempty(vertex3)
        return DataFrame(response=[], extT=[], diagram=[])
    end

    Factor = 1 / (2π)^para.loopDim

    ver3 = mergeby(vertex3, [:response, :extT]; name=name, factor=Factor,
        getid=g -> Ver3Id(para, g[1, :response], k=extK, t=g[1, :extT])
    )
    return ver3
end

ext: jl


text:
using XGPaint
using Healpix
using HDF5

## Load halos from HDF5 files, establish a CIB model and cosmology
@time halo_pos, halo_mass = read_halo_catalog_hdf5(
    "/global/cfs/cdirs/sobs/www/users/Radio_WebSky/websky_halos-light.hdf5");
cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0)
model = CIB_Planck2013{Float32}(nside=4096)

## Write one chunk to disk

function write_chunk(
    output_dir, chunk_index, model, cosmo,
    pos, mass, freqs)
    # Allocate some arrays and fill them up for centrals and satellites
    @time sources = generate_sources(model, cosmo, pos, mass);

    # Deposit the sources into maps
    fluxes_cen = Array{Float32, 1}(undef, sources.N_cen)
    fluxes_sat = Array{Float32, 1}(undef, sources.N_sat)
    m = HealpixMap{Float64,RingOrder}(model.nside)

    h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index).h5"), "w") do file
        write(file, "redshift", sources.redshift_cen)
        write(file, "theta", sources.theta_cen)
        write(file, "phi", sources.phi_cen)
    end
    h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index).h5"), "w") do file
        write(file, "redshift", sources.redshift_sat)
        write(file, "theta", sources.theta_sat)
        write(file, "phi", sources.phi_sat)
    end

    # loop over all frequencies and paint sources to appropriate freq map
    @time begin
        for freq in freqs

            XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources,
                fluxes_cen, fluxes_sat)

            # save sources with mass, redshift, angles
            h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file
                write(file, "flux", fluxes_cen)
            end
            h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file
                write(file, "flux", fluxes_sat)
            end

            filename = joinpath(output_dir, "cib_$(freq).fits")

            if chunk_index > 1
                m0 = Healpix.readHealpixMapFromFITS(filename, 1, Float32)
                m.pixels = m.pixels + m0.pixels
            end
            Healpix.saveToFITS(m, "!$(filename)", typechar="D")
        end
    end

end

##
function run_all_chunks(output_dir, halo_pos, halo_mass, freqs; N_chunks=2)
    # provide views into halo positions and masses for chunks of the halos
    N_halos = size(halo_mass, 1)
    chunksize = trunc(Integer, N_halos / N_chunks + 1)
    chunks = chunk(N_halos, chunksize)
    for chunk_index in 1:length(chunks)
    left_ind, right_ind = chunks[chunk_index]
        println("Chunk ", chunk_index, "/", length(chunks),
            " ", left_ind, " ", right_ind)
        pos = @view halo_pos[:, left_ind:right_ind]
        mass = @view halo_mass[left_ind:right_ind]
        write_chunk(output_dir, chunk_index, model, cosmo,
            pos, mass, freqs)
    end
end
## compute on all chunks, on all halos

freqs = [
    "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4",
    "63.9", "67.8", "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143",
    "153", "164", "189", "210", "217", "232", "256", "275", "294", "306", "314",
    "340", "353", "375", "409", "467", "525", "545", "584", "643", "729", "817",
    "857", "906", "994", "1080"
]
# freqs = ["143"]

scratch_dir = joinpath(ENV["SCRATCH"], "cib_sources")
println("SCRATCH: ", scratch_dir)
mkpath(scratch_dir)
mkpath(joinpath(scratch_dir, "sources"))
run_all_chunks(scratch_dir, halo_pos, halo_mass, freqs)

##

max_stars_repo_name: xzackli/XGPaint

max_issues_repo_name: xzackli/XGPaint

max_forks_repo_name: xzackli/XGPaint


text:
#
# Functions to convert between body names and NAIF ID codes,
# and vice versa. This functionality is enabled by the
# CSPICE library, CSPICE_jll, and SPICE.jl!
#
# Mass parameter and radii values are also provided through
# functions. This functionality is enabled by JPL HORIZONS
# and the IAU!
#

export naifcode, naifname

"""
Return the NAIF ID for the provided celestial body name. If no ID for
the provided body is found, a `KeyError` exception is thrown.

!!! note
    This is a simple wrapper around `SPICE.bodn2c`, which is itself
    a wrapper around the CSPICE library's `bod2nc` function! This
    integer code lookup is robust to different capitalizations.
    For example, `"Earth"`, `"earth"`, and `"eArTH"` should all return
    integer code `399`.
"""
function naifcode(name::AbstractString)
    code = SPICE.bodn2c(name)

    if code isa Integer
        return code
    else
        throw(KeyError("No NAIF ID for $name could be found!"))
    end
end

"""
Return the celestial body name for the provided NANIF ID. If no name for
the provided ID is found, a `KeyError` exception is thrown.

!!! note
    This is a simple wrapper around `SPICE.bodc2n`, which is itself
    a wrapper around the CSPICE library's `bodc2n` function!
"""
function naifname(code::Int)
    name = SPICE.bodc2n(code)

    if name isa AbstractString
        return name
    else
        throw(KeyError("No body name for $code could be found!"))
    end
end


text:
struct Senseval3{S}
    filepaths::Vector{S}
end

function Senseval3(dirpath, ver=3.0)
    @assert(isdir(dirpath), dirpath)

    if "wordnet$ver" ∈ readdir(dirpath)
        paths = joinpath.(dirpath,"wordnet$ver", readdir(joinpath(dirpath, "wordnet$ver")))
    else
        paths = joinpath.(dirpath, readdir(dirpath))
    end
    Senseval3(paths)
end

Senseval3() = Senseval3(datadep"Senseval 3")

MultiResolutionIterators.levelname_map(::Type{Senseval3}) = [
    :doc=>1, :contextfile=>1, :context=>1,
    :sent=>2, :sentence=>2,
    :word=>3, :token=>3,
    :char=>4, :character=>5
    ]

function parse_sense_annotated_word_senseval3(line::AbstractString)
    captures = match(r"<wf cmd=done.* pos=(.*) lemma=(.*) wnsn=(.*) lexsn=(\d.*:\d*|[U]|NOANSWER|NOWN).*>(.*).*</wf>", line).captures
    pos, lemma, wnsn, lexsn, word = captures
    if ';' in wnsn
        # Discard Extra Senses
        wnsn = split(wnsn, ';') |> first
        lexsn = split(lexsn, ';') |> first
    end

    SenseAnnotatedWord(pos, lemma, parse(Int, wnsn), lexsn, word)
end

function parse_senseval3file(filename)
    local sent
    lines = @NestedVector(TaggedWord,2)()
    context = Document(intern(basename(filename)), lines)

    # structure

    function new_sentence(line)
        sent = @NestedVector(TaggedWord,1)()
        push!(lines, sent)
    end

    # words
    get_tagged(line) = push!(sent, parse_tagged_word(line))
    get_sense_annotated(line) = push!(sent, parse_sense_annotated_word_senseval3(line))
    get_punc(line) = push!(sent, parse_punc(line))

    # Parse
    subparsers = [
        "<wf cmd=tag"=>          get_tagged,
        "<wf cmd=ignore"=>       get_tagged,
        "<wf cmd=done"=>         line -> contains(line,"lemma=") ? get_sense_annotated(line) : get_tagged(line),
        "<punc>"=>               get_punc,
        "<context"=>             ignore,
        "</context" =>           ignore,
        "<s" =>                  new_sentence,
        "</s" =>                 ignore
    ]

    apply_subparsers(filename,subparsers)

    return context
end

function load(corpus::Senseval3, doc_buffersize=16)
    Channel(;ctype=Document{@NestedVector(TaggedWord, 2), String}, csize=doc_buffersize) do ch
        for fn in corpus.filepaths
            doc = parse_senseval3file(fn)
            put!(ch, doc)
        end
    end
end


text:
export eicxx, efunction

eexport(s) = Expr(:export, unique(s)...)

function efunction(name, args, body)
    Expr(:function, Expr(:call, name, args...), body)
end
eequals(x,y) = Expr(Symbol("="), x, y)
function ecxxnew(constructor::Symbol, args)
    Expr(:macrocall, Symbol("@cxxnew"), Expr(:call, constructor, args...))
end

eicxx(s::String) = Expr(:macrocall, Symbol("@icxx_str"), s)
etyped(obj, T) = Expr(Symbol("::"), obj, T)
function ecxxtype(name)
    # expression for something like
    # Cxx.CppPtr{Cxx.CppValue{Cxx.CxxQualType{Cxx.CppBaseType{name},(false,false,false)},N},(false,false,false)}
    Expr(:curly, :(Cxx.CppPtr),
        Expr(:macrocall, Symbol("@vcpp_str"), name),
        (false, false, false)
    )
end

function emethodcall(obj, method, args)
    callargs = join(("\$($arg)" for arg in args), ", ")
    s = "\$($obj) -> $method($callargs);"
    eicxx(s)
end


text:
module TerminatorsTests

include("utils.jl")
include("test_withtimeout.jl")

end  # module TerminatorsTests


text:
module Instances

using Niancat.Persistence
using Niancat.Gameface
using Niancat.Http
using Niancat.GameServiceImpl

mutable struct GameInstances
    factory::Ref{Function}
    initialstate::String
    game::Union{Game, Nothing}
    httpclient::HttpClient
    persistence::GamePersistence

    function GameInstances(httpclient::HttpClient, persistence::GamePersistence)
        new(Ref{Function}(), "", nothing, httpclient, persistence)
    end
end

function registergame!(gi::GameInstances, _name::String, initialstate::String, factory::Function)
    gi.factory[] = factory
    gi.initialstate = initialstate
end

function makegameservice(
    gi::GameInstances,
    description::GameInstanceDescription) :: GameService

    ConcreteGameService(description.instanceinfo, description.gameinstanceid, gi.persistence, gi.httpclient)
end

function loadgameinstances!(gi::GameInstances)
    instancedescriptions = getgameinstances(gi.persistence)

    if !isempty(instancedescriptions)
        description = first(instancedescriptions)
        f = gi.factory[]
        gi.game = f(description.state, makegameservice(gi, description))
    end
end

getgame(gi::GameInstances, _gamename::String, _instancename::String) = gi.game

function newgame!(gi::GameInstances, instancename::String, gamename::String)
    # TODO I could imagine having two versions of Niancat, one Swedish and one English.
    # They would be registered under different names, but have the same game name in the
    # database, just different initial states.
    # Then we shouldn't use gamename here, but the underlying database game name.
    description = newgameinstance!(gi.persistence, instancename, gamename, gi.initialstate)
    f = gi.factory[]
    gi.game = f(description.state, makegameservice(gi, description))
end

export GameInstances,  getgame, registergame!, newgame!

end


text:
include("wrapper/exprs.jl")
include("wrapper/type_conversions.jl")
include("wrapper/naming_conventions.jl")
include("wrapper/conventions.jl")
include("wrapper/docs.jl")
include("wrapper/wrap.jl")

include("wrapper/dependency_resolution.jl")
include("wrapper/write.jl")

export
    # Naming Conventions
    ### Convention types
    CamelCaseLower,
    CamelCaseUpper,
    SnakeCaseLower,
    SnakeCaseUpper,

    ### Convention utilities
    detect_convention,
    nc_convert,
    remove_parts,
    remove_prefix,

    # Expr
    name,
    category,
    deconstruct,
    reconstruct,
    rmlines,
    striplines,
    unblock,
    prettify,
    concat_exs,
    broadcast_ex,
    is_broadcast,

    # Wrapping
    VulkanWrapper,
    wrap,
    add_constructor,
    add_constructors,
    extend_from_vk,
    extend_handle_constructor,
    to_expr,
    document


text:
# This file contains the implementation of the Reference Region Model
# and Constrained Reference Region Model, which are to be used with non-linear
# least-squares fitting

@everywhere function RRM(t::Vector{Float64}, p::Vector{Float64}, Crr::Vector{Float64})
  # Reference Region Model
  # Tofts modle using a reference tissue curve (Crr) instead of AIF (Cp)
	# Ct = kR * Crr + kR * (kepRR-kep) * conv(Crr,exp(-kep*t))
	# Reference: Yankeelov et al. (2005), MRI, 23(4), 519–29. doi:10.1016/j.mri.2005.02.013
  kR = p[1] # relative kTrans
  kepRR = p[2] # kep of reference tissue
  kep = p[3] # kep of TOI
  Ct = zeros(t)
  for k in 1:length(Ct)
    tk = t[k]
    @simd for j = 1:k
      @inbounds y = exp(kep*(t[j] - tk)) * Crr[j]
      @inbounds Ct[k] += ifelse((j == 1 || j == k) && k > 1, 0.5*y, y)
    end
  end
  dtp = (t[2] - t[1])
  @simd for k = 1:length(Ct)
    @inbounds Ct[k] = kR * (Crr[k] + dtp*(kepRR-kep)*Ct[k])
  end
  Ct
end

@everywhere function CRRM(tmpT::Vector{Float64}, p::Vector{Float64}, Crr::Vector{Float64})
  # Constrained Reference Region Model
  # Tofts modle using a reference tissue curve (Crr) instead of AIF (Cp)
	# Ct = kR * Crr + kR * (kepRR-kep) * conv(Crr,exp(-kep*t)) --- but kepRR is held fixed
	# Reference: Ahmed & Levesque. (2016), MRM, http://doi.org/10.1002/mrm.26530
  kR = p[1] # relative kTrans
  kep = p[2] # kep of TOI
	# This part is ugly. We needed to pass an extra fixed parameter (i.e. kepRR)
	# but I didn't want to modify the nlsfit() function...
	# ... So, I snuck kepRR as the last element of time.
  kepRR = tmpT[end] # use last element of time as kepRR (because that's how we smuggled it into this function)
  t=tmpT[1:end-1] # delete that last element from time because it doesn't represent an actual time point
  # The rest of the function is similar to RRM()
  Ct = zeros(t)
  for k in 1:length(Ct)
    tk = t[k]
    @simd for j = 1:k
      @inbounds y = exp(kep*(t[j] - tk)) * Crr[j]
      @inbounds Ct[k] += ifelse((j == 1 || j == k) && k > 1, 0.5*y, y)
    end
  end
  dtp = (t[2] - t[1])
  @simd for k = 1:length(Ct)
    @inbounds Ct[k] = kR * (Crr[k] + dtp*(kepRR-kep)*Ct[k])
  end
  Ct
end

max_stars_repo_name: MPUmri/CRRM

max_issues_repo_name: MPUmri/CRRM

max_forks_repo_name: MPUmri/CRRM


text:
struct IndirectSolver <: LinearSolver end

scsint_t(::Type{IndirectSolver}) = Int

function scs_set_default_settings(::Type{IndirectSolver}, data::ScsData{Int})
    return ccall(
        (:scs_set_default_settings, indirect),
        Cvoid,
        (Ptr{Cvoid},),
        data,
    )
end

function scs_init(
    ::Type{IndirectSolver},
    data::ScsData{Int},
    cone::ScsCone{Int},
    info::ScsInfo{Int},
)
    return ccall(
        (:scs_init, indirect),
        Ptr{Cvoid},
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        data,
        cone,
        info,
    )
end

function scs_solve(
    ::Type{IndirectSolver},
    p_work::Ptr{Cvoid},
    data::ScsData{Int},
    cone::ScsCone{Int},
    solution::ScsSolution,
    info::ScsInfo{Int},
)
    return ccall(
        (:scs_solve, indirect),
        Int,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        p_work,
        data,
        cone,
        solution,
        info,
    )
end

function scs_finish(::Type{IndirectSolver}, p_work::Ptr{Cvoid})
    return ccall((:scs_finish, indirect), Cvoid, (Ptr{Cvoid},), p_work)
end


text:
using ModelingToolkit, StaticArrays, LinearAlgebra
using DiffEqBase
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ * (y - x),
    D(y) ~ -z - y,
    D(z) ~ y - β * z]

@test ModelingToolkit.islinear(@named sys = ODESystem(eqs))

eqs2 = [D(x) ~ σ * (y - x),
    D(y) ~ -z - 1 / y,
    D(z) ~ y - β * z]

@test !ModelingToolkit.islinear(@named sys = ODESystem(eqs2))

eqs3 = [D(x) ~ σ * (y - x),
    D(y) ~ -z - y,
    D(z) ~ y - β * z + 1]

@test ModelingToolkit.isaffine(@named sys = ODESystem(eqs))


text:
using HopfNormalForm
using Test

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


text:
# Computes s_n vector from Luger et al. (2018) to ~machine
# precision for b+r > 1:
include("s2.jl")
include("area_triangle.jl")

#using GSL

function vector_sum(x::Array{T,1},y::Array{T,1},n::Int64) where {T <: Real}
# Trying to stabilize arithmetic:
sum_pos = zero(typeof(x[1]))
sum_neg = zero(typeof(x[1]))
# Sort from smallest to largest:
ind = sortperm(abs.(x.*y))
#println(x[ind].*y[ind])
# Carry out sum from smallest to largest (keeping track of positive and negative
# terms separately):
for i=1:n
#  term = x[i]*y[i]
  term = x[ind[i]]*y[ind[i]]
  if term > 0.0
    sum_pos += term
  else
    sum_neg += term
  end
end
# Add together the results:
return sum_pos+sum_neg
end

function aiuv(delta::T,u::Int64,v::Int64) where {T <: Real}
# Computes the double-binomial coefficients A_{i,u,v}:
a=zeros(typeof(delta),u+v+1)
for i=0:u+v
  j1 = maximum([0,u-i])
#  coeff = binomial(u,j1)*binomial(v,u+v-i-j1)*(-1.)^(u+j1)*delta^(v+u-i-j1)
#  a[i+1] = coeff
  for j=j1:minimum([u+v-i,u])
#  for j=j1+1:minimum([u+v-i,u])
#    coeff *= -(u-j+1)*(u+v-i-j+1)/(j*(i+j-u)*delta)
    a[i+1] += binomial(u,j)*binomial(v,u+v-i-j)*(-1.)^(u+j)*delta^(v+u-i-j)
#    a[i+1] += coeff
  end
end
return a
end

function Iv_series(k2::T,v::Int64) where {T <: Real}
# Use series expansion to compute I_v:
nmax = 100
n = 1; error = Inf; if k2 < 1; tol = eps(k2); else; tol = eps(inv(k2)); end
# Computing leading coefficient (n=0):
coeff = 2/(2v+1)
# Add leading term to I_v:
Iv = one(k2)*coeff
# Now, compute higher order terms until desired precision is reached:
while n < nmax && abs(error) > tol
  coeff *= (2.0*n-1.0)*.5*(2n+2v-1)/(n*(2n+2v+1))*k2
  Iv += coeff
  error = coeff/Iv
  n +=1
end
return Iv*k2^v*sqrt(k2)
end

# Compute I_v with hypergeometric function (this requires GSL library,
# which can't use BigFloat or ForwardDiff types):
function Iv_hyp(k2::T,v::Int64) where {T <: Real}
a = 0.5*one(k2); b=v+0.5*one(k2); c=v+1.5*one(k2);  fac = 2/(1+2v)
return sqrt(k2)*k2^v*fac*hypergeom([a,b],c,k2)
end

# Compute J_v with hypergeometric function:
function Jv_hyp(k2::T,v::Int64) where {T <: Real}
if k2 < 1
  a = 0.5; b=v+0.5; c=v+3.0;  fac = 3pi/(4*(v+1)*(v+2))
  for i=1:v
    fac *= (1.-.5/i)
  end
  return sqrt(k2)*k2^v*fac*hypergeom([a,b],c,k2)
else # k^2 >=1
  k2inv = inv(k2)
  # Found a simpler expression than the one in paper (and perhaps more stable for large k^2):
  return  sqrt(pi)*gamma(v+.5)*(sf_hyperg_2F1_renorm(-.5,v+.5,v+1.,k2inv)-(.5+v)*k2inv*sf_hyperg_2F1_renorm(-.5,v+1.5,v+2.,k2inv))
end
end

function Jv_series(k2::T,v::Int64) where {T <: Real}
# Use series expansion to compute J_v:
nmax = 100
n = 1; error = Inf; if k2 < 1; tol = eps(k2); else; tol = eps(inv(k2)); end
# Computing leading coefficient (n=0):
#coeff = 3pi/(2^(2+v)*factorial(v+2))
if k2 < 1
  coeff = 3pi/(2^(2+v)*exp(lfact(v+2)))
# multiply by (2v-1)!!
  for i=1:v
    coeff *= 2.*i-1
  end
# Add leading term to J_v:
  Jv = one(k2)*coeff
# Now, compute higher order terms until desired precision is reached:
  while n < nmax && abs(error) > tol
    coeff *= (2.0*n-1.0)*(2.0*(n+v)-1.0)*.25/(n*(n+v+2))*k2
    Jv += coeff
    error = coeff/Jv
    n +=1
  end
  return Jv*k2^v*sqrt(k2)
else # k^2 >= 1
  coeff = convert(typeof(k2),pi)
  # Compute (2v-1)!!/(2^v v!):
  for i=1:v
    coeff *= 1.-.5/i
  end
  Jv = one(k2)*coeff; n=1
  while n < nmax && abs(error) > tol
    coeff *= (1.-2.5/n)*(1.-.5/(n+v))/k2
    Jv += coeff
    error = coeff/Jv
    n +=1
  end
  return Jv
end
end

function IJv_raise!(v_max::Int64,k2::T,kc::T,Iv::Array{T,1},Jv::Array{T,1})  where {T <: Real}
# This function needs debugging. [ ]
# Compute I_v, J_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate upwards in v:
v = v_max
# Compute I_v via upward iteration on v:
if k2 < 1
# First, compute value for v=0:
  if k2 < 0.5
    Iv[1] = 2*asin(sqrt(k2))
#  Iv[1] = acos(1.-2k2)
  else
    Iv[1] = 2*acos(kc)
  end
# Try something else:
#  Iv[1] = asin(2*k*kc)
# Next, iterate upwards in v:
#  f0 = kc/k
  f0 = kc*k
  v = 1
# Loop over v, computing I_v and J_v from higher v:
  while v <= v_max
    Iv[v+1]=((2v-1)*Iv[v]/2-f0)/v
    f0 *= k2
    v += 1
  end
else # k^2 >= 1
  # Compute v=0
  Iv[1] = pi
  for v=1:v_max
    Iv[v+1]=Iv[v]*(1.0-0.5/v)
  end
end
# Need to compute J_v for v=0, 1:
v= 0
if k2 < 1
  # Use cel_bulirsch:
  if k2 > 0
    Jv[v+1]=2/(3k2*k)*cel_bulirsch(k2,kc,one(k2),k2*(3k2-1),k2*(1-k2))
    Jv[v+2]= 2/(15k2*k)*cel_bulirsch(k2,kc,one(k2),2k2*(3k2-2),k2*(4-7k2+3k2*k2))
  else
    Jv[v+1]= 0.0
    Jv[v+2]= 0.0
  end
else # k^2 >=1
  k2inv = inv(k2)
  Jv[v+1]=2/3*cel_bulirsch(k2inv,kc,one(k2),3-k2inv,3-5k2inv+2k2inv^2)
  Jv[v+2]=cel_bulirsch(k2inv,kc,one(k2),12-8*k2inv,2*(9-8k2inv)*(1-k2inv))/15
end
v=2
while v <= v_max
  Jv[v+1] = (2*(v+1+(v-1)*k2)*Jv[v]-k2*(2v-3)*Jv[v-1])/(2v+3)
  v += 1
end
return
end

function IJv_tridiag!(v_max::Int64,k2::T,kc::T,Iv::Array{T,1},Jv::Array{T,1})  where {T <: Real}
# Compute I_v, J_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate downwards in v:
v = v_max
# Add in k2 > 1 cases [ ]
# First, compute approximation for large v:
#Iv[v+1]=Iv_hyp(k2,v)
if k2 < 1
  Iv[v+1]=Iv_series(k2,v)
# Next, iterate downwards in v:
  f0 = k2^(v-1)*k*kc
# Loop over v, computing I_v and J_v from higher v:
  while v >= 1
    Iv[v] = 2/(2v-1)*(v*Iv[v+1]+f0)
    f0 /= k2
    v -= 1
  end
else # k^2 >= 1
  # Compute v=0 (no need to iterate downwards in this case):
  Iv[1] = pi
  for v=1:v_max
    Iv[v+1]=Iv[v]*(1-.5/v)
  end
end
# Try tri-diagonal solver
# Need to compute J_v for v=0 and v=v_max:
if k2 < 1
  # Use cel_bulirsch:
#  println("k2: ",k2)
  if k2 > 0
#    fe = 2*(2k2-1); fk = (1-k2)*(2-3k2)
#    Jv[1] = 2/(3k2*k)*cel_bulirsch(k2,kc,one(k2),fk+fe,fk+fe*(1-k2))
    Jv[1] = 2/(3k2*k)*cel_bulirsch(k2,kc,one(k2),k2*(3k2-1),k2*(1-k2))
  else
    Jv[1] = 0.0
  end
else # k^2 >=1
  k2inv = inv(k2)
#  fe = 2*(2-k2inv); fk=-1+k2inv
#  Jv[1]=2/3*cel_bulirsch(k2inv,kc,one(k2),fk+fe,fk+fe*(1-k2inv))
  Jv[1]=2/3*cel_bulirsch(k2inv,kc,one(k2),3-k2inv,3-5k2inv+2k2inv^2)
end
Jv[v_max+1]=Jv_series(k2,v_max)
# Now, implement tridiagonal algorithm:
c = zeros(typeof(k2),v_max-1)
d = zeros(typeof(k2),v_max-1)
# Iterate upwards in v (lower):
v = 2
fac = 2*((v+1)+(v-1)*k2)
c[1] = -(2v+3)/fac
d[1] =  (2v-3)*k2/fac*Jv[1]
for v=3:v_max-1
#  f2 = k2*(2v-3); f1 = 2*(v+1+(v-1)*k2)/f2; f3 = (2v+3)/f2
  fac = 2*((v+1)+(v-1)*k2); den = fac + (2v-3)*k2*c[v-2]
  c[v-1] = -(2v+3)/den
  d[v-1] =  (2v-3)*k2*d[v-2]/den
end
v = v_max
fac = 2*((v+1)+(v-1)*k2); den = fac + (2v-3)*k2*c[v-2]
d[v_max-1]=((2v+3)*Jv[v_max+1]+(2v-3)*k2*d[v-2])/den
# Now, back-substitution:
Jv[v_max]=d[v_max-1]
for v=v_max-1:-1:2
  Jv[v]=d[v-1]-c[v-1]*Jv[v+1]
end
return
end

function IJv_lower!(v_max::Int64,k2::T,kc::T,Iv::Array{T,1},Jv::Array{T,1})  where {T <: Real}
# Compute I_v, J_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate downwards in v:
v = v_max
# Add in k2 > 1 cases [ ]
# First, compute approximation for large v:
#Iv[v+1]=Iv_hyp(k2,v)
if k2 < 1
  Iv[v+1]=Iv_series(k2,v)
# Next, iterate downwards in v:
  f0 = k2^(v-1)*k*kc
# Loop over v, computing I_v and J_v from higher v:
  while v >= 1
    Iv[v] = 2/(2v-1)*(v*Iv[v+1]+f0)
    f0 /= k2
    v -= 1
  end
else # k^2 >= 1
  # Compute v=0 (no need to iterate downwards in this case):
  Iv[1] = pi
  for v=1:v_max
    Iv[v+1]=Iv[v]*(1-.5/v)
  end
end
v= v_max
# Need to compute top two for J_v:
#if typeof(k2) == BigFloat
  Jv[v]=Jv_series(k2,v-1); Jv[v+1]=Jv_series(k2,v)
#else
#  Jv[v]=Jv_hyp(k2,v-1); Jv[v+1]=Jv_hyp(k2,v)
#end
#if typeof(k2) == Float64
#  println("v ",v," k2 ",k2," Jv_ser ",convert(Float64,Jv_series(big(k2),v-1))," ",convert(Float64,Jv_series(big(k2),v))," Jv_hyp ",Jv_hyp(k2,v-1)," ",Jv_hyp(k2,v))
#end
# Iterate downwards in v (lower):
while v >= 2
  f2 = k2*(2v-3); f1 = 2*(v+1+(v-1)*k2)/f2; f3 = (2v+3)/f2
  Jv[v-1] = f1*Jv[v]-f3*Jv[v+1]
  v -= 1
end
# Compute first two exactly:
v= 0
if k2 < 1
#  # Use cel_bulirsch:
#  if k2 > 0
#    Jv[v+1]=2/(3*k2*k)*cel_bulirsch(k2,kc,one(k2),k2*(3k2-1),k2*(1-k2))
#    Jv[v+2]=2/(15*k2*k)*cel_bulirsch(k2,kc,one(k2),2*k2*(3k2-2),k2*(4-7k2+3k2*k2))
##    Jv[v+3] =-2/(35*k)*cel_bulirsch(k2,kc,one(k2),8-11k2+k2*k2,(-8+5k2+2k2*k2)*(1-k2))
#  else
#    Jv[v+1]= 0.0
#    Jv[v+2]= 0.0
#  end
else # k^2 >=1
  k2inv = inv(k2)
  Jv[v+1]=2/3*cel_bulirsch(k2inv,kc,one(k2),3-k2inv,3-5k2inv+2k2inv^2)
  Jv[v+2]=cel_bulirsch(k2inv,kc,one(k2),12-8*k2inv,2*(9-8k2inv)*(1-k2inv))/15
#  Jv[v+3]=2/(35*k2)*cel_bulirsch(k2inv,kc,one(k2),-k2^2+11k2-8,(k2^2+16*k2-16)*(1-k2inv))
end
return
end

function Huv_down(l_max,kappa)
# Iterate downwards in v and upwards in u:
Huv = zeros(typeof(kappa),l_max+3,3l_max+3)
kap2 = 0.5*kappa
ckap = cos(kap2); skap = sin(kap2)
skap2 = skap*skap
ckap2 = 1.0-skap2; ckapn = ckap; skapn = skap
u=0
if isodd(l_max)
  v = 3l_max+1
else
  v = 3l_max+2
end
# First, compute approximation for large v:
Huv[u+1,v+1]=kap2^(v+1)*(2/(v+1)-(u+v/3)/(v+3)*kap2^2+
          (5*(3u+v)^2-30u-2v)/360*kap2^4)
#println("Huv: ",Huv[u+1,v+1]," error: ",last_term)
#println("Huv: ",Huv[u+1,v+1])
# Next, iterate downwards in v:
while v >= 2
  Huv[u+1,v-1] = (u+v)/(v-1)*(Huv[u+1,v+1]+2/(u+v)*skap^(v-1)*ckap^(u+1))
  v -= 2
end
# Then, iterate upwards in u:
for u=2:2:l_max+2
  skapn = skap
  v = 0
  while v <= 3l_max+2
    Huv[u+1,v+1]= (2*ckapn*skapn+(u-1)*Huv[u-1,v+1])/(u+v)
    # We only need even terms:
    skapn *= skap2
    v +=2
  end
  ckapn *= ckap2
end
return Huv
end

function s_n!(l_max::Int64,r::T,b::T,sn::Array{T,1}) where {T <: Real}
Jv_check = false # Checks values of J_v
@assert(r > 0.0) # if r=0, then no occultation - can just use phase curve term.
# Computes the s_n terms up to l_max
# Find n_max:
n_max = l_max^2+2*l_max
u_max = floor(Int64,l_max/2+1)
v_max = l_max + 3
if r >= 1+b
  # full obscuration - return zeros
  return sn
end
if b == 0.0
  # Annular eclipse - integrate around the full boundary of both bodies:
  lam = pi/2; slam = one(r); clam = zero(r); clam2 = zero(r)
  phi = pi/2; sphi = one(r); cphi = zero(r); cphi2 = zero(r)
  k2 = Inf; kc = 0.0
else
  if b > abs(1.0-r) && b < 1.0+r
#    lam = asin((1.0-r^2+b^2)/(2*b)); slam = (1.0-r^2+b^2)/(2*b); clam = sqrt((1-(b-r)^2)*((b+r)^2-1))/(2b)
#    slam = (1.0-r^2+b^2)/(2*b); lam = asin(slam);  clam = sqrt((1-(b-r)^2)*((b+r)^2-1))/(2b)
#    slam = (1.0-r^2+b^2)/(2*b);  clam = sqrt((1-(b-r)^2)*((b+r)^2-1))/(2b);  clam2 = (1-(b-r)^2)*((b+r)^2-1)/(4b^2); lam = acos(clam); if lam > pi/2; lam -= pi; end
#    slam = ((1.0-r)*(1.0+r)+b^2)/(2*b);  clam = sqrt((1-b+r)*(1+b-r)*(b+r-1)*(b+r+1))/(2b);  clam2 = (1-b+r)*(1+b-r)*(b+r-1)*(b+r+1)/(4b^2); lam = acos(clam); if slam < 0.; lam = -lam; end
    slam = ((1.0-r)*(1.0+r)+b^2)/(2*b);  clam = 2*area_triangle(1.,b,r)/b;  clam2 = (1-b+r)*(1+b-r)*(b+r-1)*(b+r+1)/(4b^2); lam = acos(clam); if slam < 0.; lam = -lam; end

#    slam = ((1.0-r)*(1.0+r)+b^2)/(2*b);  clam = sqrt((1-b+r)*(1+b-r)*(b+r-1)*(b+r+1))/(2b);  clam2 = (1-b+r)*(1+b-r)*(b+r-1)*(b+r+1)/(4b^2); lam = asin(slam)
#    slam = ((1.0-r)*(1.0+r)+b^2)/(2*b);  lam = asin(slam); clam = cos(lam); clam2=clam^2
    sphi = (1.0-r^2-b^2)/(2*b*r); phi = asin(sphi); cphi = clam/r; cphi2 = clam2/r^2
  else
    lam=pi/2; slam = one(r); clam = zero(r); clam2 = zero(r)
    phi=pi/2; sphi = one(r); cphi = zero(r); cphi2 = zero(r)
  end
# Next, compute k^2 = m:
  k2 = (1.0-(b-r)^2)/(4b*r);
#  if k2 > 1
#    kc = sqrt((1-(b+r)^2)/(1-(b-r)^2))
#  else
#    kc = sqrt(((b+r)^2-1)/(4*b*r))
#  end
  if k2 > 1
    if k2 > 2.0
      kc = sqrt(1.-inv(k2))
    else
      kc2 = (1-(b+r)^2)/(1-(b-r)^2)
      kc = sqrt(kc2)
#      k2 = inv(1-kc2)
    end
  else
    if k2 > 0.5
      kc2 = ((b+r)^2-1)/(4*b*r)
      kc = sqrt(kc2)
#      k2 = 1-kc2
    else
      kc = sqrt(1.-k2)
    end
  end
end

function Hv_raise!(l_max::Int64,k2::T,kc::T,Hv::Array{T,1})  where {T <: Real}
# This Iv function is defined for Q(G_n)/H_uv calculation:
#Iv = zeros(typeof(k2),v_max+1)
# Compute H_uv,  0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate upwards in v:
v_max = l_max+3; v = v_max
# Compute I_v via upward iteration on v:
if k2 < 1
# First, compute value for v=0:
  Hv[1] = 2*asin(k)
# Next, iterate upwards in v:
#  f0 = kc/k
  f0 = kc*k
  v = 1
# Loop over v, computing I_v and J_v from higher v:
  while v <= v_max
    Hv[v+1]=((2v-1)*Hv[v]/2-f0)/v
    f0 *= k2
    v += 1
  end
else # k^2 >= 1
  # Compute v=0
  Hv[1] = pi
  for v=1:v_max
    Hv[v+1]=Hv[v]*(1-.5/v)
  end
end
return
end

function Hv_lower!(l_max::Int64,k2::T,kc::T,Hv::Array{T,1})  where {T <: Real}
# Compute H_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate downwards in v:
v_max = l_max+3; v = v_max
# Add in k2 > 1 cases [ ]
# First, compute approximation for large v:
if k2 < 1
  Hv[v+1]=Iv_series(k2,v)
# Next, iterate downwards in v:
  f0 = k2^(v-1)*k*kc
# Loop over v, computing I_v and J_v from higher v:
  while v >= 1
    Hv[v] = 2/(2v-1)*(v*Hv[v+1]+f0)
    f0 /= k2
    v -= 1
  end
else # k^2 >= 1
  # Compute v=0 (no need to iterate downwards in this case):
  Hv[1] = pi
  for v=1:v_max
    Hv[v+1]=Hv[v]*(1-.5/v)
  end
end
return
end

# First, compute Huv:
#Hv = zeros(typeof(r),v_max+1)
#k2H = ((b+1)^2-r^2)/(4b)
#if k2H < 0.5
#  Hv_lower!(l_max,k2H,sqrt(abs((r^2-(1-b)^2))/(4b)),Hv)
#else
#  Hv_raise!(l_max,k2H,sqrt(abs((r^2-(1-b)^2))/(4b)),Hv)
#end
Huv = zeros(typeof(r),l_max+3,l_max+1)
##clam = cos(lam); slam = sin(lam)
#clam2 = clam*clam; clamn = clam; slamn = slam
#for u=0:2:l_max+2
#  if u == 0
#    Huv[1,1]=  2*lam+pi
#    Huv[1,2]= -2*clam
#    slamn = slam
#    v=2
#    while v <= l_max
#      Huv[1,v+1]= (-2*clam*slamn+(v-1)*Huv[1,v-1])/(u+v)
#      slamn *= slam
#      v+=1
#    end
#  else
#    slamn = slam
#    v = 0
##    Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
##    slamn *= slam
##    v = 1
##    Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
##    slamn *= slam
##    v = 2
#    while v <= l_max
##      Huv[u+1,v+1]= (-2*clam*slamn+(v-1)*Huv[u+1,v-1])/(u+v)
#      Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
#      slamn *= slam
#      v+=1
#    end
#    clamn *= clam2
#  end
#end
clam2 = clam*clam; clamn = clam; slamn = slam
for u=0:2:l_max+2
  if u == 0
    Huv[1,1]=  2*lam+pi
    Huv[1,2]= -2*clam
    slamn = slam
    v=2
    while v <= l_max
      Huv[1,v+1]= (-2*clam*slamn+(v-1)*Huv[1,v-1])/(u+v)
      slamn *= slam
      v+=1
    end
  else
    slamn = slam
    v=0
    Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
    slamn *= slam
    v=1
#    Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
    clamn *= clam2
    Huv[u+1,v+1] = -2*clamn/(u+1)
#    println("u: ",u," v: ",v," Huv[u+1,v+1] ",Huv[u+1,v+1]," -2*clamn/(u+1): ", -2*clamn/(u+1))
    slamn = slam
    v=2
    while v <= l_max
      Huv[u+1,v+1]= (-2*clamn*slamn+(v-1)*Huv[u+1,v-1])/(u+v)
      slamn *= slam
      v+=1
    end
  end
end

Iuv = zeros(typeof(r),l_max+3,l_max+1)
cphi2 = cphi*cphi; cphin = cphi; sphin = sphi
for u=0:2:l_max+2
  if u == 0
    Iuv[1,1]=  2*phi+pi
    Iuv[1,2]= -2*cphi
    sphin = sphi
    v=2
    while v <= l_max
      Iuv[1,v+1]= (-2*cphi*sphin+(v-1)*Iuv[1,v-1])/(u+v)
      sphin *= sphi
      v+=1
    end
  else
    sphin = sphi
    v=0
    Iuv[u+1,v+1]= (2*cphin*sphin+(u-1)*Iuv[u-1,v+1])/(u+v)
    sphin *= sphi
    v=1
    cphin *= cphi2
    Iuv[u+1,v+1] = -2*cphin/(u+1)
    sphin = sphi
    v=2
    while v <= l_max
      Iuv[u+1,v+1]= (-2*cphin*sphin+(v-1)*Iuv[u+1,v-1])/(u+v)
      sphin *= sphi
      v+=1
    end
  end
end
#println("Iuv: ",Iuv)

Kuv_old = zeros(typeof(r),l_max+3,l_max+1)
bonr = b/r
for u=0:2:l_max+2
  for v=0:l_max
    for i=0:v
      fac = binomial(v,i)*bonr^(v-i)
      Kuv_old[u+1,v+1] += fac*Iuv[u+1,i+1]
    end
  end
end
#println("Kuv: ",Kuv_old)
#read(STDIN,Char)

#println("gamma: ",asin(-clam))
#Huv = Huv_down(l_max,acos(-slam))

Iv = zeros(typeof(k2),v_max+1); Jv = zeros(typeof(k2),v_max+1)
# This computes I_v for the largest v, and then works down to smaller values:
if k2 > 0
  if k2 < 0.5 || k2 > 2.0
    IJv_lower!(v_max,k2,kc,Iv,Jv)
  else
    IJv_raise!(v_max,k2,kc,Iv,Jv)
  end
#  IJv_tridiag!(v_max,k2,kc,Iv,Jv)
end
#println("v_max: ",v_max)
if Jv_check && typeof(k2) == Float64
# We can't compute Jv_hyp for values of k2 close to one and large values of v.
  for v=0:8
     Jv_tmp = Jv_hyp(k2,v)
     if abs(Jv[v+1]-Jv_tmp) > 1e-5
       println("v: ",v," k2: ",k2," Jv: ",Jv[v+1]," Jv_hyp: ",Jv_tmp)
     end
  end
end
#println("Iv: ",Iv," Jv: ",Jv)
#IJv_raise!(v_max,k2,kc,Iv,Jv)
#IJv_tridiag!(v_max,k2,kc,Iv,Jv)
#Ivr = zeros(typeof(k2),v_max+1); Jvr = zeros(typeof(k2),v_max+1)
#IJv_raise!(v_max,k2,kc,Ivr,Jvr)
#println("Jv lower: ",Jv," Jv raise: ",Jvr," diff: ",Jv-Jvr)
#Ivt = zeros(typeof(k2),v_max+1); Jvt = zeros(typeof(k2),v_max+1)
#IJv_tridiag!(v_max,k2,kc,Ivt,Jvt)
#println("Jv lower: ",Jv," Jv tridiag: ",Jvt," diff: ",Jv-Jvt)
Kuv = zeros(typeof(r),u_max+1,v_max+1)
Luv = zeros(typeof(r),u_max+1,v_max+1,2)
delta = (b-r)/(2r)
l = 0; n = 0; m = 0; pofgn = zero(typeof(r)); qofgn = zero(typeof(r))
#  k^3*(4br)^(3/2) = (1-(2r\delta)^2)^{3/2}:
#Lfac = (1-(2r*delta)^2)^1.5
Lfac = (1-(b-r)^2)^1.5
while n <= n_max
  if n == 2
    sn[n+1],Eofk,Em1mKdm = s2(r,b)
  else
    mu = l-m; nu = l+m; u=0; v=0
    pofgn = zero(typeof(r)); qofgn = zero(typeof(r))
    # Equation for P(G_n) and Q(G_n):
    if (isodd(mu) && isodd(round(Int64,(mu-1)/2))) || (iseven(mu) && isodd(round(Int64,mu/2)))
      # These cases are zero
    else
      # First, get values of u and v:
      if mod(mu,4) == 0
        u = convert(Int64,mu/4)+1; v= convert(Int64,nu/2)
      elseif iseven(l) && mu == 1
        u = convert(Int64,l/2)-1; v=0
      elseif isodd(l) && mu == 1
        u = convert(Int64,(l-1)/2)-1; v=1
      else
        u=convert(Int64,(mu-1)/4); v=convert(Int64,(nu-1)/2)
      end
      # If they haven't been computed yet, compute Kuv, Luv:
      if Kuv[u+1,v+1] == 0.0 && Luv[u+1,v+1,1] == 0.0 && Luv[u+1,v+1,2] == 0.0
        # First, compute double-binomial coefficients:
        a=aiuv(delta,u,v)
#        Kuv[u+1,v+1]   = sum(a[1:u+v+1].*Iv[u+1:2u+v+1])
        Kuv[u+1,v+1]   = vector_sum(a[1:u+v+1],Iv[u+1:2u+v+1],u+v+1)
#        println("Kuv 1 ",sum(a[1:u+v+1].*Iv[u+1:2u+v+1])," Kuv 2 ",Kuv[u+1,v+1])
#        Luv[u+1,v+1,1] = Lfac*sum(a[1:u+v+1].*Jv[u+1:v+2u+1])
        Luv[u+1,v+1,1] = Lfac*vector_sum(a[1:u+v+1],Jv[u+1:v+2u+1],u+v+1)
        if v <= 1
#          Luv[u+1,v+1,2] = Lfac*sum(a[1:u+v+1].*Jv[u+2:v+2u+2])
          Luv[u+1,v+1,2] = Lfac*vector_sum(a[1:u+v+1],Jv[u+2:v+2u+2],u+v+1)
        end
      end
      # Now, compute P(Gn) & Q(Gn):
      if mod(mu,4) == 0
        pofgn = 2*(2r)^(l+2)*Kuv[u+1,v+1]
#        println("u: ",u," v: ",v," mu: ",mu," nu: ",nu," Kuv: ",Kuv_old[2u+1,v+1]," Kuv new: ",2^(l+3)*Kuv[u+1,v+1])
#        pofgn = r^(l+2)*Kuv_old[2u+1,v+1]
# Use alternate form of Kuv (6/11/2018 notes):
#        Kuv_alt = zero(r)
#        coeff = exp(2*lgamma(u+.5)-lgamma(2u+1))
#        Kuv_alt = coeff*delta^v
#        for i=1:v
#          coeff *= (u+i-.5)/(2u+i)
#          Kuv_alt += binomial(v,i)*delta^(v-i)*coeff
#        end
#        println("u ",u," v ",v," Kuv ",Kuv[u+1,v+1]," Kuv_alt ",Kuv_alt)
##        pofgn = 2*(2r)^(l+2)*Kuv_alt
        if iseven(v) || k2 <= 1  # Q is zero for odd v & k^2 > 1
          qofgn = Huv[2u+1,v+1]
#          a=aiuv(-convert(typeof(k2),0.5),u,v)
#          qofgn = 2^(2u+v+1)*sum(a[1:u+v+1].*Hv[u+1:2u+v+1])
#          qofgn = 2^(2u+v+1)*vector_sum(a[1:u+v+1],Hv[u+1:2u+v+1],u+v+1)
#          println("u ",u," v ",v," Huv ",Huv[2u+1,v+1]," Qnew ",qofgn)
        end
      else
        pofgn = Luv[u+1,v+1,1]
        if mu == 1
          pofgn -= 2*Luv[u+1,v+1,2]
        else
          pofgn *= 2
        end
        pofgn *= (2r)^(l-1)
      end
#      elseif iseven(l) && mu == 1
#        pofgn = -(2r)^(l-1)*(2*Luv[u+1,v+1,2]-Luv[u+1,v+1,1])
#      elseif isodd(l) && mu == 1
#        pofgn = -(2r)^(l-1)*(2*Luv[u+1,v+1,2]-Luv[u+1,v+1,1])
#      else
#        pofgn = 2*(2r)^(l-1)*Luv[u+1,v+1,1]
    end
    sn[n+1] = qofgn-pofgn
#    sn[n+1] = -pofgn
#    sn[n+1] = qofgn
  end
  m +=1
  if m > l
    l += 1
    m = -l
  end
  n += 1
end
# Return the vector of coefficients:
#return sn
return
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using Pkg
Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
using Plots, MLCourse, CSV, MLJ, DataFrames, MLJLinearModels, Random, Flux, MLJFlux, MLJXGBoostInterface, MLJDecisionTreeInterface
include("./Data_Processing.jl")

#Load the data
drop_std = generate(option = "drop", std = "true", valid = "false", test = "true");
med_std = generate(option = "med", std = "true", valid = "false", test = "true");
drop = generate(option = "drop", std = "false", valid = "false", test = "true");
med = generate(option = "med", std = "false", valid = "false", test = "true");

#Random Forest Classifier
m_forest_drop = machine(RandomForestClassifier(n_trees = 500), select(drop.train[:,:], Not(:precipitation_nextday)), drop.train.precipitation_nextday) |> fit!
m_forest_med = machine(RandomForestClassifier(n_trees = 1200), select(med.train[:,:], Not(:precipitation_nextday)), med.train.precipitation_nextday) |> fit!
m_forest_drop_std = machine(RandomForestClassifier(n_trees = 500), select(drop_std.train[:,:], Not(:precipitation_nextday)), drop_std.train.precipitation_nextday) |> fit!
m_forest_med_std = machine(RandomForestClassifier(n_trees = 1200), select(med_std.train[:,:], Not(:precipitation_nextday)), med_std.train.precipitation_nextday) |> fit!

auc_forest_drop = MLJ.auc(predict(m_forest_drop, select(drop.test[:,:], Not(:precipitation_nextday))), drop.test.precipitation_nextday) # 0.909
auc_forest_med = MLJ.auc(predict(m_forest_med, select(med.test[:,:], Not(:precipitation_nextday))), med.test.precipitation_nextday) # 0.927
auc_forest_drop_std = MLJ.auc(predict(m_forest_drop_std, select(drop_std.test[:,:], Not(:precipitation_nextday))), drop_std.test.precipitation_nextday) # 0.908
auc_forest_med_std = MLJ.auc(predict(m_forest_med_std, select(med_std.test[:,:], Not(:precipitation_nextday))), med_std.test.precipitation_nextday) # 0.926

#The med data seems to have the best AUC (not a big difference between standardized or not)

function TunedModel_forest(data)
    Random.seed!(2711)
    model = RandomForestClassifier()
    self_tuning_model = TunedModel(model = model,
                                resampling = CV(nfolds = 5),
                                tuning = Grid(),
                                range = range(model, :n_trees, values=100:50:1200),
                                measure = auc)
    self_tuning_mach = machine(self_tuning_model,
                            select(data, Not(:precipitation_nextday)),
                            data.precipitation_nextday) |> fit!
    self_tuning_mach
end

#Test the tuned model on the drop data set and med data set with both standardized and not-standardized data:
rep1 = report(TunedModel_forest(drop_std.train));
rep1.best_history_entry.measurement
rep1.best_model.n_trees
scatter(reshape(rep1.plotting.parameter_values, :), rep1.plotting.measurements, xlabel = "N", ylabel = "AUC")

rep2 = report(TunedModel_forest(med_std.train));
rep2.best_history_entry.measurement
rep2.best_model.n_trees
scatter(reshape(rep2.plotting.parameter_values, :), rep2.plotting.measurements, xlabel = "N", ylabel = "AUC")

rep3 = report(TunedModel_forest(drop.train));
rep3.best_history_entry.measurement
rep3.best_model.n_trees
scatter(reshape(rep3.plotting.parameter_values, :), rep3.plotting.measurements, xlabel = "N", ylabel = "AUC")

rep4 = report(TunedModel_forest(med.train));
rep4.best_history_entry.measurement
rep4.best_model.n_trees
scatter(reshape(rep4.plotting.parameter_values, :), rep4.plotting.measurements, xlabel = "N", ylabel = "AUC")

#Test the best model over the test set with the best hyper-parameters to pick the best one:
best_mach1 = machine(RandomForestClassifier(n_trees = rep1.best_model.n_trees), select(drop_std.train[:,:], Not(:precipitation_nextday)), drop_std.train.precipitation_nextday)|> fit!
best_mach2 = machine(RandomForestClassifier(n_trees = rep2.best_model.n_trees), select(med_std.train[:,:], Not(:precipitation_nextday)), med_std.train.precipitation_nextday)|> fit!
best_mach3 = machine(RandomForestClassifier(n_trees = rep3.best_model.n_trees), select(drop.train[:,:], Not(:precipitation_nextday)), drop.train.precipitation_nextday)|> fit!
best_mach4 = machine(RandomForestClassifier(n_trees = rep4.best_model.n_trees), select(med.train[:,:], Not(:precipitation_nextday)), med.train.precipitation_nextday)|> fit!
losses(best_mach1, select(drop_std.test, Not(:precipitation_nextday)), drop_std.test.precipitation_nextday)
losses(best_mach2, select(med_std.test, Not(:precipitation_nextday)), med_std.test.precipitation_nextday)
losses(best_mach3, select(drop.test, Not(:precipitation_nextday)), drop.test.precipitation_nextday)
losses(best_mach4, select(med.test, Not(:precipitation_nextday)), med.test.precipitation_nextday)

#Write in the submission file with the best machine trained on all data
train, test = generate(option = "med", std = "true", valid = "false", test = "false");
best_mach = machine(RandomForestClassifier(n_trees = rep2.best_model.n_trees), select(train[:,:], Not(:precipitation_nextday)), train.precipitation_nextday) |> fit!

sample = CSV.read(joinpath(@__DIR__, "data", "sample_submission.csv"), DataFrame);
pred = pdf.(predict(best_mach, test), true);
sample.precipitation_nextday = pred;
CSV.write(joinpath(@__DIR__, "results", "RandomForest_submission.csv"), sample);

max_stars_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_issues_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_forks_repo_name: clara-rossignol/Will-it-rain-tomorrow


text:
include("../src/components/emissions_component.jl")

@testset "emissions" begin

    Precision = 1.0e-11
    T = length(MimiDICE2010.model_years)
    f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

    m = Model()

    set_dimension!(m, :time, MimiDICE2010.model_years)

    add_comp!(m, emissions, :emissions)

# Set the parameters that would normally be internal connection from their Excel values
    set_param!(m, :emissions, :YGROSS, read_params(f, "B92:BI92", T))

# Load the rest of the external parameters
    p = dice2010_excel_parameters(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))
    set_param!(m, :emissions, :sigma, p[:sigma])
    set_param!(m, :emissions, :MIU, p[:MIU])
    set_param!(m, :emissions, :etree, p[:etree])

# Run the one-component model
    run(m)

# Extract the generated variables
    CCA     = m[:emissions, :CCA]
    E       = m[:emissions, :E]
    EIND    = m[:emissions, :EIND]

# Extract the true values
    True_CCA    = read_params(f, "B117:BI117", T)
    True_E      = read_params(f, "B109:BI109", T)
    True_EIND   = read_params(f, "B110:BI110", T)

# Test that the values are the same
    @test maximum(abs, CCA .- True_CCA) ≈ 0. atol = Precision
    @test maximum(abs, E .- True_E) ≈ 0. atol = Precision
    @test maximum(abs, EIND .- True_EIND) ≈ 0. atol = Precision

end


text:
function initplayers(lines)
    sep =  findfirst(==(""), lines)
    p1 = parse.(Int, lines[2:sep-1])
    p2 = parse.(Int, lines[sep+2:end])
    return (p1, p2)
end

function oneround!(p1, p2)
    c1, c2 = popfirst!.([p1, p2])
    c1 > c2 ? push!(p1, c1, c2) : push!(p2,c2, c1)
    return (p1, p2)
end

function score(p)
    l = length(p)
    return sum(a*b for (a, b) in zip(p, l:-1:1))
end

function play!(p1, p2)
    while !(any(isempty.([p1, p2])))
        oneround!(p1, p2)
    end
end

function part1(data)
    p1, p2 = initplayers(data)
    play!(p1, p2)
    p = isempty(p2) ? p1 : p2
    return score(p)
end

const data = readlines("day22.txt")

println("Part 1: $(part1(data))")

function oneround2!(p1, p2, state)
    if (p1,p2) ∈ state
        return 1
    else
        push!(state, (copy(p1),copy(p2)))
    end
    if isempty(p1)
        return 2
    elseif isempty(p2)
        return 1
    end
    c1, c2 = popfirst!.([p1, p2])
    if c1 > length(p1) || c2 > length(p2)
        c1 > c2 ? push!(p1, c1, c2) : push!(p2,c2, c1)
    else
        winner = game!(copy(p1[1:c1]), copy(p2[1:c2]))
	winner == 1 ? push!(p1, c1, c2) : push!(p2, c2, c1)
    end
    return nothing
end

function game!(p1, p2)
    state = []
    w = nothing
    while w === nothing
        w = oneround2!(p1, p2, state)
    end
    return w
end

function part2(data)
    p1, p2 = initplayers(data)
    game!(p1, p2)
    p = isempty(p2) ? p1 : p2
    return score(p)
end

println("Part 2: $(part2(data))")


text:
# Given a direction, what is the preceding direction
function preceding_direction(dir::Direction)
    prec_dir = (dir - 1) % 8
    if prec_dir < 0
        7
    else
        prec_dir
    end
end
## A macro is a way of generating a new output expression, given an unevaluated input expression.
## Julia first parses and evaluates the macro, and the processed code produced by the macro is eventually evaluated like an ordinary expression.
macro check_pairs(vals, a, b, c, d)
    return :((round(getvalue($vals[$a]), 2) - round(getvalue($vals[$c]), 2) <= -d_min &&
    round(getvalue($vals[$a]), 2) - round(getvalue($vals[$d]), 2) <= -d_min &&
    round(getvalue($vals[$b]), 2) - round(getvalue($vals[$c]), 2) <= -d_min &&
    round(getvalue($vals[$b]), 2) - round(getvalue($vals[$d]), 2) <= -d_min))
end

function build_model!(model::JuMP.Model, transit_map::InputGraph,
                        faces::Set{Set{GenericEdge}},
                        planarity_constraints::Integer = 0,
                        elen_factor::Integer = 3,
                        bc_factor::Integer = 3,
                        rpos_factor::Integer = 1
                        )

    # model parameters
    station_list = stations(transit_map)
    edge_list = edges(transit_map)
    N = nstations(transit_map)
    M = nedges(transit_map)
    d_min = 1
    # max canvas should be as small as possible for bigM formulation
    max_canvas = sum(map(x -> x.min_length, edge_list)) * 2
    get_degree = x -> out_degree(transit_map, x)
    max_degree = maximum(map(get_degree, station_list))

    # now we need a mapping from a station id to a number from 1:N
    # the same is needed is needed for edges
    # Todo: Make this run in constant time
    # TODO: Also redo all those helper functions
    function get_station_index(station_id)
        p = 1
        for i in 1:N
            if station_list[i].id == station_id
                return p
            end
            p = p + 1
        end
        -1
    end

    function get_neighbors(station)
        # try only the out-degree
        neighbors_from = map(x -> x.to.id, filter(x -> x.from == station && !x.is_single_label_edge, edge_list))
        setdiff(unique(neighbors_from), [station.id])
    end

    # TODO: Is this really correct?
    function sort_node_by_direction(root_node)
        function (target_node)
            t1 = map(x -> x.direction, filter(x -> x.from.id == root_node && x.to.id == target_node, edge_list))
            first(t1)
        end
    end

    # We need a list of incident edges for the bend costs
    indicent_edge_list = incident_edges(transit_map)
    n_incident_edges = length(indicent_edge_list)

    # also for the planarity constraints we need a list of
    # non_incident_edges edges
    non_indicent_edge_list = non_incident_edges(transit_map, faces)
    n_non_incident_edges = length(non_indicent_edge_list)

    # TODO: Only edge variables for non dummy edges / nodes
    # where sensible
    # setup variables
    @variable(model, 0 <= x[1:N] <= max_canvas)
    @variable(model, 0 <= y[1:N] <= max_canvas)
    @variable(model, 0 <= z1[1:N] <= max_canvas)
    @variable(model, 0 <= z2[1:N] <= max_canvas)
    @variable(model, 1 <= elen[1:M] <= max_canvas)
    @variable(model, 0 <= a[1:M, 1:3] <= 1, Bin)

    @variable(model, 0 <= b[1:N, 1:max_degree] <= 1, Bin)

    # bend cost
    @variable(model, 0 <= bc[1:n_incident_edges] <= 7)
    @variable(model, 0 <= bc_s1[1:n_incident_edges] <= 1, Bin)
    @variable(model, 0 <= bc_s2[1:n_incident_edges] <= 1, Bin)

    # the direction variables
    # only add those that are really needed
    function is_edge(i, j)
        u = station_list[i]
        v = station_list[j]
        length(filter(x -> x.from == u && x.to == v || x.to == u && x.from == v, edge_list)) > 0
    end
    @variable(model, 0 <= d[i = 1:N, j = 1:N; is_edge(i, j)] <= 7)

    # relative postions
    @variable(model, 0 <= rpos[1:M] <= 1, Bin)

    # planarity
    if planarity_constraints > 0
        @variable(model, 0 <= g[1:n_non_incident_edges, 0:7] <= 1, Bin)
    end

    # The objective is a weighted sum that needs to be minimized
    @objective(model, Min, elen_factor * sum(elen[i] for i = 1:M) +
                      bc_factor * sum(bc[i] for i = 1:n_incident_edges) +
                      rpos_factor * sum(rpos[i] for i = 1:M))

    # first some general constraints to model the relationship with
    # xy and z variables
    # For each node
    for i in 1:N
        @constraint(model, z1[i] == x[i] + y[i])
        @constraint(model, z2[i] == x[i] - y[i])
    end

    # For each edge
    for i in 1:M

        # We need to ensure that only one sector is chosen per edge
        @constraint(model, sum(a[i, s] for s = 1:3) == 1)

        edge = edge_list[i]
        min_length = edge.min_length
        orig_dir = edge.direction
        u = get_station_index(edge.from.id)
        v = get_station_index(edge.to.id)
        is_single_label_edge = edge.is_single_label_edge
        is_standard_edge = !is_single_label_edge

        # Ensure that the direction variable has the correct value
        # This block also defines the d variables
        # as I did not figure out how to create a selective group of variables
        prec_dir = preceding_direction(orig_dir)
        succ_dir = (orig_dir + 1) % 8
        orig_dir_rev = (orig_dir + 4) % 8
        prec_dir_rev = preceding_direction(orig_dir_rev)
        succ_dir_rev = (orig_dir_rev + 1) % 8
        dir = [prec_dir, orig_dir, succ_dir]
        dir_rev = [prec_dir_rev, orig_dir_rev, succ_dir_rev]
        if is_standard_edge
            @constraint(model, d[u, v] == prec_dir * a[i, 1] +
                orig_dir * a[i, 2] + succ_dir * a[i, 3])
            @constraint(model, d[v, u] == prec_dir_rev * a[i, 1] +
                orig_dir_rev * a[i, 2] + succ_dir_rev * a[i, 3])
        end

        # constraints that ensure that the relative position cost are
        # computed correctly
        if is_standard_edge
            @constraint(model, -8 * rpos[i] <= d[u, v] - orig_dir)
            @constraint(model, 8 * rpos[i] >= d[u, v] - orig_dir)
        end
        # now we need to set the edge length variable
        # elen is an upperbound on the edgelength
        @constraint(model, x[u] - x[v] <= elen[i])
        @constraint(model, -x[u] + x[v] <= elen[i])
        @constraint(model, y[u] - y[v] <= elen[i])
        @constraint(model, -y[u] + y[v] <= elen[i])

        # the following block ensures the octoliniarity
        if is_standard_edge
            for s in 1:3
                @constraint(model, d[u, v] - dir[s] <= 7 * (1 - a[i, s]))
                @constraint(model, d[u, v] - dir[s] >= -7 * (1 - a[i, s]))
                @constraint(model, d[v, u] - dir_rev[s] <= 7 * (1 - a[i, s]))
                @constraint(model, d[v, u] - dir_rev[s] >= -7 * (1 - a[i, s]))
            end
        end

        # For dummy edges for single labels only selected directions are possible
        # direction 2 and 6 are not possible
        if is_single_label_edge
            allowed_dirs = [1, 3, 4, 5, 7]
            ae = @variable(model, [i, allowed_dirs], Bin)
            @constraint(model, sum(ae[i, d] for d = allowed_dirs) <= 1)
            # Do I need to set the reverse edge?
            @constraint(model, d[u, v] == sum(d * ae[i, d] for d = allowed_dirs))
        end

        bM = max_canvas
        # Todo: Extract this into a macro
        if orig_dir == 0
            # direction == 7
            @constraint(model, z1[v] - z1[u] <= bM * (1 - a[i, 1]))
            @constraint(model, z1[v] - z1[u] >= -1 * bM * (1 - a[i, 1]))
            @constraint(model, z2[v] - z2[u] >= -bM * (1 - a[i, 1]) + 2 * min_length)

            # direction == 0
            @constraint(model, y[v] - y[u] <= bM * (1 - a[i, 2]))
            @constraint(model, y[v] - y[u] >= -bM * (1 - a[i, 2]))
            @constraint(model, x[v] - x[u] >= -1 * bM * (1 - a[i, 2]) + min_length)

            # direction == 1
            @constraint(model, z2[v] - z2[u] <= bM * (1 - a[i, 3]))
            @constraint(model, z2[v] - z2[u] >= -1 * bM * (1 - a[i, 3]))
            @constraint(model, z1[v] - z1[u] >= -1 * bM * (1 - a[i, 3]) + 2 * min_length)
        elseif orig_dir == 1
            # direction == 0
            @constraint(model, y[v] - y[u] <= bM * (1 - a[i, 1]))
            @constraint(model, y[v] - y[u] >= -bM * (1 - a[i, 1]))
            @constraint(model, x[v] - x[u] >= -bM * (1 - a[i, 1]) + min_length)

            # direction == 1
            @constraint(model, z2[v] - z2[u] <= 1 * bM * (1 - a[i, 2]))
            @constraint(model, z2[v] - z2[u] >= -1 * bM * (1 - a[i, 2]))
            @constraint(model, z1[v] - z1[u] >= -1 * bM * (1 - a[i, 2]) + 2 * min_length)

            # direction == 2
            @constraint(model, x[v] - x[u] <= bM * (1 - a[i, 3]))
            @constraint(model, x[v] - x[u] >= -bM * (1 - a[i, 3]))
            @constraint(model, y[v] - y[u] >= -bM * (1 - a[i, 3]) + min_length)
        elseif orig_dir == 2
            # direction == 1
            @constraint(model, z2[v] - z2[u] <= 1 * bM * (1 - a[i, 1]))
            @constraint(model, z2[v] - z2[u] >= -1 * bM * (1 - a[i, 1]))
            @constraint(model, z1[v] - z1[u] >= -1 * bM * (1 - a[i, 1]) + 2 * min_length)

            # direction == 2
            @constraint(model, x[u] - x[v] <= bM * (1 - a[i, 2]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 2]))
            @constraint(model, y[v] - y[u] >= -bM * (1 - a[i, 2]) + min_length)

            # direction == 3
            @constraint(model, z1[u] - z1[v] <= 1 * bM * (1 - a[i, 3]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 3]))
            @constraint(model, z2[u] - z2[v] >= -1 * bM * (1 - a[i, 3]) + 2 * min_length)
        elseif orig_dir == 3
            # direction == 2
            @constraint(model, x[u] - x[v] <= bM * (1 - a[i, 1]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 1]))
            @constraint(model, y[v] - y[u] >= -bM * (1 - a[i, 1]) + min_length)

            # direction == 3
            @constraint(model, z1[u] - z1[v] <= 1 * bM * (1 - a[i, 2]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 2]))
            @constraint(model, z2[u] - z2[v] >= -1 * bM * (1 - a[i, 2]) + 2 * min_length)

            # direction == 4
            @constraint(model, y[u] - y[v] <= bM * (1 - a[i, 3]))
            @constraint(model, y[u] - y[v] >= -bM * (1 - a[i, 3]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 3]) + min_length)
        elseif orig_dir == 4
            # direction == 3
            @constraint(model, z1[u] - z1[v] <= 1 * bM * (1 - a[i, 1]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 1]))
            @constraint(model, z2[u] - z2[v] >= -1 * bM * (1 - a[i, 1]) + 2 * min_length)

            # direction == 4
            @constraint(model, y[u] - y[v] <= bM * (1 - a[i, 2]))
            @constraint(model, y[u] - y[v] >= -bM * (1 - a[i, 2]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 2]) + min_length)

            # direction == 5
            @constraint(model, z2[u] - z2[v] <= 1 * bM * (1 - a[i, 3]))
            @constraint(model, z2[u] - z2[v] >= -1 * bM * (1 - a[i, 3]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 3]) + 2 * min_length)
        elseif orig_dir == 5

            # direction == 4
            @constraint(model, y[u] - y[v] <= bM * (1 - a[i, 1]))
            @constraint(model, y[u] - y[v] >= -bM * (1 - a[i, 1]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 1]) + min_length)

            # direction == 5
            @constraint(model, z2[u] - z2[v] <= 1 * bM * (1 - a[i, 2]))
            @constraint(model, z2[u] - z2[v] >= -1 * bM * (1 - a[i, 2]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 2]) + 2 * min_length)

            # direction == 6
            @constraint(model, x[u] - x[v] <= bM * (1 - a[i, 3]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 3]))
            @constraint(model, y[u] - y[v] >= -1 * bM * (1 - a[i, 3]) + min_length)
        elseif orig_dir == 6
            # direction == 5
            @constraint(model, z2[u] - z2[v] <= 1 * bM * (1 - a[i, 1]))
            @constraint(model, z2[u] - z2[v] >= -1 * bM * (1 - a[i, 1]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 1]) + 2 * min_length)

            # direction == 6
            @constraint(model, x[u] - x[v] <= bM * (1 - a[i, 2]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 2]))
            @constraint(model, y[u] - y[v] >= -1 * bM * (1 - a[i, 2]) + min_length)

            # direction == 7
            @constraint(model, z1[u] - z1[v] <= 1 * bM * (1 - a[i, 3]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 3]))
            @constraint(model, z2[v] - z2[u] >= -1 * bM * (1 - a[i, 3]) + 2 * min_length)
        elseif orig_dir == 7
            # direction == 6
            @constraint(model, x[u] - x[v] <= bM * (1 - a[i, 1]))
            @constraint(model, x[u] - x[v] >= -bM * (1 - a[i, 1]))
            @constraint(model, y[u] - y[v] >= -1 * bM * (1 - a[i, 1]) + min_length)

            # direction == 7
            @constraint(model, z1[u] - z1[v] <= 1 * bM * (1 - a[i, 2]))
            @constraint(model, z1[u] - z1[v] >= -1 * bM * (1 - a[i, 2]))
            @constraint(model, z2[v] - z2[u] >= -1 * bM * (1 - a[i, 2]) + 2 * min_length)

            # direction == 0
            @constraint(model, y[u] - y[v] <= bM * (1 - a[i, 3]))
            @constraint(model, y[u] - y[v] >= -bM * (1 - a[i, 3]))
            @constraint(model, x[v] - x[u] >= -1 * bM * (1 - a[i, 3]) + min_length)
        else
            error("Direction not between 0 and 7")
        end
    end

    for i in 1:N
        #continue
        # we also need to add the circular order constraints
        station = station_list[i]
        station_id = station.id
        deg = out_degree(transit_map, station)
        neighbors = get_neighbors(station)
        if deg >= 2 && length(neighbors) >= 2
            neighbor_indexes = map(get_station_index,
                                    sort(neighbors, by = sort_node_by_direction(station_id)))
            @constraint(model, sum(b[i, j] for j = 1:length(neighbors)) == 1)

            # add circle constraints
            for j in 1:length(neighbor_indexes)
                jIdx = neighbor_indexes[j]
                if j + 1 <= length(neighbor_indexes)
                    jIdxP = neighbor_indexes[j + 1]
                else
                    jIdxP = neighbor_indexes[1]
                end
                @constraint(model, d[i, jIdx] <= d[i, jIdxP] - 1 + 8*b[i, j])
            end
        end
    end

    # TODO: Generate this code to reduce duplication
    if planarity_constraints == 1 # with callback
        # last step is to register a callback to preserve planarity
        function check_edge_spacing(cb)
            bM2 = max_canvas + d_min

            #println("Enter callback")
            for i in 1:n_non_incident_edges
                t = non_indicent_edge_list[i]
                e1 = t[1]
                e2 = t[2]
                u1 = get_station_index(e1.from.id)
                v1 = get_station_index(e1.to.id)
                u2 = get_station_index(e2.from.id)
                v2 = get_station_index(e2.to.id)
                if !(@check_pairs(x, u1, v1, u2, v2) ||
                    @check_pairs(z1, u1, v1, u2, v2) ||
                    @check_pairs(y, u1, v1, u2, v2) ||
                    @check_pairs(z2, u1, v1, u2, v2) ||

                    @check_pairs(z2, u2, v2, u1, v1) ||
                    @check_pairs(x, u2, v2, u1, v1) ||
                    @check_pairs(z1, u2, v2, u1, v1) ||
                    @check_pairs(y, u2, v2, u1, v1))
                    #println("Add constraints for edge", i)
                    # now we need to add a constraint that prevents this
                    @constraint(cb, sum(g[i, d] for d = 0:7) >= 1)

                    # direction 0
                    @constraint(cb, x[u1] - x[u2] <= bM2 * (1 - g[i, 0]) - d_min)
                    @constraint(cb, x[u1] - x[v2] <= bM2 * (1 - g[i, 0]) - d_min)
                    @constraint(cb, x[v1] - x[u2] <= bM2 * (1 - g[i, 0]) - d_min)
                    @constraint(cb, x[v1] - x[v2] <= bM2 * (1 - g[i, 0]) - d_min)

                    # direction 1
                    @constraint(cb, z1[u1] - z1[u2] <= bM2 * (1 - g[i, 1]) - d_min)
                    @constraint(cb, z1[u1] - z1[v2] <= bM2 * (1 - g[i, 1]) - d_min)
                    @constraint(cb, z1[v1] - z1[u2] <= bM2 * (1 - g[i, 1]) - d_min)
                    @constraint(cb, z1[v1] - z1[v2] <= bM2 * (1 - g[i, 1]) - d_min)

                    # direction 2
                    @constraint(cb, y[u1] - y[u2] <= bM2 * (1 - g[i, 2]) - d_min)
                    @constraint(cb, y[u1] - y[v2] <= bM2 * (1 - g[i, 2]) - d_min)
                    @constraint(cb, y[v1] - y[u2] <= bM2 * (1 - g[i, 2]) - d_min)
                    @constraint(cb, y[v1] - y[v2] <= bM2 * (1 - g[i, 2]) - d_min)

                    # direction 3
                    @constraint(cb, z2[u2] - z2[u1] <= bM2 * (1 - g[i, 3]) - d_min)
                    @constraint(cb, z2[u2] - z2[v1] <= bM2 * (1 - g[i, 3]) - d_min)
                    @constraint(cb, z2[v2] - z2[u1] <= bM2 * (1 - g[i, 3]) - d_min)
                    @constraint(cb, z2[v2] - z2[v1] <= bM2 * (1 - g[i, 3]) - d_min)

                    # direction 4
                    @constraint(cb, x[u2] - x[u1] <= bM2 * (1 - g[i, 4]) - d_min)
                    @constraint(cb, x[u2] - x[v1] <= bM2 * (1 - g[i, 4]) - d_min)
                    @constraint(cb, x[v2] - x[u1] <= bM2 * (1 - g[i, 4]) - d_min)
                    @constraint(cb, x[v2] - x[v1] <= bM2 * (1 - g[i, 4]) - d_min)

                    # direction 5
                    @constraint(cb, z1[u2] - z1[u1] <= bM2 * (1 - g[i, 5]) - d_min)
                    @constraint(cb, z1[u2] - z1[v1] <= bM2 * (1 - g[i, 5]) - d_min)
                    @constraint(cb, z1[v2] - z1[u1] <= bM2 * (1 - g[i, 5]) - d_min)
                    @constraint(cb, z1[v2] - z1[v1] <= bM2 * (1 - g[i, 5]) - d_min)

                    # direction 6
                    @constraint(cb, y[u2] - y[u1] <= bM2 * (1 - g[i, 6]) - d_min)
                    @constraint(cb, y[u2] - y[v1] <= bM2 * (1 - g[i, 6]) - d_min)
                    @constraint(cb, y[v2] - y[u1] <= bM2 * (1 - g[i, 6]) - d_min)
                    @constraint(cb, y[v2] - y[v1] <= bM2 * (1 - g[i, 6]) - d_min)

                    # direction 7
                    @constraint(cb, z2[u1] - z2[u2] <= bM2 * (1 - g[i, 7]) - d_min)
                    @constraint(cb, z2[u1] - z2[v2] <= bM2 * (1 - g[i, 7]) - d_min)
                    @constraint(cb, z2[v1] - z2[u2] <= bM2 * (1 - g[i, 7]) - d_min)
                    @constraint(cb, z2[v1] - z2[v2] <= bM2 * (1 - g[i, 7]) - d_min)
                end
            end
        end
        addlazycallback(model, check_edge_spacing)
    elseif planarity_constraints == 2 # no callbacks
        bM2 = max_canvas + d_min

        #println("Enter callback")
        for i in 1:n_non_incident_edges
            t = non_indicent_edge_list[i]
            e1 = t[1]
            e2 = t[2]
            u1 = get_station_index(e1.from.id)
            v1 = get_station_index(e1.to.id)
            u2 = get_station_index(e2.from.id)
            v2 = get_station_index(e2.to.id)
            @constraint(model, sum(g[i, d] for d = 0:7) >= 1)

            # direction 0
            @constraint(model, x[u1] - x[u2] <= bM2 * (1 - g[i, 0]) - d_min)
            @constraint(model, x[u1] - x[v2] <= bM2 * (1 - g[i, 0]) - d_min)
            @constraint(model, x[v1] - x[u2] <= bM2 * (1 - g[i, 0]) - d_min)
            @constraint(model, x[v1] - x[v2] <= bM2 * (1 - g[i, 0]) - d_min)

            # direction 1
            @constraint(model, z1[u1] - z1[u2] <= bM2 * (1 - g[i, 1]) - d_min)
            @constraint(model, z1[u1] - z1[v2] <= bM2 * (1 - g[i, 1]) - d_min)
            @constraint(model, z1[v1] - z1[u2] <= bM2 * (1 - g[i, 1]) - d_min)
            @constraint(model, z1[v1] - z1[v2] <= bM2 * (1 - g[i, 1]) - d_min)

            # direction 2
            @constraint(model, y[u1] - y[u2] <= bM2 * (1 - g[i, 2]) - d_min)
            @constraint(model, y[u1] - y[v2] <= bM2 * (1 - g[i, 2]) - d_min)
            @constraint(model, y[v1] - y[u2] <= bM2 * (1 - g[i, 2]) - d_min)
            @constraint(model, y[v1] - y[v2] <= bM2 * (1 - g[i, 2]) - d_min)

            # direction 3
            @constraint(model, z2[u2] - z2[u1] <= bM2 * (1 - g[i, 3]) - d_min)
            @constraint(model, z2[u2] - z2[v1] <= bM2 * (1 - g[i, 3]) - d_min)
            @constraint(model, z2[v2] - z2[u1] <= bM2 * (1 - g[i, 3]) - d_min)
            @constraint(model, z2[v2] - z2[v1] <= bM2 * (1 - g[i, 3]) - d_min)

            # direction 4
            @constraint(model, x[u2] - x[u1] <= bM2 * (1 - g[i, 4]) - d_min)
            @constraint(model, x[u2] - x[v1] <= bM2 * (1 - g[i, 4]) - d_min)
            @constraint(model, x[v2] - x[u1] <= bM2 * (1 - g[i, 4]) - d_min)
            @constraint(model, x[v2] - x[v1] <= bM2 * (1 - g[i, 4]) - d_min)

            # direction 5
            @constraint(model, z1[u2] - z1[u1] <= bM2 * (1 - g[i, 5]) - d_min)
            @constraint(model, z1[u2] - z1[v1] <= bM2 * (1 - g[i, 5]) - d_min)
            @constraint(model, z1[v2] - z1[u1] <= bM2 * (1 - g[i, 5]) - d_min)
            @constraint(model, z1[v2] - z1[v1] <= bM2 * (1 - g[i, 5]) - d_min)

            # direction 6
            @constraint(model, y[u2] - y[u1] <= bM2 * (1 - g[i, 6]) - d_min)
            @constraint(model, y[u2] - y[v1] <= bM2 * (1 - g[i, 6]) - d_min)
            @constraint(model, y[v2] - y[u1] <= bM2 * (1 - g[i, 6]) - d_min)
            @constraint(model, y[v2] - y[v1] <= bM2 * (1 - g[i, 6]) - d_min)

            # direction 7
            @constraint(model, z2[u1] - z2[u2] <= bM2 * (1 - g[i, 7]) - d_min)
            @constraint(model, z2[u1] - z2[v2] <= bM2 * (1 - g[i, 7]) - d_min)
            @constraint(model, z2[v1] - z2[u2] <= bM2 * (1 - g[i, 7]) - d_min)
            @constraint(model, z2[v1] - z2[v2] <= bM2 * (1 - g[i, 7]) - d_min)

        end
    end

    # bend cost constraints
    for i in 1:n_incident_edges
        t = indicent_edge_list[i]
        e1 = t[1]
        e2 = t[2]
        u = get_station_index(e1.from.id)
        v = get_station_index(e1.to.id)
        w = get_station_index(e2.to.id)
        @constraint(model, -1 * bc[i] <= d[u, v] - d[v, w] - 8 * bc_s1[i] + 8 * bc_s2[i])
        @constraint(model, bc[i] >= d[u, v] - d[v, w] - 8 * bc_s1[i] + 8 * bc_s2[i])
    end

    ModelVariables(x, y)
end

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
abstract type AbstractLoadFormulation <: AbstractDeviceFormulation end
abstract type AbstractControllablePowerLoadFormulation <: AbstractLoadFormulation end
struct StaticPowerLoad <: AbstractLoadFormulation end
struct InterruptiblePowerLoad <: AbstractControllablePowerLoadFormulation end
struct DispatchablePowerLoad <: AbstractControllablePowerLoadFormulation end

########################### dispatchable load variables ####################################
function activepower_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
) where {L <: PSY.ElectricLoad}
    add_variable(
        psi_container,
        devices,
        variable_name(ACTIVE_POWER, L),
        false,
        :nodal_balance_active,
        -1.0;
        ub_value = x -> PSY.get_maxactivepower(x),
        lb_value = x -> 0.0,
    )
    return
end

function reactivepower_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
) where {L <: PSY.ElectricLoad}
    add_variable(
        psi_container,
        devices,
        variable_name(REACTIVE_POWER, L),
        false,
        :nodal_balance_reactive,
        -1.0;
        ub_value = x -> PSY.get_maxreactivepower(x),
        lb_value = x -> 0.0,
    )
    return
end

function commitment_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
) where {L <: PSY.ElectricLoad}
    add_variable(psi_container, devices, variable_name(ON, L), true)
    return
end

####################################### Reactive Power Constraints #########################
"""
Reactive Power Constraints on Controllable Loads Assume Constant PowerFactor
"""
function reactivepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
    model::DeviceModel{L, <:AbstractControllablePowerLoadFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {L <: PSY.ElectricLoad}
    time_steps = model_time_steps(psi_container)
    constraint = JuMPConstraintArray(undef, (PSY.get_name(d) for d in devices), time_steps)
    assign_constraint!(psi_container, REACTIVE, L, constraint)

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(atan((PSY.get_maxreactivepower(d) / PSY.get_maxactivepower(d))))
        reactive = get_variable(psi_container, REACTIVE_POWER, L)[name, t]
        real = get_variable(psi_container, ACTIVE_POWER, L)[name, t] * pf
        constraint[name, t] = JuMP.@constraint(psi_container.JuMPmodel, reactive == real)
    end
    return
end

function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
    model::DeviceModel{L, DispatchablePowerLoad},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {L <: PSY.ElectricLoad}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    if !parameters && !use_forecast_data
        constraint_data = Vector{DeviceRange}(undef, length(devices))
        for (ix, d) in enumerate(devices)
            name = PSY.get_name(d)
            ub = PSY.get_activepower(d)
            limits = (min = 0.0, max = ub)
            range_data = DeviceRange(name, limits)
            add_device_services!(range_data, d, model)
            constraint_data[ix] = range_data
        end
        device_range(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE_RANGE, L),
            variable_name(ACTIVE_POWER, L),
        )
        return
    end

    forecast_label = "get_maxactivepower"
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, x -> PSY.get_maxactivepower(x), ts_vector)
        add_device_services!(timeseries_data, d, model)
        constraint_data[ix] = timeseries_data
    end

    if parameters
        device_timeseries_param_ub(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE, L),
            UpdateRef{L}(ACTIVE_POWER, forecast_label),
            variable_name(ACTIVE_POWER, L),
        )
    else
        device_timeseries_ub(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE, L),
            variable_name(ACTIVE_POWER, L),
        )
    end
    return
end

function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
    model::DeviceModel{L, InterruptiblePowerLoad},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {L <: PSY.ElectricLoad}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    if !parameters && !use_forecast_data
        constraint_data = Vector{DeviceRange}(undef, length(devices))
        for (ix, d) in enumerate(devices)
            name = PSY.get_name(d)
            ub = PSY.get_active(d)
            limits = (min = 0.0, max = ub)
            range_data = DeviceRange(name, limits)
            add_device_services!(range_data, d, model)
            constraint_data[ix] = range_data
        end
        device_semicontinuousrange(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE_RANGE, L),
            variable_name(ACTIVE_POWER, L),
            variable_name(ON, L),
        )
        return
    end

    forecast_label = "get_maxactivepower"
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, x -> PSY.get_maxactivepower(x), ts_vector)
        add_device_services!(timeseries_data, d, model)
        constraint_data[ix] = timeseries_data
    end

    if parameters
        device_timeseries_ub_bigM(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE, L),
            variable_name(ACTIVE_POWER, L),
            UpdateRef{L}(ON, forecast_label),
            constraint_name(ON, L),
        )
    else
        device_timeseries_ub_bin(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE, L),
            variable_name(ACTIVE_POWER, L),
            variable_name(ON, L),
        )
    end
    return
end

########################## Addition of to the nodal balances ###############################
function nodal_expression!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
    ::Type{<:PM.AbstractPowerModel},
) where {L <: PSY.ElectricLoad}
    #Run the Active Power Loop
    nodal_expression!(psi_container, devices, PM.AbstractActivePowerModel)
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)
    if use_forecast_data
        forecast_label = "get_maxactivepower"
        peak_value_function = x -> PSY.get_maxreactivepower(x)
    else
        forecast_label = ""
        peak_value_function = x -> PSY.get_reactivepower(x)
    end
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, peak_value_function, ts_vector)
        constraint_data[ix] = timeseries_data
    end
    if parameters
        include_parameters(
            psi_container,
            constraint_data,
            UpdateRef{L}(REACTIVE_POWER, forecast_label),
            :nodal_balance_reactive,
            -1.0,
        )
        return
    else
        for t in model_time_steps(psi_container)
            for device in constraint_data
                add_to_expression!(
                    psi_container.expressions[:nodal_balance_reactive],
                    device.bus_number,
                    t,
                    -device.multiplier * device.timeseries[t],
                )
            end
        end
    end
    return
end

function nodal_expression!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
    ::Type{<:PM.AbstractActivePowerModel},
) where {L <: PSY.ElectricLoad}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)
    if use_forecast_data
        forecast_label = "get_maxactivepower"
        peak_value_function = x -> PSY.get_maxactivepower(x)
    else
        forecast_label = ""
        peak_value_function = x -> PSY.get_activepower(x)
    end
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, peak_value_function, ts_vector)
        constraint_data[ix] = timeseries_data
    end

    if parameters
        include_parameters(
            psi_container,
            constraint_data,
            UpdateRef{L}(ACTIVE_POWER, forecast_label),
            :nodal_balance_active,
            -1.0,
        )
        return
    else
        for t in model_time_steps(psi_container)
            for device in constraint_data
                add_to_expression!(
                    psi_container.expressions[:nodal_balance_active],
                    device.bus_number,
                    t,
                    -device.multiplier * device.timeseries[t],
                )
            end
        end
    end
    return
end

############################## FormulationControllable Load Cost ###########################
function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
    ::Type{DispatchablePowerLoad},
    ::Type{<:PM.AbstractPowerModel},
) where {L <: PSY.ControllableLoad}
    add_to_cost(psi_container, devices, variable_name(ACTIVE_POWER, L), :variable, -1.0)
    return
end

function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{L},
    ::Type{InterruptiblePowerLoad},
    ::Type{<:PM.AbstractPowerModel},
) where {L <: PSY.ControllableLoad}
    add_to_cost(psi_container, devices, variable_name(ON, L), :fixed, -1.0)
    return
end


text:
using SpikingNetworks
@unit
parseParameters"""
N=100
area = 20000*um^2
Cm = 1*ufarad*cm^(-2)*area
gl = 5e-5*siemens*cm^(-2)*area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm^(-2)*area
g_kd = 30*msiemens*cm^(-2)*area
VT = -63*mV
I=[Float64(0.7*nA*i/N) for i=0:N-1]"""

parseEquations"""
dv/dt=(gl*(El-v)-g_na*m^3*h*(v-ENa)-g_kd*n^4*(v-EK)+I)/Cm
dm/dt=0.32*(13-v+VT)/(exp((13-v+VT)/4)-1)*(1-m)-0.28*(v-VT-40)/(exp((v-VT-40)/5)-1)*m
dn/dt=0.032*(15-v+VT)/(exp((15-v+VT)/5)-1)*(1-n)-0.5*exp((10-v+VT)/40)*n
dh/dt=0.128*exp((17-v+VT)/18)*(1-h)-4/(1+exp((40-v+VT)/5))*h"""

net=init"""
N=N
v=1.0*El*ones(N)
m=0.*ones(N)
n=0.*ones(N)
h=1.*ones(N)"""
spike"""
v.>-40
3"""

@time run!(net,100/net.dt)
t,neuron=rasterPlot(net)
using MatlabPlot
mplot(SpikingNetwork.I,map(x->length(x),net.tSpike))


text:
module GeoMakie

using Reexport
using GeometryBasics, Colors, Makie.MakieLayout, Makie, ImageMagick
using GeoJSON

using Makie.FileIO
using Makie.MakieLayout: Formatting
using Makie.MakieLayout.GridLayoutBase
using Makie.MakieLayout.GridLayoutBase: Side
using GeoInterface: GeoInterface, coordinates, AbstractPolygon, AbstractMultiPolygon, features, geometry

import Makie: convert_arguments, convert_attribute
using Makie.DocStringExtensions

# fix conflicts
const AbstractGeometry = GeometryBasics.AbstractGeometry
const Point = Makie.Point
const attributes = Makie.attributes
const volume = Makie.volume
const Mesh = GeometryBasics.Mesh
const Text = Makie.Text

include("conversions.jl")
include("data.jl")

@reexport using Colors, Makie

end # module


text:
#=
mdFnn.jl:

a Julia based alternative for running multi-dimensional recurrence analysis. Based on Dan Mønster's code: https://github.com/danm0nster/mdembedding

requires:
using Pkg
Pkg.add("LinearAlgebra") # diag caculation
Pkg.add("Distances") # matrix pairwise distances
Pkg.add("ElasticArrays") # efficent matrices

using Distances, ElasticArrays, LinearAlgebra, StatsBase;
=#

function mdFnn(data, tau, maxEmb, numSamples, Rtol, Atol)
    if typeof(data) <: Dataset
        data = data |> Matrix
    end

    dims = size(data,2);
    N = size(data,1);

    fnnPerc = [100.0]; # first FNN

    if dims == 1
        Ra = sum(cov(data));
    else
        Ra = sum(diag(cov(data))); # estimate of attractor size
    end

    if (N-tau*(maxEmb-1)) < numSamples # check whether enough data points exist for random sampling
        numSamples = N-tau*(maxEmb-1);
        samps = 1:1:numSamples;
    else
        #numSamples = N-tau*(maxEmb-1);
        samps = sortperm(rand(numSamples));
    end

    embData = ElasticArray{Float64}(undef, N-tau*(maxEmb-1), 0);

    for i = 1:maxEmb # embed data
        append!(embData,data[1+(i-1)*tau:end-(maxEmb-i)*tau, 1:dims]);
        dists =  pairwise(Euclidean(),embData.^2; dims=1);
        r2d1 = Array{Float64}(undef, 0);
        yRd1 = Array{Int}(undef, 0);
        for j = 1:numSamples # get nearest neighbors and distances
            temp = sort(dists[:,samps[j]]);
            coord = sortperm(dists[:,samps[j]]);
            push!(r2d1, temp[2]);
            push!(yRd1, coord[2] |> Int);
        end
        if i == 1
            global r2d = r2d1;
            global yRd = r2d1;
        else
            fnnTemp = Array{Float64}(undef, 0);
            for j = 1:length(r2d1)
                temp = dists[:,samps[j]];
                # check whether neighbors are false
                if sqrt(abs((temp[yRd1[j]] - r2d[j])/r2d[j])) > Rtol || abs(temp[yRd1[j]] - r2d[j])/Ra > Atol
                    append!(fnnTemp, 1);
                else
                    append!(fnnTemp, 0);
                end
            end
            push!(fnnPerc, 100*sum(fnnTemp)/length(fnnTemp)); # compute percentage of FNN
            r2d = r2d1;
            yRd = r2d1;
        end
    end

    return(fnnPerc)
end


text:
using NonstationaryProcesses
using Plots

sim = gaussianBimodalSim(
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = ((5.0,), (1.0,), (0.2)),  # (μ, σ, α)
    transient_t0 = 0.0,
    t0 = 0.0,
    savedt = 1,
    tmax = 1000.0,
    solver_rng=rand(UInt))

S = simulate(sim)

p1 = tsdensity(S, densityoffset=-15)

S = simulate(sim(parameter_profile_parameters=((5.0,), (1.0,), (0.2))))
p2 = tsdensity(S, densityoffset=-15)


text:
export Sphere1Sphere1, PriorSphere1, PackedSphere1Sphere1, PackedPriorSphere1

"""
$(TYPEDEF)
"""
mutable struct Sphere1Sphere1{T} <: IncrementalInference.FunctorPairwise where {T <: SamplableBelief}
  Z::T
  Sphere1Sphere1{T}() where {T <: SamplableBelief} = new{T}()
  Sphere1Sphere1{T}(z1::T) where {T <: SamplableBelief} = new{T}(z1)
end
Sphere1Sphere1(z::T) where {T <: SamplableBelief} = Sphere1Sphere1{T}(z)

getSample(s::Sphere1Sphere1{<: SamplableBelief}, N::Int=1) = (rand(s.Z, N), )

function (s::Sphere1Sphere1{<: SamplableBelief})(res::Vector{Float64},
                                                 userdata::FactorMetadata,
                                                 idx::Int,
                                                 meas::Tuple,
                                                 wxi::Array{Float64,2},
                                                 wxj::Array{Float64,2}  )
  #
  wXjhat = addtheta(wxi[1,idx], meas[1][1,idx])
  res[1] = difftheta(wxj[1,idx], wXjhat)  # jXjhat =
  nothing
end

"""
$(TYPEDEF)

Introduce direct observations on all dimensions of a Sphere1 variable:

Example:
--------
```julia
PriorSphere1( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))) )
```
"""
mutable struct PriorSphere1{T} <: IncrementalInference.FunctorSingleton  where {T <: IncrementalInference.SamplableBelief}
    Z::T
    PriorSphere1{T}() where T = new{T}()
    PriorSphere1{T}(x::T) where {T <: IncrementalInference.SamplableBelief}  = new{T}(x)
end
PriorSphere1(x::T) where {T <: IncrementalInference.SamplableBelief} = PriorSphere1{T}(x)
function PriorSphere1(mu::Array{Float64}, cov::Array{Float64,2}, W::Vector{Float64})
  @warn "PriorSphere1(mu,cov,W) is deprecated in favor of PriorSphere1(T(...)) -- use for example PriorSphere1(MvNormal(mu, cov))"
  PriorSphere1(MvNormal(mu[:], cov))
end
function getSample(p2::PriorSphere1, N::Int=1)
  return (rand(p2.Z,N), )
end

"""
$(TYPEDEF)

Serialized object for storing PriorSphere1.
"""
mutable struct PackedPriorSphere1  <: IncrementalInference.PackedInferenceType
  datastr::String
  PackedPriorSphere1() = new()
  PackedPriorSphere1(x::String) = new(x)
end
function convert(::Type{PackedPriorSphere1}, d::PriorSphere1)
  return PackedPriorSphere1(string(d.Z))
end
function convert(::Type{PriorSphere1}, d::PackedPriorSphere1)
  distr = extractdistribution(d.datastr)
  return PriorSphere1{typeof(distr)}(distr)
end

# --------------------------------------------

"""
$(TYPEDEF)

Serialized object for storing Sphere1Sphere1.
"""
mutable struct PackedSphere1Sphere1  <: IncrementalInference.PackedInferenceType
  datastr::String
  PackedSphere1Sphere1() = new()
  PackedSphere1Sphere1(x::String) = new(x)
end
function convert(::Type{Sphere1Sphere1}, d::PackedSphere1Sphere1)
  return Sphere1Sphere1(extractdistribution(d.datastr))
end
function convert(::Type{PackedSphere1Sphere1}, d::Sphere1Sphere1)
  return PackedSphere1Sphere1(string(d.z))
end

# --------------------------------------------


text:
function matchlexem(item::String)

	if (item == "EOF")
		return Int(EOF)
	end

	#Checks if the lexem has invalid characters
	invalidpart =
	filter(lexem) do x
			!contains(==,vcat(digits_only,chars),x)
	end

	if (legnth(invalidpart > 0))
		return Int(LEX_ERR)
	end

	if (contains(==,digits_only,item[1]) && !checkisnumber(item)) #if the lexem starts with a number, but the lexem itself is not a number
		return Int(LEX_ERR)
	end

	if checkisnumber(item)
		if (checkisint(item))
			return Int(CT_INT)
		else
			return Int(CT_FLOAT)
		end
	end

	if item[1] == '"' && item[length(item)] == '"' #if it starts and ends with "
		return Int(CT_STRING)
	end

	if item[1] == '\'' && item[length(item)] == '\'' #if it starts and ends with \)
		return Int(CT_CHAR)
	end

	#From this point we are sure the token is not a number neither is a string
	return	@match item begin
				"~" => Int(OPR_UN_NEG)
				"[" => Int(O_BRCKT)
				"]" => Int(C_BRCKT)
				"(" => Int(O_PRTSIS)
				")" => Int(C_PRTSIS)
				"," => Int(COMMA)
				";" => Int(SMCL)
				"(" => Int(O_PRTSIS)
				"{" => Int(O_C_BRCKT)
				"}" => Int(C_C_BRCKT)
				"=" => Int(OPR_ATR)
				">" => Int(OPRLR_LGT)
				"<" => Int(OPRLR_LGT)
				">=" => Int(OPRLR_LGT_EQ)
				"<=" => Int(OPRLR_LGT_EQ)
				"int" => Int(IDT_INT)
				"float" => Int(IDT_FLOAT)
				"char" => Int(IDT_CHAR)
				"string" => Int(IDT_STRING)
				"print" => Int(FN_PRINT)
				"main" => Int(FN_MAIN)
				"::" => Int(VEC_IN)  #From this point we are sure the token is not an operator or a reserved symbol
				"EOF" => Int(EOF)
				 _ => false
			end

end

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
module Re

link1 = """HackerRank"""
link2 = """
HackerRank
 """

"""HackerRank"""
ptn = r"(<.+></.+>)"

"""HackerRank"""
ptn = r"(.+)"

@show match(ptn, link1)
@show match(ptn, link2)

ptn = r"(.+)</.+>"

#=
t_ptr = r"(a a)(?:123)"
str = "Aa a123b"
@show r = match(t_ptr,str)
=#

end