text:
using DataDeps

register(DataDep(
    "Stanford Sentiment Treebank",
    """
    Stanford Sentiment Treebank V1.0

    It includes fine grained sentiment labels for 215,154 phrases in the parse trees of 11,855 sentences.

    title       :   Parsing With Compositional Vector Grammars
    author      :   Richard Socher and Alex Perelygin and Jean Wu and Jason Chuang and Christopher Manning and Andrew Ng and Christopher Potts
    booktitle   :   EMNLP
    year        :   2013
    Website     :   https://nlp.stanford.edu/sentiment/index.html

    Cite this paper if this dataset is used:
    https://nlp.stanford.edu/~socherr/EMNLP2013_RNTN.pdf
    """,
    "http://nlp.stanford.edu/~socherr/stanfordSentimentTreebank.zip",
    "3f5209483b46bbf129cacbbbe6ae02fe780407034f61cf6342b7833257c3f1db",
    post_fetch_method =function (fn)
        unpack(fn)
        dir = "stanfordSentimentTreebank"
        files = readdir(dir)
        mv.(joinpath.(dir, files), files)
        rm(dir), rm("__MACOSX", recursive=true)

        local dict
        open("dictionary.txt", "r") do fdict
            dict = split.(split(read(fdict, String), '\n'), '|')[1:end-1]
            dict = permutedims(reshape(hcat(dict...), (length( dict[1]), length(dict))))
        end

        local sentiments
        open("sentiment_labels.txt", "r") do fsenti
            sentiments = split.(split(read(fsenti, String), '\n'), '|')[2:end-1]
            sentiments = permutedims(reshape(hcat(sentiments...), (length(sentiments[1]), length(sentiments))))
        end

        open("sa_dataset.txt", "w") do fdataset
            for id=1:size(dict)[1]
                line = join([dict[id, 1], sentiments[parse(Int, dict[id, 2])+1, 2]], '|')
                write(fdataset, line, "\n")
            end
        end
    end
))


text:
struct LidDrivenCavityFlow{T <: Real, Int <: Integer} <: FluidFlowProblem
    rho_0::T
    u_max::T
    ν::T
    NX::Int
    NY::Int
    domain_size::Tuple{T, T}
end

function LidDrivenCavityFlow(
    ν = 1.0 / 6.0,
    scale = 2,
    NX = 16 * scale,
    NY = NX,
    domain_size = (1.0, 1.0),
)
    # # Since our velocity profile is constant over x it is sufficient to only take
    # # 3 nodes in the x direction
    # NX = 5
    u_max = 0.01 / scale

    return LidDrivenCavityFlow(1.0, u_max, ν, NX, NY, domain_size)
end

function density(
    q::Quadrature,
    problem::LidDrivenCavityFlow,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    return 1.0
end

function pressure(
    q::Quadrature,
    problem::LidDrivenCavityFlow,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    return 1.0
end

function velocity(
    problem::LidDrivenCavityFlow,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    return [
        0.0
        0.0
    ]
end

"""
Apply a moving wall to the top of the domain and bounce back to all other
sides of the domain
"""
boundary_conditions(problem::LidDrivenCavityFlow) = [
    BounceBack(East(), 1:(problem.NX), 1:(problem.NY)),
    BounceBack(South(), 1:(problem.NX), 1:(problem.NY)),
    BounceBack(West(), 1:(problem.NX), 1:(problem.NY)),
    MovingWall(North(), 1:(problem.NX), 1:(problem.NY), [problem.u_max, 0]),
]

ext: jl


text:
using .Nodes:
    NodeNil, is_nil

mutable struct NodeStack
    nodes::Array{Ranger.AbstractNode,1}

    next_node::Ranger.AbstractNode
    running_node::Ranger.AbstractNode

    function NodeStack()
        o = new()

        o.nodes = Array{Ranger.AbstractNode,1}[]
        o.next_node = NodeNil()
        o.running_node = NodeNil()

        o
    end
end

function is_empty(stack::NodeStack)
    length(stack.nodes) == 0
end

function has_next_node(stack::NodeStack)
    !is_nil(stack.next_node)
end

function has_running_node(stack::NodeStack)
    !is_nil(stack.running_node)
end

function set_next_node_nil(stack::NodeStack)
    # Setting node to stack means nil
    stack.next_node = NodeNil();
end

function set_running_node_nil(stack::NodeStack)
    # Setting node to stack means nil
    stack.running_node = NodeNil();
end

function set_running_node(stack::NodeStack)
    stack.running_node = stack.next_node;
end

function push(stack::NodeStack, node::Ranger.AbstractNode)
    stack.next_node = node

    println("Pushing node: ", stack.next_node)

    push!(stack.nodes, node);
end

function pop(stack::NodeStack)
    if !is_empty(stack)
        top = pop!(stack.nodes)
        stack.next_node = top
        println("NodeStack: popped: ", top);
    else
        println("NodeStack -- no nodes to pop.");
    end
end

function replace(stack::NodeStack, replacement::Ranger.AbstractNode)
    stack.next_node = replacement

    # Replacement is the act of popping and pushing. i.e. replacing
    # the stack top with the new node.
    if !is_empty(stack)
        top = pop!(stack.nodes)
        println("NodeStack: popped ", top, " for replacement with ", replacement);
    else
        println("NodeStack: replace WARNING, nothing popped");
    end

    push!(stack.nodes, replacement);
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
############################################
###########################################
##########################################
# Main course:
##########################################
#############################################
##############################################
include("./Params_And_Functions.jl")

NRMSE_RLS = Matrix{Float64}(undef, NoIter,NoExp);
NRMSE_ORRLSl1 = Matrix{Float64}(undef, NoIter,NoExp);
NRMSE_ORRLSMCP = Matrix{Float64}(undef, NoIter,NoExp);
NRMSE_ORRLSCDl1 = Matrix{Float64}(undef, NoIter,NoExp);
NRMSE_ORRLSCDMCP = Matrix{Float64}(undef, NoIter,NoExp);
NRMSE_HRLS_ADMM = Matrix{Float64}(undef, NoIter,NoExp);
NRMSE_HRLS_GIST = Matrix{Float64}(undef, NoIter,NoExp);
NRMSE_HRLS_FMHSDM = Matrix{Float64}(undef, NoIter,NoExp);

TimeRLS = Vector{Float64}(undef, NoExp);
TimeORRLSl1 = Vector{Float64}(undef, NoExp);
TimeORRLSMCP = Vector{Float64}(undef, NoExp);
TimeORRLSCDl1 = Vector{Float64}(undef, NoExp);
TimeORRLSCDMCP = Vector{Float64}(undef, NoExp);
TimeORHRLSADMM = Vector{Float64}(undef, NoExp);
TimeORHRLSGIST = Vector{Float64}(undef, NoExp);
TimeORHRLSFMHSDM = Vector{Float64}(undef, NoExp);

for iEXP in 1:NoExp

    @printf("********************** Experiment %d/%d **********************\n", iEXP, NoExp);
    print("-- Initializing..\n");

    System, Rnoise, invRnoise, MaxEigInvRnoise, VarWhiteNoise, ARmatrix = GenerateModelMatrices(L, ARcoeffU, ARcoeffL, SNR, SystemIsSparse, SystemSparsityLevel);
    Input0, Output0, LastNoiseVec = InputOutput(P,#
                                                L,#
                                                K0,#
                                                System,#
                                                ARmatrix,#
                                                pBernoulli,#
                                                VarOutliers,#
                                                VarWhiteNoise,#
                                                sqrt(VarWhiteNoise)*randn(L));
    print("-- Warm start via Gurobi..\n");
    F0, o0 = WarmStart(P, L, K0+1, 1, invRnoise, lambdaORRLSl1, Input0, Output0);
    F0HRLS, o0HRLS = WarmStart(P, L, K0+1, 2*(K0+1), invRnoise, lambdao, Input0, Output0);
    print("-- Generating input-output data..\n");
    Input, Output, _ = InputOutput(P,#
                                   L,#
                                   NoIter,#
                                   System,#
                                   ARmatrix,#
                                   pBernoulli,#
                                   VarOutliers,#
                                   VarWhiteNoise,#
                                   LastNoiseVec);
    ProblemParams = GenerateProblemParams(NoIter,#
                                          P,#
                                          L,#
                                          SystemIsSparse,#
                                          Input,#
                                          System,#
                                          Output,#
                                          Input0[:, K0+1],#
                                          Output0[:, K0+1],#
                                          Rnoise,#
                                          invRnoise,#
                                          MaxEigInvRnoise,#
                                          eVarPi,#
                                          alpha,#
                                          Radius,#
                                          gamma,#
                                          epsilonORRLS,#
                                          F0,#
                                          o0,#
                                          F0HRLS,#
                                          o0HRLS,#
                                          lambdaF,#
                                          lambdao,#
                                          rho,#
                                          lambda,#
                                          MultiplyGradBy,#
                                          lambdaFMHSDM,#
                                          alphaFMHSDM,#
                                          JFMHSDM,#
                                          UseBatchMethod,#
                                          lambdaORRLSl1,#
                                          rhoADMM,#
                                          JADMM,#
                                          thetaMCP,#
                                          lambdaMCP,#
                                          tminMCP,#
                                          tmaxMCP,#
                                          etaMCP,#
                                          sigmaMCP,#
                                          mBuffMCP,#
                                          JGIST);

    # Methods:

    NRMSE_RLS[:, iEXP], TimeRLS[iEXP] = RLS(ProblemParams);
    @printf("-- RLS: %esecs/iteration.\n", TimeRLS[iEXP]);

    NRMSE_ORRLSl1[:, iEXP], TimeORRLSl1[iEXP] = ORRLSl1(ProblemParams);
    @printf("-- OR-RLS(LASSO): %esecs/iteration.\n", TimeORRLSl1[iEXP]);

    NRMSE_ORRLSMCP[:,iEXP], TimeORRLSMCP[iEXP] = ORRLSMCP(ProblemParams);
    @printf("-- OR-RLS(MCP): %esecs/iteration.\n", TimeORRLSMCP[iEXP]);

    NRMSE_ORRLSCDl1[:,iEXP], TimeORRLSCDl1[iEXP] =  ORRLSCD_L1(ProblemParams);
    @printf("-- OR-RLS(CDl1): %esecs/iteration.\n", TimeORRLSCDl1[iEXP]);

    NRMSE_ORRLSCDMCP[:,iEXP], TimeORRLSCDMCP[iEXP] = ORRLSCD_MCP(ProblemParams);
    @printf("-- OR-RLS(CDMCP): %esecs/iteration.\n", TimeORRLSCDMCP[iEXP]);

    NRMSE_HRLS_ADMM[:,iEXP], TimeORHRLSADMM[iEXP] = HRLS_ADMM(ProblemParams);
    @printf("-- OR-HRLS(ADMM): %esecs/iteration.\n", TimeORHRLSADMM[iEXP]);

    NRMSE_HRLS_GIST[:,iEXP], TimeORHRLSGIST[iEXP] = HRLS_GIST(ProblemParams);
    @printf("-- OR-HRLS(GIST): %esecs/iteration.\n", TimeORHRLSGIST[iEXP]);

    NRMSE_HRLS_FMHSDM[:,iEXP],  TimeORHRLSFMHSDM[iEXP] = HRLS_FMHSDM(ProblemParams);
    @printf("-- OR-HRLS(FMHSDM): %esecs/iteration.\n", TimeORHRLSFMHSDM[iEXP]);

end

NRMSE_RLSav =  Array{Int64}(undef, NoIter)
NRMSE_ORRLSl1av =  Array{Int64}(undef, NoIter)
NRMSE_ORRLSMCPav =  Array{Int64}(undef, NoIter)
NRMSE_ORRLSCDl1av =  Array{Int64}(undef, NoIter)
NRMSE_ORRLSCDMCPav =  Array{Int64}(undef, NoIter)
NRMSE_HRLS_ADMMav =  Array{Int64}(undef, NoIter)
NRMSE_HRLS_GISTav =  Array{Int64}(undef, NoIter)
NRMSE_HRLS_FMHSDMav =  Array{Int64}(undef, NoIter)

NRMSE_RLSav = sum(NRMSE_RLS, dims=2)/NoExp;
NRMSE_ORRLSl1av = sum(NRMSE_ORRLSl1, dims=2)/NoExp;
NRMSE_ORRLSMCPav = sum(NRMSE_ORRLSMCP, dims=2)/NoExp;
NRMSE_ORRLSCDl1av = sum(NRMSE_ORRLSCDl1, dims=2)/NoExp;
NRMSE_ORRLSCDMCPav = sum(NRMSE_ORRLSCDMCP, dims=2)/NoExp;
NRMSE_HRLS_ADMMav = sum(NRMSE_HRLS_ADMM, dims=2)/NoExp;
NRMSE_HRLS_GISTav = sum(NRMSE_HRLS_GIST, dims=2)/NoExp;
NRMSE_HRLS_FMHSDMav = sum(NRMSE_HRLS_FMHSDM, dims=2)/NoExp;

TimeRLSav = sum(TimeRLS)/NoExp;
TimeORRLSl1av = sum(TimeORRLSl1)/NoExp;
TimeORRLSMCPav = sum(TimeORRLSMCP)/NoExp;
TimeORRLSCDl1av = sum(TimeORRLSCDl1)/NoExp;
TimeORRLSCDMCPav = sum(TimeORRLSCDMCP)/NoExp;
TimeORHRLSADMMav = sum(TimeORHRLSADMM)/NoExp;
TimeORHRLSGISTav = sum(TimeORHRLSGIST)/NoExp;
TimeORHRLSFMHSDMav = sum(TimeORHRLSFMHSDM)/NoExp;

Sq_TimeRLS = Vector{Float64}(undef, NoExp);
Sq_TimeORRLSl1 = Vector{Float64}(undef, NoExp);
Sq_TimeORRLSMCP = Vector{Float64}(undef, NoExp);
Sq_TimeORRLSCDl1 = Vector{Float64}(undef, NoExp);
Sq_TimeORRLSCDMCP = Vector{Float64}(undef, NoExp);
Sq_TimeORHRLSADMM = Vector{Float64}(undef, NoExp);
Sq_TimeORHRLSGIST = Vector{Float64}(undef, NoExp);
Sq_TimeORHRLSFMHSDM = Vector{Float64}(undef, NoExp);

for iEXP in 1:NoExp
    Sq_TimeRLS[iEXP] = (TimeRLSav - TimeRLS[iEXP])^2
    Sq_TimeORRLSl1[iEXP] = (TimeORRLSl1av - TimeORRLSl1[iEXP])^2
    Sq_TimeORRLSMCP[iEXP] = (TimeORRLSMCPav - TimeORRLSMCP[iEXP])^2
    Sq_TimeORRLSCDl1[iEXP] = (TimeORRLSCDl1av - TimeORRLSCDl1[iEXP])^2
    Sq_TimeORRLSCDMCP[iEXP] = (TimeORRLSCDMCPav - TimeORRLSCDMCP[iEXP])^2
    Sq_TimeORHRLSADMM[iEXP] = (TimeORHRLSADMMav - TimeORHRLSADMM[iEXP])^2
    Sq_TimeORHRLSGIST[iEXP] = (TimeORHRLSGISTav - TimeORHRLSGIST[iEXP])^2
    Sq_TimeORHRLSFMHSDM[iEXP] = (TimeORHRLSFMHSDMav - TimeORHRLSFMHSDM[iEXP])^2
end

SD_TimeRLS = sum(Sq_TimeRLS)/NoExp
SD_TimeORRLSl1 = sum(Sq_TimeORRLSl1)/NoExp
SD_TimeORRLSMCP = sum(Sq_TimeORRLSMCP)/NoExp
SD_TimeORRLSCDl1 = sum(Sq_TimeORRLSCDl1)/NoExp
SD_TimeORRLSCDMCP = sum(Sq_TimeORRLSCDMCP)/NoExp
SD_TimeORHRLSADMM = sum(Sq_TimeORHRLSADMM)/NoExp
SD_TimeORHRLSGIST = sum(Sq_TimeORHRLSGIST)/NoExp
SD_TimeORHRLSFMHSDM = sum(Sq_TimeORHRLSFMHSDM)/NoExp

# Data = SharedArray{Int64};
Data = [NRMSE_RLSav NRMSE_ORRLSCDl1av NRMSE_ORRLSCDMCPav NRMSE_ORRLSl1av NRMSE_ORRLSMCPav NRMSE_HRLS_ADMMav NRMSE_HRLS_GISTav NRMSE_HRLS_FMHSDMav];

print("-- Plotting..\n");

# gr()
plotlyjs()
plot(1:NoIter, Data, #
     label = ["OR-RLS(RLS)" "OR-RLS(CD-L1)" "OR-RLS(CD-MCP)" "OR-RLS(LASSO)" "OR-RLS(GIST)" "OR-HRLS(ADMM)" "OR-HRLS(GIST)" "OR-HRLS(FMHSDM)"],#
     xlabel = "Iteration/time index",#
     ylabel = "NRMSE",#
     yaxis = :log10,#
     ylims = (-Inf, 10),#
     lw = 2#
     );

# # display(p)
# # print(Data)
# # savefig("fig")

#  open("VarLinParam_1.txt", "w") do io
#            print("###################################################\r\n")
#            print(io, "Problems Parameters\r\n")
#            print(io,"####################################################\r\n")
#            print(io,"SNR (DB):")
#            writedlm(io,SNR)
#            print(io,"No. of experiments:")
#            writedlm(io,NoExp)
#            print(io,"Number of iterations:")
#            writedlm(io,NoIter)
#            print(io,"Variance for outliers:")
#            writedlm(io,VarOutliers)
#            print(io,"pBernoulli:")
#            writedlm(io,pBernoulli)
#            print(io,"K0:")
#            writedlm(io,K0)
#            print(io,"lambda:")
#            writedlm(io,lambda)
#            print(io,"lambdao:")
#            writedlm(io,lambdao)
#            print(io,"lambdaFMHSDM:")
#            writedlm(io,lambdaFMHSDM)
#            print(io,"JFMHSDM:")
#            writedlm(io,JFMHSDM)
#            print(io,"epsilonORRLS:")
#            writedlm(io,epsilonORRLS)
#            print(io,"lambdaORRLSl1:")
#            writedlm(io,lambdaORRLSl1)
#            print(io,"rhoADMM:")
#            writedlm(io,rhoADMM)
#            print(io,"JADMM:")
#            writedlm(io,JADMM)
#            print(io,"thetaMCP:")
#            writedlm(io,thetaMCP)
#            print(io,"lambdaMCP:")
#            writedlm(io,lambdaMCP)
#            print(io,"JGIST:")
#            writedlm(io,JGIST)
#            print(io,"SD_TimeRLS:")
#            writedlm(io,SD_TimeRLS)

#        end;

#  writedlm("VarLinData_1.txt", Data)

#  open("VarLinTime_1.txt", "w") do inou
#     print("###################################################\r\n")
#     print(inou, "Average Time\r\n")
#     print(inou,"####################################################\r\n")
#     print(inou, "Average_TimeRLS:")
#     writedlm(inou,TimeRLSav)
#     print(inou,"Average_TimeORRLSl1:")
#     writedlm(inou,TimeORRLSl1av)
#     print(inou,"Average_TimeORRLSMCP:")
#     writedlm(inou,TimeORRLSMCPav)
#     print(inou,"Average_TimeORRLSCDl1:")
#     writedlm(inou,TimeORRLSCDl1av)
#     print(inou,"Average_TimeORRLSCDMCP:")
#     writedlm(inou,TimeORRLSCDMCPav)
#     print(inou,"Average_TimeORHRLSADMM:")
#     writedlm(inou,TimeORHRLSADMMav)
#     print(inou,"Average_TimeORHRLSGIST:")
#     writedlm(inou,TimeORHRLSGISTav)
#     print(inou,"Average_TimeORHRLSFMHSDM:")
#     writedlm(inou,TimeORHRLSFMHSDMav)
#     print("###################################################\r\n")
#     print(inou, "Standard Deviation\r\n")
#     print(inou,"####################################################\r\n")
#     print(inou, "SD_TimeRLS:")
#     writedlm(inou,SD_TimeRLS)
#     print(inou, "SD_TimeORRLSl1:")
#     writedlm(inou,SD_TimeORRLSl1)
#     print(inou, "SD_TimeORRLSMCP:")
#     writedlm(inou,SD_TimeORRLSMCP)
#     print(inou, "SD_TimeORRLSCDl1:")
#     writedlm(inou,SD_TimeORRLSCDl1)
#     print(inou, "SD_TimeORRLSCDMCP:")
#     writedlm(inou,SD_TimeORRLSCDMCP)
#     print(inou, "SD_TimeORHRLSADMM:")
#     writedlm(inou,SD_TimeORHRLSADMM)
#     print(inou, "SD_TimeORHRLSGIST:")
#     writedlm(inou,SD_TimeORHRLSGIST)
#     print(inou, "SD_TimeORHRLSFMHSDM:")
#     writedlm(inou,SD_TimeORHRLSFMHSDM)
# end;

max_stars_repo_name: KostasSlavakis/OR-HO-RLS

max_issues_repo_name: KostasSlavakis/OR-HO-RLS

max_forks_repo_name: KostasSlavakis/OR-HO-RLS


text:
# This is the getindex function for Assoc, the Associate Array.
import Base.getindex
StringOrNumArray  = Union{AbstractString,Array,Number}

#The Base getindex function which most higher level function would call upon.
function getindex(A::Assoc, i::Array{Int64}, j::Array{Int64})
    #Check if A is empty
    if nnz(A.A) == 0
    return Assoc([1],[1],0,(+))
    end

    return deepCondense(Assoc(A.row[i],A.col[j],A.val,A.A[i,j]))
    end

#Singular Case
getindex(A::Assoc,i::Any)                  = getindex(A,i,:)

PreviousTypes = Array{Int64}

#Variations by derivations
#Each of the additional type would call upon the base functions above.
#It does so by building up row and column combinations of the new type with all previous types.

#Get index with basic addressing.
#Variations between Element, Array, Colon, Range
getindex(A::Assoc,i::Array{Union{AbstractString,Number}},j::PreviousTypes)                                       = getindex(A,find(x-> x in i,A.row),j)
getindex(A::Assoc,i::PreviousTypes,j::Array{Union{AbstractString,Number}})                                       = getindex(A,i,find(x-> x in j,A.col))
getindex(A::Assoc,i::Array{Union{AbstractString,Number}},j::Array{Union{AbstractString,Number}})                 = getindex(A,find(x-> x in i,A.row),find(x-> x in j,A.col))

PreviousTypes = Union{PreviousTypes,Array{Union{AbstractString,Number}}}

getindex(A::Assoc,i::Int64,j::PreviousTypes)         = getindex(A,[i],j)
getindex(A::Assoc,i::PreviousTypes,j::Int64)         = getindex(A,i,[j])
getindex(A::Assoc,i::Int64,j::Int64)                 = getindex(A,[i],[j])

PreviousTypes = Union{PreviousTypes,Int64}

getindex(A::Assoc,i::Colon,j::PreviousTypes)         = getindex(A,1:size(A.row,1),j)
getindex(A::Assoc,i::PreviousTypes,j::Colon)         = getindex(A,i,1:size(A.col,1))
getindex(A::Assoc,i::Colon,j::Colon)                 = getindex(A,1:size(A.row,1),1:size(A.col,1))

PreviousTypes = Union{PreviousTypes,Colon}

getindex(A::Assoc,i::Range,j::PreviousTypes)         = getindex(A,collect(i),j)
getindex(A::Assoc,i::PreviousTypes,j::Range)         = getindex(A,i,collect(j))
getindex(A::Assoc,i::Range,j::Range)                 = getindex(A,collect(i),collect(j))

PreviousTypes = Union{PreviousTypes,Range}

#Variations of Single Sequence Strings that is separated by a single character separator.
getindex(A::Assoc, i::AbstractString, j::PreviousTypes)   = getindex(A, find( x -> in(x,StrUnique(i)[1]),A.row), j)
getindex(A::Assoc, i::PreviousTypes ,j::AbstractString)   = getindex(A, i ,find( x -> in(x,StrUnique(j)[1]),A.col))
getindex(A::Assoc, i::AbstractString ,j::AbstractString)  = getindex(A, find( x -> in(x,StrUnique(i)[1]),A.row) ,find( x -> in(x,StrUnique(j)[1]),A.col))

PreviousTypes = Union{PreviousTypes,AbstractString}

#Variations by Regex
getindex(A::Assoc, i::Regex, j::PreviousTypes)  = getindex(A, find( x -> ismatch(i,x),A.row), j)
getindex(A::Assoc, i::PreviousTypes, j::Regex)  = getindex(A, i, find( x -> ismatch(j,x),A.col))
getindex(A::Assoc, i::Regex, j::Regex)          = getindex(A, find( x -> ismatch(i,x),A.row), find( x -> ismatch(j,x),A.col))

PreviousTypes = Union{PreviousTypes,Regex}

#Variation with StartsWith
type StartsWith
    inputString::AbstractString
end

function StartsWithHelper(Ar::Array{Union{AbstractString,Number}},S::StartsWith)
    str_list = []
    if S.inputString[end] == ','
        str_list,~ = StrUnique(S.inputString)
    else
        str_list = [S.inputString]
    end
    result_indice = Array{Int64,1}()
    for str in str_list
        str_result_first = searchsortedfirst(Ar,str)
        str_result_last = searchsortedlast(Ar,str*string(Char(255)))
        if str_result_first <= str_result_last
            [push!(result_indice,x) for x in str_result_first:str_result_last]
        end
    end
    return result_indice
end

getindex(A::Assoc,i::PreviousTypes,j::StartsWith) = getindex(A,i,StartsWithHelper(Col(A),j))

getindex(A::Assoc,i::StartsWith,j::PreviousTypes) = getindex(A,StartsWithHelper(Row(A),i),j)

getindex(A::Assoc,i::StartsWith,j::StartsWith) = getindex(A,StartsWithHelper(Row(A),i),StartsWithHelper(Col(A),j))

PreviousTypes = Union{PreviousTypes,StartsWith}

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
# ## Basis Set Data
# Note use of curly braces here. Julia assumes that if you have square braces, you want
# things flattened as much as possible (to be as fast as possible, I guess). Curlys
# preserve the list structure the way I would expect from Python

const sto3g = [
    # H
    [('S',
      [(3.4252509099999999, 0.15432897000000001),
       (0.62391373000000006, 0.53532813999999995),
       (0.16885539999999999, 0.44463454000000002)])],
    # He
    [('S',
      [(6.3624213899999997, 0.15432897000000001),
       (1.1589229999999999, 0.53532813999999995),
       (0.31364978999999998, 0.44463454000000002)])],
    # Li
    [('S',
      [(16.119575000000001, 0.15432897000000001),
       (2.9362007000000001, 0.53532813999999995),
       (0.79465050000000004, 0.44463454000000002)]),
     ('S',
      [(0.63628969999999996, -0.099967230000000004),
       (0.14786009999999999, 0.39951282999999999),
       (0.048088699999999998, 0.70011546999999996)]),
     ('P',
      [(0.63628969999999996, 0.15591627),
       (0.14786009999999999, 0.60768372000000004),
       (0.048088699999999998, 0.39195739000000002)])],
    # Be
    [('S',
      [(30.167871000000002, 0.15432897000000001),
       (5.4951153000000001, 0.53532813999999995),
       (1.4871927, 0.44463454000000002)]),
     ('S',
      [(1.3148331, -0.099967230000000004),
       (0.3055389, 0.39951282999999999),
       (0.099370700000000006, 0.70011546999999996)]),
     ('P',
      [(1.3148331, 0.15591627),
       (0.3055389, 0.60768372000000004),
       (0.099370700000000006, 0.39195739000000002)])],
    # B
    [('S',
      [(48.791113000000003, 0.15432897000000001),
       (8.8873622000000001, 0.53532813999999995),
       (2.4052669999999998, 0.44463454000000002)]),
     ('S',
      [(2.2369561, -0.099967230000000004),
       (0.51982050000000002, 0.39951282999999999),
       (0.16906180000000001, 0.70011546999999996)]),
     ('P',
      [(2.2369561, 0.15591627),
       (0.51982050000000002, 0.60768372000000004),
       (0.16906180000000001, 0.39195739000000002)])],
    # C
    [('S',
      [(71.616837000000004, 0.15432897000000001),
       (13.045095999999999, 0.53532813999999995),
       (3.5305122, 0.44463454000000002)]),
     ('S',
      [(2.9412493999999998, -0.099967230000000004),
       (0.68348310000000001, 0.39951282999999999),
       (0.22228990000000001, 0.70011546999999996)]),
     ('P',
      [(2.9412493999999998, 0.15591627),
       (0.68348310000000001, 0.60768372000000004),
       (0.22228990000000001, 0.39195739000000002)])],
    # N
    [('S',
      [(99.106168999999994, 0.15432897000000001),
       (18.052312000000001, 0.53532813999999995),
       (4.8856602000000002, 0.44463454000000002)]),
     ('S',
      [(3.7804559000000002, -0.099967230000000004),
       (0.87849659999999996, 0.39951282999999999),
       (0.28571439999999998, 0.70011546999999996)]),
     ('P',
      [(3.7804559000000002, 0.15591627),
       (0.87849659999999996, 0.60768372000000004),
       (0.28571439999999998, 0.39195739000000002)])],
    # O
    [('S',
      [(130.70931999999999, 0.15432897000000001),
       (23.808861, 0.53532813999999995),
       (6.4436083000000002, 0.44463454000000002)]),
     ('S',
      [(5.0331513000000001, -0.099967230000000004),
       (1.1695960999999999, 0.39951282999999999),
       (0.38038899999999998, 0.70011546999999996)]),
     ('P',
      [(5.0331513000000001, 0.15591627),
       (1.1695960999999999, 0.60768372000000004),
       (0.38038899999999998, 0.39195739000000002)])],
    # F
    [('S',
      [(166.67912999999999, 0.15432897000000001),
       (30.360811999999999, 0.53532813999999995),
       (8.2168206999999995, 0.44463454000000002)]),
     ('S',
      [(6.4648032000000004, -0.099967230000000004),
       (1.5022812000000001, 0.39951282999999999),
       (0.48858849999999998, 0.70011546999999996)]),
     ('P',
      [(6.4648032000000004, 0.15591627),
       (1.5022812000000001, 0.60768372000000004),
       (0.48858849999999998, 0.39195739000000002)])],
    # Ne
    [('S',
       [(207.01561000000001, 0.15432897000000001),
        (37.708151000000001, 0.53532813999999995),
        (10.205297, 0.44463454000000002)]),
      ('S',
       [(8.2463151000000003, -0.099967230000000004),
        (1.9162661999999999, 0.39951282999999999),
        (0.62322929999999999, 0.70011546999999996)]),
      ('P',
       [(8.2463151000000003, 0.15591627),
        (1.9162661999999999, 0.60768372000000004),
            (0.62322929999999999, 0.39195739000000002)])]
]
basis_set_data = Dict("sto3g" => sto3g)

mutable struct BasisSet # list of CGBFs
    bfs::Array{CGBF,1}
end

basisset() = BasisSet(CGBF[])

function push!(basis::BasisSet,cbf::CGBF)
    Base.push!(basis.bfs,cbf)
end

function build_basis(mol::Molecule,name="sto3g")
    data = basis_set_data[name]
    basis_set = basisset()
    for atom in mol.atomlist
        for btuple in data[atom.atno]
            sym,primlist = btuple
            for (I,J,K) in sym2power[sym]
                cbf = cgbf(atom.x,atom.y,atom.z,I,J,K)
                push!(basis_set,cbf)
                for (expn,coef) in primlist
                    push!(cbf,expn,coef)
                end
            end
        end
    end
    return basis_set
end

const sym2power = Dict(
    'S' => [(0,0,0)],
    'P' => [(1,0,0),(0,1,0),(0,0,1)],
    'D' => [(2,0,0),(0,2,0),(0,0,2),(1,1,0),(1,0,1),(0,1,1)]
    )

function test_geo_basis()
    @testset "test_geo_basis" begin
    @test nuclear_repulsion(h2) ≈ 0.7223600367
    @test nel(h2) == 2
    @test nel(h2o) == 10
    @test length(sto3g)==10
    bfs = build_basis(h2)
    @test length(bfs.bfs)==2
    l,r = bfs.bfs
    @test overlap(l,l)      ≈ 1
    @test overlap(r,r)      ≈ 1
    @test overlap(l,r)      ≈ 0.6647387449282997
    @test kinetic(l,l)      ≈ 0.76003188
    @test kinetic(r,r)      ≈ 0.76003188
    @test kinetic(l,r)      ≈ 0.24141861181119084
    @test coulomb(l,l,l,l)  ≈ 0.7746059439196398
    @test coulomb(r,r,r,r)  ≈ 0.7746059439196398
    @test coulomb(l,l,r,r)  ≈ 0.5727937653511646
    @test coulomb(l,l,l,r)  ≈ 0.4488373301593464
    @test coulomb(l,r,l,r)  ≈ 0.3025451156654606
    bfs = build_basis(h2o)

    s1,s2,px,py,pz,hl,hr = bfs.bfs
    @test coulomb(s1,s2,hl,hr) ≈ 0.03855344493645537
    @test coulomb(s1,pz,hl,hr) ≈ -0.0027720110485359053
    @test coulomb(s1,hl,pz,hr) ≈ -0.010049491284827426
    @test coulomb(s1,py,hl,hr)==0
    @test coulomb(s1,hl,py,hr)==0
    end #testset
end


text:
function createabcfunction(prob, t, distancefunction, alg; save_idxs = nothing, sample_u0 = false, kwargs...)
    function simfunc(params, constants, data)
        local u0
        if sample_u0
            u0 = save_idxs === nothing ? params[1:length(prob.u0)] : params[1:length(save_idxs)]
            if length(u0) < length(prob.u0)
                for i in length(u0):length(prob.u0)
                    push!(u0,prob.u0[i])
                end
            end
        else
            u0 = prob.u0
        end
        sol = solve(prob, alg, u0=u0, p=params, saveat = t, save_idxs = save_idxs, kwargs...)
        if size(sol, 2) < length(t)
            return Inf,nothing
        else
            simdata = convert(Array, sol)
            return distancefunction(data, simdata), nothing
        end
    end
end

function abc_inference(prob::DiffEqBase.DEProblem, alg, t, data, priors; ϵ=0.001,
                       distancefunction = euclidean, ABCalgorithm = ABCSMC, progress = false,
                       num_samples = 500, maxiterations = 10^5, save_idxs = nothing, sample_u0 = false, parallel = false, kwargs...)

    abcsetup = ABCalgorithm(createabcfunction(prob, t, distancefunction, alg; save_idxs = save_idxs, sample_u0 = sample_u0, kwargs...),
                            length(priors),
                            ϵ,
                            ApproxBayes.Prior(priors);
                            nparticles = num_samples,
                            maxiterations = maxiterations
                            )

    abcresult = runabc(abcsetup, data, progress = progress, parallel = parallel)
    return abcresult
end


text:
################################################################################
#
#          EllCrv/EllCrv.jl : Elliptic curves over general fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2016 Tommy Hofmann
# (C) 2016 Robin Ammon
# (C) 2016 Sofia Brenner
# (C) 2022 Jeroen Hanselman
#
################################################################################

export EllCrv, EllCrvPt

export base_field, division_polynomial, division_polynomial_univariate, EllipticCurve, infinity,
       isfinite, isinfinite, isweierstrassmodel, is_on_curve, j_invariant,
       short_weierstrass_model, +, *, a_invars, b_invars, c_invars, equation

################################################################################
#
#  Types
#
################################################################################

mutable struct EllCrv{T}
  base_field::Ring
  short::Bool
  a_invars::Tuple{T, T, T, T, T}
  b_invars::Tuple{T, T, T, T}
  c_invars::Tuple{T,T}
  disc::T
  j::T
  coeff::Vector{T} # remove usage in Oscar and then remove this here

  torsion_points#::Vector{EllCrvPt}
  torsion_structure#Tuple{Vector{Int}, Vector{EllCrvPt}}

  function EllCrv{T}(coeffs::Vector{T}, check::Bool = true) where {T}
    if length(coeffs) == 2
      if check
        d = -16*(4*coeffs[1]^3 + 27*coeffs[2]^2)
        if d != 0
          E = new{T}()
          E.short = true
          # fixed on Nemo master
          K = parent(coeffs[1])
          E.base_field = K
          E.a_invars = (zero(K),zero(K),zero(K),coeffs[1],coeffs[2])

          E.disc = d
        else
          error("discriminant is zero")
        end
      else
        E = new{T}()
        E.short = true
        K = parent(coeffs[1])
        E.base_field = K
        E.a_invars = (zero(K),zero(K),zero(K),coeffs[1],coeffs[2])
      end
    elseif length(coeffs) == 5 # coeffs = [a1, a2, a3, a4, a6]
      if check
        a1 = coeffs[1]
        a2 = coeffs[2]
        a3 = coeffs[3]
        a4 = coeffs[4]
        a6 = coeffs[5]

        b2 = a1^2 + 4*a2
        b4 = a1*a3 + 2*a4
        b6 = a3^2 + 4*a6
        b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

        d = (-b2^2*b8 - 8*b4^3 - 27*b6^2 + 9*b2*b4*b6)

        if d != 0
          E = new{T}()
          if a1==a2==a3==0
            E.short = true
          else
            E.short = false
          end
          E.a_invars = (a1, a2, a3, a4, a6)
          E.b_invars = (b2, b4, b6, b8)
          E.disc = d
          E.base_field = parent(coeffs[1])
        else
          error("Discriminant is zero")
        end
      else
        E = new{T}()
        a1 = coeffs[1]
        a2 = coeffs[2]
        a3 = coeffs[3]
        a4 = coeffs[4]
        a6 = coeffs[5]
        if a1==a2==a3==0
            E.short = true
          else
            E.short = false
          end
        E.a_invars = (a1, a2, a3, a4, a6)
        E.base_field = parent(coeffs[1])
      end
    else
      error("Length of coefficient array must be either 2 or 5")
    end
    E.coeff = coeffs
    return E
  end
end

mutable struct EllCrvPt{T}
  coordx::T
  coordy::T
  isinfinite::Bool
  parent::EllCrv{T}

  function EllCrvPt{T}(E::EllCrv{T}, coords::Vector{T}, check::Bool = true) where {T}
    if check
      if is_on_curve(E, coords)
        P = new{T}(coords[1], coords[2], false, E)
        return P
      else
        error("Point is not on the curve")
      end
    else
      P = new{T}(coords[1], coords[2], false, E)
      return P
    end
  end

  function EllCrvPt{T}(E::EllCrv{T}) where {T}
    z = new{T}()
    z.parent = E
    z.isinfinite = true
    return z
  end
end

################################################################################
#
#  Constructors for Elliptic Curve
#
################################################################################

function EllipticCurve(x::Vector{T}, check::Bool = true) where T <: RingElem
  E = EllCrv{T}(x, check)
  return E
end

function EllipticCurve(K::Field, x::Vector{T}, check::Bool = true) where T
  EllipticCurve([ K(z) for z in x], check)
end

#  Implicit promotion in characterstic 0
function EllipticCurve(x::Vector{Int}, check::Bool = true)
  return EllipticCurve(fmpq[ FlintQQ(z) for z in x], check)
end

function EllipticCurve(x::Vector{BigInt}, check::Bool = true)
  return EllipticCurve(fmpq[ FlintQQ(z) for z in x], check)
end

function EllipticCurve(x::Vector{fmpz}, check::Bool = true)
  return EllipticCurve(fmpq[ FlintQQ(z) for z in x], check)
end

function EllipticCurve(x::Vector{Rational{Int}}, check::Bool = true)
  return EllipticCurve(fmpq[ FlintQQ(z) for z in x], check)
end

################################################################################
#
#  Constructors for Point on Elliptic Curve
#
################################################################################

function (E::EllCrv{T})(coords::Vector{S}, check::Bool = true) where {S, T}
  if length(coords) != 2
    error("Need two coordinates")
  end

  if S == T
    parent(coords[1]) != base_field(E) &&
        error("Objects must be defined over same field")
    return EllCrvPt{T}(E, coords, check)
  else
    return EllCrvPt{T}(E, map(base_field(E), coords), check)
  end
end

################################################################################
#
#  Equality of Models
#
################################################################################

@doc Markdown.doc"""
    ==(E::EllCrv, F::EllCrv) -> Bool

Return true if $E$ and $F$ are given by the same model over the same field.
"""
function ==(E::EllCrv, F::EllCrv) where T
  return a_invars(E) == a_invars(F) && base_field(E) == base_field(F)
end

################################################################################
#
#  Field access
#
################################################################################

@doc Markdown.doc"""
    base_field(E::EllCrv) -> Field

Return the base field over which E is defined.
"""
function base_field(E::EllCrv{T}) where T
  return E.base_field::parent_type(T)
end

function Base.deepcopy_internal(E::EllCrv, dict::IdDict)
    return EllipticCurve(E.a_invars)
end

function parent(P::EllCrvPt)
  return P.parent
end

@doc Markdown.doc"""
    isfinite(E::EllCrvPt) -> Bool

Return true if P is not the point at infinity.
"""
function isfinite(P::EllCrvPt)
  return !P.isinfinite
end

@doc Markdown.doc"""
    isinfinite(E::EllCrvPt) -> Bool

Return true if P is the point at infinity.
"""
function isinfinite(P::EllCrvPt)
  return P.isinfinite
end

@doc Markdown.doc"""
    isweierstrassmodel(E::EllCrv) -> Bool

Return true if E is in short Weierstrass form.
"""
function isweierstrassmodel(E::EllCrv)
  return E.short
end

@doc Markdown.doc"""
    a_invars(E::EllCrv{T}) -> Tuple{T, T, T, T, T}

Return the Weierstrass coefficients of E as a tuple (a1, a2, a3, a4, a6)
such that E is given by y^2 + a1xy + a3y = x^3 + a2x^2 + a4x + a6.
"""
function a_invars(E::EllCrv)
  return E.a_invars
end

@doc Markdown.doc"""
    coefficients(E::EllCrv{T}) -> Tuple{T, T, T, T, T}

Return the Weierstrass coefficients of E as a tuple (a1, a2, a3, a4, a6)
such that E is given by y^2 + a1xy + a3y = x^3 + a2x^2 + a4x + a6.
"""
coefficients(E::EllCrv) = a_invars(E)

@doc Markdown.doc"""
    b_invars(E::EllCrv{T}) -> Tuple{T, T, T, T}

Return the b-invariants of E as a tuple (b2, b4, b6, b8).
"""
function b_invars(E::EllCrv)
  if isdefined(E, :b_invars)
    # fixed on Nemo master
    return deepcopy(E.b_invars)
  else
    a1, a2, a3, a4, a6 = a_invars(E)

    b2 = a1^2 + 4*a2
    b4 = a1*a3 + 2*a4
    b6 = a3^2 + 4*a6
    b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

    E.b_invars = (b2, b4, b6, b8)
    return (b2, b4, b6, b8)
  end
end

@doc Markdown.doc"""
    c_invars(E::EllCrv{T}) -> Tuple{T, T, T, T}

Return the c-invariants of E as a tuple (c4, c6).
"""
function c_invars(E::EllCrv)
  if isdefined(E, :c_invars)
    # fixed on Nemo master
    return deepcopy(E.c_invars)
  else
    b2,b4,b6,b8 = b_invars(E)
    c4 = b2^2 - 24*b4
    c6 = -b2^3 + 36*b2*b4 - 216*b6
    return (c4,c6)
  end
end

################################################################################
#
# Changing the model
#
################################################################################

@doc Markdown.doc"""
    short_weierstrass_model(E::EllCrv{fmpq}) ->
      (EE::EllCrv, function(EllCrvPt), function(EllCrvPt))

Transform a curve given in long Weierstrass form over QQ to short Weierstrass
form. Return short form and both transformations for points on the curve;
first transformation from E (long form) to EE (short form),
second transformation the other way round
"""
function short_weierstrass_model(E::EllCrv)
  return _short_weierstrass_model(E)
end

function _short_weierstrass_model(E::EllCrv{T}) where T
  R = base_field(E)
  p = characteristic(R)

  if (p == 2) || (p == 3)
      error("Converting to short form not possible in characteristic 2 and 3")
  end

  a1, _, a3= a_invars(E)

  b2, b4, b6, b8 = b_invars(E)

  c4, c6 = c_invars(E)

  Anew = -divexact(c4, 48)
  Bnew = -divexact(c6, 864)

  EE = EllipticCurve([Anew, Bnew])::EllCrv{T}

  # we are hitting https://github.com/JuliaLang/julia/issues/15276

  _b2 = deepcopy(b2)
  _a1 = deepcopy(a1)
  _a3 = deepcopy(a3)

  # transforms a point on E (long form) to a point on EE (short form)
  trafo = function(P::EllCrvPt)

    if P.isinfinite
      return infinity(EE)
    end

    xnew = P.coordx + divexact(_b2, 12)
    ynew = P.coordy + divexact(_a1*P.coordx + _a3, 2)
    Q = EE([xnew, ynew])::EllCrvPt{T}
    return Q
  end

  # transforms a point on EE (short form) back to a point on E (long form)
  ruecktrafo = function(R::EllCrvPt)
    if R.isinfinite
        return infinity(E)
    end

    xnew = R.coordx - divexact(_b2, 12)
    ynew = R.coordy - divexact(_a1*xnew + _a3, 2)
    S = E([xnew, ynew])
    return S::EllCrvPt{T}
  end

  # type annotation necessary due to #15276
  return EE::EllCrv{T}, trafo, ruecktrafo
end

################################################################################
#
#  Equation
#
################################################################################

@doc Markdown.doc"""
    equation(E::EllCrv) -> Poly

Return the equation defining the elliptic curve E.
"""
function equation(E::EllCrv)
  K = base_field(E)
  Kx, x = PolynomialRing(K,"x")
  Kxy, y = PolynomialRing(Kx,"y")

  a1, a2, a3, a4, a6 = a_invars(E)

  return y^2 + a1*x*y + a3*y - (x^3 + a2*x^2 + a4*x + a6)
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, E::EllCrv)
  print(io, "Elliptic curve with equation\n")
  a1, a2, a3, a4, a6 = a_invars(E)
  sum = Expr(:call, :+)
  push!(sum.args, Expr(:call, :^, :y, 2))
  c = a1
  if !iszero(c)
    if isone(c)
      push!(sum.args, Expr(:call, :*, :x, :y))
    else
      push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), :x, :y))
    end
  end
  c = a3
  if !iszero(c)
    if isone(c)
      push!(sum.args, Expr(:call, :*, :y))
    else
      push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), :y))
    end
  end
  print(io, AbstractAlgebra.expr_to_string(AbstractAlgebra.canonicalize(sum)))

  print(io, " = ")
  sum = Expr(:call, :+)
  push!(sum.args, Expr(:call, :^, :x, 3))

  c = a2
  if !iszero(c)
    if isone(c)
      push!(sum.args, Expr(:call, :*, Expr(:call, :^, :x, 2)))
    else
      push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), Expr(:call, :^, :x, 2)))
    end
  end

  c = a4
  if !iszero(c)
    if isone(c)
      push!(sum.args, Expr(:call, :*, :x))
    else
      push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), :x))
    end
  end

  c = a6
  if !iszero(c)
    push!(sum.args, AbstractAlgebra.expressify(c, context = io))
  end

  print(io, AbstractAlgebra.expr_to_string(AbstractAlgebra.canonicalize(sum)))
end

function show(io::IO, P::EllCrvPt)
    if P.isinfinite
        print(io, "Point at infinity of $(P.parent)")
    else
        print(io, "Point $(P.coordx),$(P.coordy) of $(P.parent)")
    end
end

################################################################################
#
#  Point at infinity
#
################################################################################

@doc Markdown.doc"""
    infinity(E::EllCrv) -> EllCrvPt

Return the point at infinity.
"""
function infinity(E::EllCrv{T}) where T
  infi = EllCrvPt{T}(E)
  return infi
end

################################################################################
#
#  Test for inclusion
#
################################################################################

@doc Markdown.doc"""
    is_on_curve(E::EllCrv{T}, coords::Vector{T}) -> Bool

Return true if `coords` defines a point on $E$ and false otherwise. The array
`coords` must have length 2.
"""
function is_on_curve(E::EllCrv{T}, coords::Vector{T}) where T
  length(coords) != 2 && error("Array must be of length 2")
  a1, a2, a3, a4, a6 = a_invars(E)
  x = coords[1]
  y = coords[2]

  if E.short == true
    if y^2 == x^3 + a4*x + a6
      return true
    else
      return false
    end
  else
    if (y^2 + a1*x*y + a3*y ==
            x^3 + a2*x^2+a4*x + a6)
      return true
    else
      return false
    end
  end
end

################################################################################
#
#  Discriminant
#
################################################################################

@doc Markdown.doc"""
    discriminant(E::EllCrv{T}) -> T

Compute the discriminant of $E$.
"""
function discriminant(E::EllCrv{T}) where T
  if isdefined(E, :disc)
    return E.disc
  end
  if E.short == true
    _, _, _, a4, a6 = a_invars(E)
    d = -16*(4*a4^3 + 27*a6^2)
    E.disc = d
    return d::T
  else
    b2, b4, b6, b8 = b_invars(E)
    d = -b2^2*b8 - 8*b4^3 - 27*b6^2 + 9*b2*b4*b6
    E.disc = d
    return d::T
  end
end

################################################################################
#
#  j-invariant
#
################################################################################

# p. 46 Washington, p. 72 Cohen
@doc Markdown.doc"""
    j(E::EllCrv{T}) -> T

Compute the j-invariant of $E$.
"""
function j_invariant(E::EllCrv{T}) where T
  if isdefined(E, :j)
    return E.j
  end

  if E.short == true

    R = base_field(E)
    a1, a2, a3, a4, a6 = a_invars(E)
    j = divexact(-1728*(4*a4)^3,discriminant(E))
    E.j = j
    return j::T
  else

    c4, c6 = c_invars(E)

    j = divexact(c4^3, discriminant(E))
    E.j = j
    return j::T
  end
end

################################################################################
#
#  Addition of Points
#
################################################################################

# washington p. 14, cohen p. 270
@doc Markdown.doc"""
    +(P::EllCrvPt, Q::EllCrvPt) -> EllCrvPt

Add two points on an elliptic curve.
"""
function +(P::EllCrvPt{T}, Q::EllCrvPt{T}) where T
  parent(P) != parent(Q) && error("Points must live on the same curve")

  # Is P = infinity or Q = infinity?
  if P.isinfinite
      return Q
  elseif Q.isinfinite
      return P
  end

  E = P.parent

  # Distinguish between long and short form
  if E.short == true
    if P.coordx != Q.coordx
        m = divexact(Q.coordy - P.coordy, Q.coordx - P.coordx)
        x = m^2 - P.coordx - Q.coordx
        y = m * (P.coordx - x) - P.coordy
    elseif P.coordy != Q.coordy
        return infinity(E)
    elseif P.coordy != 0
        _, _, _, a4 = a_invars(E)
        m = divexact(3*(P.coordx)^2 + a4, 2*P.coordy)
        x = m^2 - 2*P.coordx
        y = m* (P.coordx - x) - P.coordy
    else
        return infinity(E)
    end

    Erg = E([x, y], false)

  else
  a1, a2, a3, a4, a6 = a_invars(E)

    # Use [Cohen, p. 270]
    if P.coordx == Q.coordx
      if Q.coordy == -a1*P.coordx - a3 - P.coordy # then P = -Q
        return infinity(E)
      elseif P.coordy == Q.coordy # then P = Q
        m = divexact(3*((P.coordx)^2) + 2*a2*P.coordx + a4 - a1*P.coordy, 2*P.coordy + a1*P.coordx + a3)
        x = -P.coordx - Q.coordx - a2 + a1*m + m^2
        y = -P.coordy - m*(x - P.coordx) - a1*x - a3
      else # then P != +-Q
        m = divexact(Q.coordy - P.coordy, Q.coordx - P.coordx)
        x = -P.coordx - Q.coordx - a2 + a1*m + m^2
        y = -P.coordy - m*(x - P.coordx) - a1*x - a3
      end
    else # now P != +-Q
      m = divexact(Q.coordy - P.coordy, Q.coordx - P.coordx)
      x = -P.coordx - Q.coordx - a2 + a1*m + m^2
      y = -P.coordy - m*(x - P.coordx) - a1*x - a3
    end

    Erg = E([x, y], false)

  end
  return Erg
end

################################################################################
#
#  Inverse
#
################################################################################

@doc Markdown.doc"""
    -(P::EllCrvPt) -> EllCrvPt

Compute the inverse of the point $P$ on an elliptic curve.
"""
function -(P::EllCrvPt)
  E = P.parent

  if !isfinite(P)
    return infinity(E)
  end

  if E.short == true
    Q = E([P.coordx, -P.coordy], false)
  else
    a1,_, a3 = a_invars(E)
    Q = E([P.coordx, -a1*P.coordx - a3 - P.coordy], false)
  end

  return Q
end

@doc Markdown.doc"""
    ==(P::EllCrvPt, Q::EllCrvPt) -> Bool

Return true if $P$ and $Q$ are equal and live over the same elliptic curve
$E$.
"""
function ==(P::EllCrvPt{T}, Q::EllCrvPt{T}) where T
  # both are infinite
  if P.isinfinite && Q.isinfinite
    return true
  end

  # one of them is infinite
  if xor(P.isinfinite, Q.isinfinite)
    return false
  end

  # Otherwise, compare coordinates
  if P.coordx == Q.coordx && P.coordy == Q.coordy
    return true
  else
    return false
  end
end

################################################################################
#
#  Scalar multiplication
#
################################################################################

# algorithm 'integer times a point', [Washington, p. 18]
@doc Markdown.doc"""
    *(n::Int, P::EllCrvPt) -> EllCrvPt

Compute the point $nP$.
"""
function *(n::Int, P::EllCrvPt)
  B = infinity(P.parent)
  C = P

  if n >= 0
      a = n
  else
      a = -n
  end

  while a != 0
    if mod(a,2) == 0
      a = div(a,2)
      C = C + C
    else
      a = a - 1
      B = B + C
    end
  end

  if n < 0
    B = -B
  end

  return B
end

################################################################################
#
#  Division polynomials
#
################################################################################
"""
    division_polynomial_univariate(E::EllCrv, n::Int, [x]) -> Poly

Compute the n-th univariate division polynomial of an elliptic curve defined
over a field k following Mazur and Tate. By default the result is a univariate polynomial over the base ring of `E`.
When `x` is given, the output is evaluated using the given value for `x`.

A triple of objects is returned:
- The n-th division polynomial as a univariate polynomial with a mutiplicity
   of 2 at the non-zero two-torsion points.
- The n-th division polynomial as a univariate polynomial divided by the
  univariate 2-torsion polynomial when n is even.
- The complementary factor, i.e. the first output divided by the second output.
"""
function division_polynomial_univariate(E::EllCrv, n::Int, x = PolynomialRing(base_field(E),"x")[2])

  R = parent(x)
  if isweierstrassmodel(E)
    poly = divpol_g_short(E,n,x)
    if mod(n,2) == 0
      _, _, _, A, B = a_invars(E)
      twotorsfactor = 4*(x^3+A*x+B)
    else
      twotorsfactor = one(R)
    end
  else
    poly = divpol_g(E,n,x)
      if mod(n,2) == 0
        b2, b4, b6 = b_invars(E)
        twotorsfactor = 4*x^3+b2*x^2+2*b4*x+b6
      else
        twotorsfactor = one(R)
      end
  end
  return twotorsfactor*poly, poly, twotorsfactor
end

"""
    division_polynomial(E::EllCrv, n::Int, x, y) -> Poly

Compute the n-th division polynomial of an elliptic curve defined over a field
k following Mazur and Tate. When x and or y are given the output is
automatically evaluated using the given values.
"""
function division_polynomial(E::EllCrv, n::Int,x = PolynomialRing(base_field(E),"x")[2], y = PolynomialRing(parent(x),"y")[2])
  R = parent(y)
  if isweierstrassmodel(E)
     if mod(n,2) == 0
      return 2*y*divpol_g_short(E,n,x)
    else
      return R(divpol_g_short(E,n,x))
    end
  else
    a1, _, a3 = a_invars(E)
    if mod(n,2) == 0
      return (2*y + a1*x + a3)*divpol_g(E,n,x)
    else
      return R(divpol_g(E,n,x))
    end
  end
end

function divpol_g_short(E, n::Int, x = PolynomialRing(base_field(E),"x")[2])

  Kx = parent(x)
  _, _, _, A, B = a_invars(E)

  B6sqr = (4*x^3+4*A*x+4*B)^2

  if n == 1 || n == 2
    return one(parent(x))
  elseif n == 3
    return 3*x^4 + 6*(A)*x^2 + 12*(B)*x - (A)^2
  elseif n == 4
    return 2*(x^6 + 5*(A)*x^4 + 20*(B)*x^3 - 5*(A)^2*x^2 - 4*(A)*(B)*x - 8*(B)^2 - (A)^3)
  elseif mod(n,2) == 0
    m = div(n,2)
    return (divpol_g_short(E,m,x))*(divpol_g_short(E,m+2, x)*divpol_g_short(E,m-1,x)^2 - divpol_g_short(E,m-2,x)*divpol_g_short(E,m+1,x)^2)
  else m = div(n-1,2)
    m = div(n-1,2)
    part1 = divpol_g_short(E,m+2,x)  * divpol_g_short(E,m,x)^3
    part2 = divpol_g_short(E,m-1,x) * divpol_g_short(E,m+1,x)^3
    if mod(m,2) == 0
      return B6sqr * part1 - part2
    else
      return part1 - B6sqr * part2
    end
  end
end

function divpol_g(E, n::Int, x = PolynomialRing(base_field(E),"x")[2])

  Kx = parent(x)

  b2, b4, b6, b8 = E.b_invars
  B4 = 6*x^2+b2*x+b4
  B6sqr = (4*x^3+b2*x^2+2*b4*x+b6)^2
  B8 = 3*x^4 + b2*x^3 + 3*b4*x^2 + 3*b6*x + b8

  if n == 1 || n ==2
    return one(Kx)
  elseif n == 3
    return B8
  elseif n == 4
    return -B6sqr+B4*B8
  elseif mod(n,2) == 0
    m = div(n-2,2)
    return divpol_g(E,m+3,x)*divpol_g(E,m+1,x)*divpol_g(E,m,x)^2 - divpol_g(E,m+1,x)*divpol_g(E,m-1,x)*divpol_g(E,m+2,x)^2
  else
    m = div(n-1,2)
    part1 = divpol_g(E,m+2,x)  * divpol_g(E,m,x)^3
    part2 = divpol_g(E,m-1,x) * divpol_g(E,m+1,x)^3
    if mod(m,2) == 0
      return B6sqr * part1 - part2
    else
      return part1 - B6sqr * part2
    end
  end
end

function replace_all_squares_modulo(f, g, F)
    # assumes that f is in Z[x,y^2] and g in Z[x]. Replaces y^2 with g.
    # the result will be in Z[x]
    z = zero(parent(g)) # this is the zero in Z[x]
    d = div(degree(f), 2) # degree of f in y^2 should be even
    for i in 0:d
        z = z + coeff(f, 2*i)*powermod(g, i, F)
    end
    return z
end


text:
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
using Dates

using PowerSystems
using PowerSimulations

using Cbc #solver
solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

PowerSystems.download(PowerSystems.TestData; branch = "master") # *note* add `force=true` to get a fresh copy
base_dir = dirname(dirname(pathof(PowerSystems)));

TAMU_DIR = joinpath(base_dir, "data", "ACTIVSg2000");
sys = TamuSystem(TAMU_DIR)
horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval);

sim_folder = mkpath(joinpath(pkgpath, "TAMU-sim"))
stages_definition = Dict(
    "UC" => Stage(
        GenericOpProblem,
        template_unit_commitment(network = CopperPlatePowerModel),
        sys,
        solver,
    ),
)
order = Dict(1 => "UC")
horizons = Dict("UC" => 24)
intervals = Dict("UC" => (Hour(24), Consecutive()))
DA_sequence = SimulationSequence(
    step_resolution = Hour(24),
    order = order,
    horizons = horizons,
    intervals = intervals,
    ini_cond_chronology = IntraStageChronology(),
)

sim = Simulation(
    name = "TAMU-test",
    steps = 3,
    stages = stages_definition,
    stages_sequence = DA_sequence,
    simulation_folder = sim_folder,
)

build!(sim)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using BenchmarkTools

exps = create_exponents(6, 10)

@btime rand($exps, 232)

binomial(4+3, 3)

using Homotopies

import DynamicPolynomials: @polyvar
import FixedPolynomials
const FP = FixedPolynomials

@polyvar x y
f = x^3*y-x*y
g = x+3y*x^2*y+y^3

H = homogenize(GeodesicOnTheSphere{Complex128}([f, g], [2*g, f]))
H2 = homogenize(StraightLineHomotopy{Complex128}([f, g], [2*g, f]))

JH = jacobian!(H)

u = zeros(Complex128, 2)
w = rand(Complex128, 3)
t = rand()
@benchmark evaluate!($u, $H, $w, $t)
@benchmark evaluate!($u, $H2, $w, $t)

v = zeros(Complex128, 2, 2)

@benchmark JH($v, $w, $t)

x = rand(Complex128, 2)

JH = jacobian!(H)

u = zeros(Complex128, 2, 2)

@benchmark JH($u, $x, $one(Complex128))

J_start = [FP.differentiate(f, i) for f in H.start, i=1:FP.nvariables.(H.start[1])]
p = J_start[1, 1]
y = rand(2)

eval2(f) = FP.evaluate(f, x)

@benchmark map!($eval2, $u, $J_start)

@benchmark FP.evaluate($p, $y)


text:
module mna1  #jl

using ReachabilityAnalysis, MAT, SparseArrays
using ReachabilityModels: @relpath

file = matopen(@relpath "mna1.mat")

# system matrix
A = sparse(read(file, "A"))

# affine term
b = sparsevec(570:578, [fill(-0.1, 5); fill(-0.2, 4)], 578)

function model(X0)
    S = @system(x' = Ax + b)
    return IVP(S, X0)
end

end  #jl

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
module AbelSand

output = include("Output.jl") # supports output functionality (e.g. array to csv)

const N_size = 700;     # the radius of the lattice Z^2, the actual size becomes (2*N+1)x(2*N+1)
const dx = [1,0,-1,0];  # for a given (x,y) in Z^2, (x + dx, y + dy) for all (dx,dy) covers the neighborhood of (x,y)
const dy = [0,1,0,-1];

struct L_coord
  # represents a lattice coordinate
  x::Int
  y::Int
end

function FindCoordinate(Z::Array{L_coord,1}, a::Int, b::Int)
   # in the given array Z of coordinates finds the (first) index of the tuple (a,b);
   # if no match, returns -1

  for i=1:length(Z)
    if (Z[i].x==a)&&(Z[i].y==b)
      return i;
    end
  end

  return -1;
end

function move(N)
    # the main function moving the pile sand grains of size N at the origin of Z^2 until the sandpile becomes stable

    Z_lat = zeros(UInt8, 2*N_size+1,2*N_size+1);  # models the integer lattice Z^2, we will have at most 4 sands on each vertex
    V_sites = falses(2*N_size+1,2*N_size+1);  # all sites which are visited by the sandpile process, are being marked here
    Odometer = zeros(UInt64, 2*N_size+1, 2*N_size+1); # stores the values of the odometer function

    walking = L_coord[];    # the coordinates of sites which need to move

    V_sites[N_size + 1,N_size + 1] = true;

    # i1, ... j2  -> show the boundaries of the box which is visited by the sandpile process
    i1, i2, j1, j2 = N_size + 1, N_size + 1, N_size + 1, N_size + 1 ;
    n = N

    t1 = time_ns();

    while n > 0
        n -= 1;

        Z_lat[N_size + 1, N_size + 1] += 1;
        if ( Z_lat[N_size + 1, N_size + 1] >= 4 )
            push!(walking, L_coord(N_size + 1, N_size + 1));
        end

        while ( length(walking) > 0 )
            w = pop!(walking);
            x = w.x;
            y = w.y;

            Z_lat[x, y] -= 4;
            Odometer[x,y] += 4;

            for k = 1:4
               Z_lat[x + dx[k], y + dy[k]] += 1; V_sites[x + dx[k], y + dy[k]] = true;
               if Z_lat[x + dx[k], y + dy[k]] >= 4
                  if FindCoordinate(walking, x + dx[k] , y + dy[k]) == -1
                      push!(walking, L_coord( x + dx[k], y + dy[k] ) );
                  end
               end
            end

            i1 = min(i1, x - 1);
            i2 = max(i2, x + 1);
            j1 = min(j1, y - 1);
            j2 = max(j2, y + 1);
        end

        #if n%200 == 0
        #    print(string("The boundaries are:: ", (i2-i1+1),"x",(j2-j1+1)), "\n");
        #    #flush(stdout);
        #end

    end #end of the main while
    t2 = time_ns();

    println("The final boundaries are:: ", (i2-i1+1),"x",(j2-j1+1), "\n");
    print( "time elapsed: " , (t2 - t1)/1.0e9, "\n" );

    Z_lat = output.printIntoFile(Z_lat,  0, string("Abel_Z_", N) );
    Odometer = output.printIntoFile(Odometer, 1,string("Abel_OD_", N)  );

    output.saveAsGrayImage(Z_lat, string("Abel_Z_", N), 20, 0);
    color_code = Dict(1=>[255, 128, 255], 2=>[255, 0, 0],3=>[0,128,255])
    output.saveAsRGBImage(Z_lat, string("Abel_Z_color_", N), color_code, 20, 0);

    # for the total elapsed time, it's better to use the @time macros on the main call

    return Z_lat, Odometer; # these are trimmed in output module

end # end of function move

end # end of the module


text:
# Autogenerated wrapper script for catch22_jll for armv7l-linux-musleabihf
export ccatch22

JLLWrappers.@generate_wrapper_header("catch22")
JLLWrappers.@declare_library_product(ccatch22, "libcatch22.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        ccatch22,
        "lib/libcatch22.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using PNML: firstpage

header("PAGES")
@testset "pages" begin
    str = """
    <?xml version="1.0"?>
    <pnml xmlns="http://www.pnml.org/version-2009/grammar/pnml">
        <net id="net0" type="pnmlcore">
            <page id="page1">
                <place id="p1"/>
                <transition id ="t1"/>
                <arc id="a1" source="p1" target="t1"/>
                <arc id="a12" source="t1" target="rp1"/>
                <referencePlace id="rp1" ref="p2"/>
                <page id="page11">
                    <place id="p11" />
                    <page id="page111">
                        <place id="p111" />
                    </page>
                </page>
                <page id="page12" />
                <page id="page13" />
                <page id="page14" />
            </page>
            <page id="page2">
            <place id="p2"/>
            <transition id ="t2"/>
                <arc id="a2" source="t2" target="p2"/>
                <arc id="a22" source="t2" target="rp2"/>
                <referencePlace id="rp2" ref="p3"/>
                <referenceTransition id="rt2" ref="t3"/>
            </page>
            <page id="page3">
                <place id="p3"/>
                <transition id ="t3"/>
                <arc id="a3" source="t3" target="p4"/>
                <page id="page4">
                    <place id="p4"/>
                    <transition id ="t4"/>
                    <arc id="a4" source="t4" target="p1"/>
                    <page id="page41">
                        <place id="p41" />
                        <page id="page411">
                            <place id="p411" />
                        </page>
                    </page>
                    <page id="page42" />
                    <page id="page43" />
                    <page id="page44" />
                </page>
            </page>
        </net>
    </pnml>
    """
    model = parse_str(str)
    net = PNML.first_net(model)

    @test net isa PNML.PnmlNet
    @test typeof(net) <: PNML.PnmlNet
    @test typeof(firstpage(net)) <: PNML.Page

    printnode(net; label="\n----------------\n Multiple nested pages")

    @testset "flatten" begin
        PNML.flatten_pages!(net)
        printnode(net; label="\n----------------\n Flattened & dereferenced to 1 page")

        #@show PNML.arc_ids(net)
        #@show PNML.place_ids(net)
        #@show PNML.transition_ids(net)
        #@show PNML.refplace_ids(net)
        #@show PNML.reftransition_ids(net)
        #println()

        expected_a = [:a1, :a12, :a2, :a22, :a3, :a4]
        @test PNML.arc_ids(net) == expected_a
        @test PNML.arc_ids(firstpage(net)) == expected_a
        @test PNML.arc_ids(net) == PNML.arc_ids(firstpage(net))
        @test_call PNML.arc_ids(net)
        @test_call PNML.arc_ids(firstpage(net))

        for a ∈ expected_a
            @test a ∈ PNML.arc_ids(net)
        end

        expected_p = [:p1, :p11, :p111, :p2, :p3, :p4, :p41, :p411]
        @test PNML.place_ids(net) == expected_p
        @test PNML.place_ids(firstpage(net)) == expected_p
        @test PNML.place_ids(net) == PNML.place_ids(firstpage(net))
        @test_call PNML.place_ids(net)
        @test_call PNML.place_ids(firstpage(net))

        for p ∈ expected_p
            @test p ∈ PNML.place_ids(net)
        end

        expected_t = [:t1, :t2, :t3, :t4]
        @test PNML.transition_ids(net) == expected_t
        @test PNML.transition_ids(firstpage(net)) == expected_t
        @test PNML.transition_ids(net) == PNML.transition_ids(firstpage(net))
        @test_call PNML.transition_ids(net)
        @test_call PNML.transition_ids(firstpage(net))

        for t ∈ expected_t
            @test t ∈ PNML.transition_ids(net)
        end

        expected_rt = []#:rt2]
        @test PNML.reftransition_ids(net) == expected_rt
        @test PNML.reftransition_ids(firstpage(net)) == expected_rt
        @test PNML.reftransition_ids(net) == PNML.reftransition_ids(firstpage(net))
        @test_call PNML.reftransition_ids(net)
        @test_call PNML.reftransition_ids(firstpage(net))

        for rt ∈ expected_rt
            @test rt ∈ PNML.reftransition_ids(net)
        end

        expected_rp = []#:rp1, :rp2]
        @test PNML.refplace_ids(net) == expected_rp
        @test PNML.refplace_ids(firstpage(net)) == expected_rp
        @test PNML.refplace_ids(net) == PNML.refplace_ids(firstpage(net))
        @test_call PNML.refplace_ids(net)
        @test_call PNML.refplace_ids(firstpage(net))

        for rp ∈ expected_rp
            @test rp ∈ PNML.refplace_ids(net)
        end
    end
end # pages


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 1.0
x, z = .1, .1
b1 = Box(x, l1, z, l1 * 100, color = RGBA(1., 1., 0.))
b2 = Box(x, l1, z, l1 * 1, color = RGBA(1., 0., 0.))

vert11 = [0.;l1 / 2;0.]
vert12 = -vert11

# Links
origin = Origin{Float64}()
link1 = Body(b1)
link2 = Body(b2)

# Constraints
joint0to1e = EqualityConstraint(Floating(origin, link1))
# joint0to1 = EqualityConstraint(Fixed(origin,link1))
joint1to2e = EqualityConstraint(Revolute(link1, link2, ex; p1=vert12, p2=vert11))

links = [link1; link2]
constraints = [joint0to1e;joint1to2e]
shapes = [b1;b2]

function joint_torque_control!(mechanism, k)
    τ = SA[0.05]
    setForce!(mechanism, joint1to2e, τ * 0)
    return
end

mech = Mechanism(origin, links, constraints, shapes = shapes, g = 0.)
setPosition!(link1,link2,p1 = vert12,p2 = vert11)

# setForce!(mech, joint1to2e, [0.05])


text:
#=
This file contains main functions for calculating the expansion entropy
of dynamical systems.
Expansion entropy is defined in (Hunt and Ott, 2015) [1] as a quantitative measure
of chaos. For details, read the docstrings of the functions below.

[1] : B. & E. Ott, ‘Defining Chaos’, [Chaos 25.9 (2015)](https://doi.org/10/gdtkcf)
=#
export boxregion, expansionentropy, expansionentropy_batch, expansionentropy_sample
using LinearAlgebra
using Statistics

"""
    expansionentropy(ds::DynamicalSystem, sampler, restraining; kwargs...)

Calculate the expansion entropy[^Hunt2015] of `ds`, in the restraining region ``S`` defined by
`restraining`, by estimating the slope of the biggest linear region
of the curve ``\\log E_{t0+T, t0}(f, S)`` versus ``T`` (using [`linear_region`](@ref)).
This is an approximation of the expansion entropy ``H_0``, according to[^Hunt2015].

`sampler` is a 0-argument function that generates a random initial condition (a sample)
of `ds`. `restraining` is a 1-argument function `restraining(u)` that given the state
`u` it returns `true` if the state is inside the restraining region ``S``.

Use [`boxregion`](@ref) for an easy way to define `sampler` and `restraining` on a
multidimension box.

## Keyword Arguments
* `N = 1000` : Number of samples taken at each batch (same as ``N`` of [1]).
* `steps = 40` : The maximal steps for which the system will be run.
* `Ttr = 0` : Transient time to evolve each initial condition before starting to comute ``E``.
  This is `t0` of [1] and of the following notation.
* `batches = 100` : Number of batches to run the calculation, see below.
* `Δt = 1` : Integration step size.
* `diffeq` is a `NamedTuple` (or `Dict`) of keyword arguments propagated into
  `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.

## Description
`N` samples are initialized and propagated forwards in time (along with their tangent space).
At every time ``t`` in `[t0+Δt, t0+2dt, ... t0+steps*Δt]` we calculate ``H``:
```math
H[t] = \\log E_{t0+T, t0}(f, S),
```
with
```math
E_{t0+T, t0}(f, S) = \\frac 1 N \\sum_{i'} G(Df_{t0+t, t0}(x_i))
```
(using same notation as [^Hunt2015]). In principle ``E`` is the average largest possible
growth ratio within the restraining region (sampled by the initial conditions).
The summation is only over ``x_i`` that stay inside the region ``S``
defined by the boolean function `restraining`.
This process is done by the [`expansionentropy_sample`](@ref) function.

Then, this is repeated for `batches` amount of times, as recommended in[^Hunt2015].
From all these batches, the mean and std of ``H`` is computed at every time point.
This is done by the [`expansionentropy_batch`](@ref) function.
When plotted versus ``t``, these create the curves and error bars of e.g. Figs 2, 3 of [1].

This function `expansionentropy` simply returns the slope of the biggest linear region
of the curve ``H`` versus ``t``, which approximates the expansion entropy ``H_0``.
It is therefore *recommended* to use [`expansionentropy_batch`](@ref) directly and
evaluate the result yourself, as this step is known to be inaccurate for
non-chaotic systems (where ``H`` fluctuates strongly around 0).

[^Hunt2015]: B. Hunt & E. Ott, ‘Defining Chaos’, [Chaos 25.9 (2015)](https://doi.org/10/gdtkcf)
"""
function expansionentropy(system, sampler, restraining; kwargs...)
    times, means, stds = expansionentropy_batch(system, sampler, restraining; kwargs...)
    if any(isnan, stds)
        i = findfirst(isnan, stds)
        @warn "All (or all except one) samples have escaped the given region at time = $(times[i])."
        times = times[1:i-1]
        means = means[1:i-1]
    end
    _, slope = linear_region(times, means)
    return slope
end

"""
    boxregion(as, bs) -> sampler, restraining

Define a box in ``\\mathbb{R}^d`` with edges the `as` and `bs` and then
return two functions: `sampler`, which generates a random initial condition in that box
and `restraining` that returns `true` if a given state is in the box.
"""
function boxregion(as, bs)
    @assert length(as) == length(bs) > 0
    gen() = [rand()*(bs[i]-as[i]) + as[i] for i in 1:length(as)]
    restraining(x) = all(as .< x .< bs)
    return gen, restraining
end

# Specialized 1-d version
function boxregion(a::Real, b::Real)
    a, b = extrema((a, b))
    gen() = rand()*(b-a) + a
    restraining = x -> a < x < b
    return gen, restraining
end

#####################################################################################
# Actual implementation of expansion entropy
#####################################################################################
"""
    maximalexpansion(M)

Calculates the maximal expansion rate of M,
i.e. the product of all singular values of M that are greater than 1. In the
notation of [1], it is the function ``G``.
"""
maximalexpansion(M) = prod(filter(x -> x > 1.0, svdvals(M)))

"""
    expansionentropy_sample(ds, sampler, restraining; kwargs...)

Return `times, H` for one sample of `ds` (see [`expansionentropy`](@ref)).
Accepts the same argumets as `expansionentropy`, besides `batches`.
"""
function expansionentropy_sample(system::DynamicalSystem, sampler, restraining;
    N=1000, steps=40, Δt=1, Ttr=0, diffeq = NamedTuple(), kwargs...)

    if !isempty(kwargs)
        @warn DIFFEQ_DEP_WARN
        diffeq = NamedTuple(kwargs)
    end

    D = dimension(system)
    M = zeros(steps)
    # M[t] will be Σᵢ G(Dfₜ₀,ₜ₀₊ₜ(xᵢ))
    # The summation is over all sampled xᵢ that stay inside S during [t0, t0 + t].

    times = @. (system.t0+Ttr)+Δt*(1:steps)
    t_identity = SMatrix{D, D, Float64}(I)
    t_integ = tangent_integrator(system; diffeq)
    Ttr > 0 && (u_integ = integrator(system))

    for i ∈ 1:N
        u = sampler() # New sample point.
        if Ttr > 0 # Evolve through transient time.
            reinit!(u_integ, u)
            step!(u_integ, Ttr, true) # stepping must be exact here for correct time vector
            u = u_integ.u
        end
        reinit!(t_integ, u, t_identity; t0 = t_integ.t0 + Ttr)
        for i ∈ 1:steps # Evolve the sample point for the duration [t0, t0+steps*Δt]
            step!(t_integ, Δt, true)
            u = get_state(t_integ)
            !restraining(u) && break # Stop the integration if the orbit leaves the region.
            Df = get_deviations(t_integ)
            M[i] += maximalexpansion(Df)
        end
    end
    return times, log.(M./N)
end

# This version only deals with 1-dimensional discrete dynamical systems, but does
# it really fast, by avoiding `tangent_integrator` and `maximalexpansion`.
function expansionentropy_sample(system::DiscreteDynamicalSystem{IIP, S, 1}, sampler, restraining;
    N=1000, steps=40, Δt=1, Ttr=0, kwargs...) where {IIP, S}
    f = system.f
    p = system.p
    jacob = system.jacobian
    t0 = system.t0
    times = @. (t0+Ttr)+Δt*(1:steps)

    Δt = max(Int(floor(Δt)), 1)
    Ttr = max(Int(floor(Ttr)), 0)
    M = zeros(steps)

    for i ∈ 1:N
        x = sampler()
        t = t0

        # Evolve for a transient time.
        for _ ∈ 1:Ttr
            x = f(x, p, t)
            t += 1
        end

        Df = 1.0
        for step ∈ 1:steps # Evolve point x for `steps` steps.
            for _ ∈ 1:Δt # Evolve Δt steps at a time.
                Df = jacob(x, p, t) * Df
                x = f(x, p, t)
                t += 1
            end

            !restraining(x) && break
            M[step] += max(1.0, abs(Df))
        end
    end
    return times, log.(M./N)
end

"""
    expansionentropy_batch(ds, sampler, restraining; kwargs...)

Run [`expansionentropy_sample`](@ref) `batch` times, and return
`times, mean(H), std(H)` for all resulting `H`, see [`expansionentropy`](@ref).

Accepts the same arguments as `expansionentropy`.
"""
function expansionentropy_batch(system, sampler, restraining; batches=100, steps=40, kwargs...)
    # TODO: It is a mistake that `expansionentropy_batch` doen't create an integrator that
    # us just reinited inside `expansionentropy_sample`...
    means = fill(NaN, steps)
    stds = fill(NaN, steps)
    eesamples = zeros(batches, steps)
    # eesamples[k, t] = The k-th sample of expansion entropy from t0 to (t0 + t)

    times = undef
    # Collect all the samples
    for k in 1:batches
        times, eesamples[k, :] = expansionentropy_sample(
            system, sampler, restraining; steps, kwargs...
        )
    end

    # Calculate the mean and standard deviations
    for t in 1:steps
        entropysamples = filter(isfinite, @view eesamples[:, t])
        # remove -Inf entries, which indicate all samples failed to stay inside the given region.
        if length(entropysamples) ≤ 1
        @warn "All (or all except one) samples have escaped the given region. "*
        "Consider increasing sample or batch number. Terminating at step = $(t)."
            break
        end
        means[t] = mean(entropysamples)
        stds[t] = std(entropysamples)
    end
    return times, means, stds
end


text:
function gausslegendre( n::Int32 )
# GAUSSLEGENDRE(n)  COMPUTE THE GAUSS-LEGENDRE NODES AND WEIGHTS IN O(n) time.

if n < 0
    x = (Float64[],Float64[])
elseif n == 0
    x = (Float64[],Float64[])
elseif n == 1
    x = ([0.0],[2.0])
elseif n == 2
    x = Array(Float64,2)
    w = Array(Float64,2)
    x[1] = -1/sqrt(3)
    x[2] = 1/sqrt(3)
    w[1] = 1.0
    w[2] =  1.0
    x = (x,w)
elseif n == 3
    x = Array(Float64,3)
    w = Array(Float64,3)
    x[1] = -sqrt(3/5)
    x[2] = 0.0
    x[3] = sqrt(3/5)
    w[1] = 5/9
    w[2] = 8/9
    w[3] = 5/9
    x = (x,w)
elseif n == 4
    x = Array(Float64,4)
    w = Array(Float64,4)
    const a::Float64 = 2/7*sqrt(6/5)
    x[1] = -sqrt(3/7+a)
    x[2] = -sqrt(3/7-a)
    x[3] =  sqrt(3/7-a)
    x[4] =  sqrt(3/7+a)
    w[1] = (18 - sqrt(30))/36;
    w[2] = (18 + sqrt(30))/36;
    w[3] = (18 + sqrt(30))/36;
    w[4] = (18 - sqrt(30))/36;
    x = (x, w)
elseif n == 5
    x = Array(Float64, 5)
    w = Array(Float64, 5)
    const b::Float64 = 2sqrt(10/7)
    x[1] = -sqrt(5+b)/3
    x[2] = -sqrt(5-b)/3
    x[3] = 0.0
    x[4] = sqrt(5-b)/3
    x[5] = sqrt(5+b)/3
    w[1] = (322-13sqrt(70))/900
    w[2] = (322+13sqrt(70))/900
    w[3] = 128/225
    w[4] = (322+13sqrt(70))/900
    w[5] = (322-13sqrt(70))/900
    x = (x, w);
elseif n <= 60
# NEWTON'S METHOD WITH THREE-TERM RECURRENCE:
    x = rec( n::Int32 )
else
# USE ASYMPTOTIC EXPANSIONS:
    x = asy( n::Int32 )
end
    return x   # nodes and weights in tuple.
end

function asy( n::Int32 )
# COMPUTE GAUSS-LEGENDRE NODES AND WEIGHTS USING ASYMPTOTIC EXPANSIONS. COMPLEXITY O(n).
    # Nodes and weights:
    m = (mod(n,2)==0) ? n>>1 : (n+1)>>1
    a = besselZeroRoots(m)
    scale!(a,1/(n + 0.5))
    x = legpts_nodes(n, a);
    w = legpts_weights(n, m, a)
    # Use symmetry to get the others:
    if ( mod(n, 2) == 1 )
        x = [-x[1:end-1] ; 0.0 ; x[end-1:-1:1] ]
        w = [w ; w[end-1:-1:1] ]
    else
        x = [-x ; x[end:-1:1] ]
        w = [w ; w[end:-1:1] ]
    end
    return x, w
end

function legpts_nodes(n::Int32, a::Array{Float64})
# ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE NODES.
    vn = 1/(n + 0.5)
    m = length(a)
    u = Array(Float64,m); nodes = Array(Float64,m)
    for jj=1:m u[jj] = cot(a[jj]) end
    for jj=1:m nodes[jj] = a[jj] + ((u[jj]-1/a[jj])/8*vn^2) end
    nodes += ((n <= 3950) ? (6(1+u.^2)./a + 25./a.^3 - u.*(31u.^2+33))/384 : 0.0)*vn^4 +
            ((n <= 255) ? u.*(2595 + 6350*u.^2 + 3779*u.^4)/15360 -1073/5120./a.^5 +
            (1+u.^2).*(-(31*u.^2 + 11)/1024./a + u/512./a.^2 + -25/3072./a.^3) : 0.0)*vn^6
    for jj=1:m nodes[jj] = cos(nodes[jj]) end
    return nodes
end

function legpts_weights(n::Int32, m::Int32, a::Array{Float64})
# ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE WEIGHTS.
    vn = 1/(n + 0.5);
    u = Array(Float64,m); ua = Array(Float64,m); weights = Array(Float64,m)
    for jj=1:m u[jj] = cot(a[jj]) end
    for jj=1:m ua[jj] = u[jj]*a[jj] end
     # Assemble coefficients:
    W1 = Array(Float64,m); W2 = Array(Float64,m); W3 = Array(Float64,m)
    for jj=1:m W1[jj] = (n <= 850000) ? (ua[jj] + a[jj]^2 - 1)/a[jj]^2/8 : 0.0 end
    for jj=1:m W2[jj] = (n <= 1500) ? ( 81 - 31*ua[jj] - 3*(1-2*u[jj]^2)*a[jj]^2 + 6*u[jj]*a[jj]^3 -
                                        (27 + 84*u[jj]^2 + 56*u[jj]^4)*a[jj]^4 )/a[jj]^4/384 : 0.0 end
    for jj=1:m W3[jj] = ( n <= 170 ) ? (187/96*u[jj]^4 + 295/256*u[jj]^2 + 151/160*u[jj]^6 + 153/1024) +
                                    (-119/768*u[jj]^2 -35/384*u[jj]^4 - 65/1024)*u[jj]/a[jj] +
                                    (5/512 + 7/384*u[jj]^4 + 15/512*u[jj]^2)/a[jj]^2 +
                                    (u[jj]^3/512 - 13/1536*u[jj])/a[jj]^3 +
                                    (-7/384*u[jj]^2 + 53/3072)/a[jj]^4 +
                                    (3749/15360*u[jj])/a[jj]^5 -1125/1024/a[jj]^6 : 0.0 end
    bJ1 = besselJ1(m)
    for jj=1:m weights[jj] = 2/(bJ1[jj]*(a[jj]/sin(a[jj]))*(1/vn^2 + W1[jj] + W2[jj]*vn^2 + W3[jj]*vn^4)) end
    return weights
end

function rec( n::Int32 )
# COMPUTE GAUSS-LEGENDRE NODES AND WEIGHTS USING NEWTON'S METHOD. THREE-TERM RECURENCE
# IS USED FOR EVALUATION. COMPLEXITY O(n^2).
    hN = mod(n,2)
    # Initial guesses:
    x0 = asy(n)[1]; x = x0[(n-hN)/2+1:n]
    # Perform Newton to find zeros of Legendre polynomial:
    PP = innerRec( n, x ); dx = -PP[1]./PP[2]; x += dx
    # One more Newton for derivatives:
    PP = innerRec( n, x ); dx = -PP[1]./PP[2]; x += dx
    #PP = (n<45)? innerRec( n, x ) : PP
    # Use symmetry to get the other Legendre nodes and weights:
    nodes = vcat( -x[(n+hN)/2:-1:hN+1] , x )
    weights = 2./((1-nodes.^2).*vcat( PP[2][(n+hN)/2:-1:1+hN] , PP[2] ).^2)
    return nodes, weights
end

function innerRec( n::Int32, x )
# EVALUATE LEGENDRE AND ITS DERIVATIVE USING THREE-TERM RECURRENCE RELATION.
    N = size(x,1)
    myPm1 = Array(Float64,N); myPPm1 = Array(Float64,N)
    for j = 1:N
        xj = x[j]; Pm2 = 1.0; Pm1 = xj; PPm1 = 1.0; PPm2 = 0.0
        for k = 1:n-1
            Pm2, Pm1 = Pm1, ((2k+1)*Pm1*xj - k*Pm2)/(k+1)
            PPm2, PPm1 = PPm1, ((2k+1)*(Pm2 + xj*PPm1)-k*PPm2)/(k+1)
        end
        myPm1[j] = Pm1
        myPPm1[j] = PPm1
    end
    return myPm1, myPPm1
end

function besselZeroRoots(m::Int32)
#BESSEL0ROOTS ROOTS OF BESSELJ(0,x). USE ASYMPTOTICS.
# Use McMahon's expansion for the remainder (NIST, 10.21.19):
    jk = Array(Float64,m); ak = Array(Float64,m); ak82 = Array(Float64,m)
    p = [1071187749376/315 0.0 -401743168/105 0.0 120928/15 0.0 -124/3 0.0 1.0 0.0]
    for jj=1:m ak[jj] = pi*( jj -.25 ) end
    for jj=1:m ak82[jj] = (.125/ak[jj])^2 end
    # First 20 are precomputed:
    twenty = [2.4048255576957728, 5.5200781102863106, 8.6537279129110122,
            11.791534439014281, 14.930917708487785, 18.071063967910922,
            21.211636629879258, 24.352471530749302, 27.493479132040254,
            30.634606468431975, 33.775820213573568, 36.917098353664044,
            40.058425764628239, 43.199791713176730, 46.341188371661814,
            49.482609897397817, 52.624051841114996, 55.765510755019979,
            58.906983926080942, 62.048469190227170]
    for jj=1:min(m,20) jk[jj] = twenty[jj] end
    for jj=21:min(m,47) jk[jj] = ak[jj] + .125/ak[jj]*(1 + ak82[jj]*(p[7] + ak82[jj]*(p[5] + ak82[jj]*p[3]))) end
    for jj=48:min(m,344) jk[jj] = ak[jj] + .125/ak[jj]*(1 + ak82[jj]*(p[7]+ ak82[jj]*p[5])) end
    for jj=345:min(m,13191) jk[jj] = ak[jj] + .125/ak[jj]*(1 + ak82[jj]*p[7]) end
    for jj=13192:m jk[jj] = ak[jj] + .125/ak[jj] end
    return jk
end

function besselJ1( m::Int32 )
# BESSELJ1 EVALUATE BESSELJ(1,x)^2 AT ROOTS OF BESSELJ(0,x). USE ASYMPTOTICS.
# Use Taylor series of (NIST, 10.17.3) and McMahon's expansion (NIST, 10.21.19):
    Jk2 = Array(Float64,m); ak = Array(Float64,m); ak2 = Array(Float64,m)
    c = [-171497088497/15206400, 461797/1152, -172913/8064, 151/80, -7/24, 0.0, 2.0]
    for jj=1:m ak[jj] = pi*(jj-.25) end
    for jj=1:m ak2[jj] = (1/ak[jj])^2 end
    # First 10 are precomputed:
    ten = [ 0.2695141239419169, 0.1157801385822037, 0.07368635113640822,
            0.05403757319811628, 0.04266142901724309, 0.03524210349099610,
            0.03002107010305467, 0.02614739149530809, 0.02315912182469139, 0.02078382912226786]
    for jj=1:min(m,10) Jk2[jj] = ten[jj] end
    for jj=11:min(m,15) Jk2[jj] = 1/(pi*ak[jj])*(c[7] + ak2[jj]^2*(c[5] + ak2[jj]*(c[4] +
                ak2[jj]*(c[3] + ak2[jj]*(c[2]+ak2[jj]*c[1]))))) end
    for jj=16:min(m,21) Jk2[jj] = 1/(pi*ak[jj])*(c[7] + ak2[jj]^2*(c[5] + ak2[jj]*(c[4] + ak2[jj]*(c[3] + ak2[jj]*c[2])))) end
    for jj=22:min(m,55) Jk2[jj] =  1/(pi*ak[jj])*(c[7] + ak2[jj]^2*(c[5] + ak2[jj]*(c[4] + ak2[jj]*c[3]))) end
    for jj=56:min(m,279) Jk2[jj] = 1/(pi*ak[jj])*(c[7] + ak2[jj]^2*(c[5] + ak2[jj]*c[4])) end
    for jj=280:min(m,2279) Jk2[jj] = 1/(pi*ak[jj])*(c[7] + ak2[jj]^2*c[5]) end
    for jj=2280:m Jk2[jj] = 1/(pi*ak[jj])*c[7] end
    return Jk2
end


text:
using IcanProj
using DataFrames
using PyPlot

"""
Load input variables directly from the forcing files.
"""
function load_from_forcings(watersheds, resolutions, opt)

    meteo_wsh = Dict()

    for watershed in watersheds

        meteo_res = Dict()

        for resolution in resolutions

            # Read soil parameter file with data about latitude and longitude

            path_param = joinpath(opt["target_folder"], watershed, resolution)

            soil_param = read_soil_params(path_param)

            # Path to model forcings for the given resolution

            if resolution == "1km"
                path_forcing = joinpath(opt["base_folder"])
            else
                path_forcing = joinpath(opt["target_folder"], watershed, resolution)
            end

            # Read one forcing file to get the size

            row = soil_param[1, :]

            lat = @sprintf("%0.5f", row[:lat][1])
            lon = @sprintf("%0.5f", row[:lon][1])

            file_src = joinpath(path_forcing, "forcing/data_$(lat)_$(lon)")

            prec, tmin, tmax, wind = read_vic_forcing(file_src)

            # Allocate arrays

            prec_array = zeros(Float64, length(prec), size(soil_param,1))
            tmin_array = zeros(Float64, length(prec), size(soil_param,1))
            tmax_array = zeros(Float64, length(prec), size(soil_param,1))
            wind_array = zeros(Float64, length(prec), size(soil_param,1))

            for i in 1:size(soil_param,1)

                row = soil_param[i, :]

                lat = @sprintf("%0.5f", row[:lat][1])
                lon = @sprintf("%0.5f", row[:lon][1])

                file_src = joinpath(path_forcing, "forcing/data_$(lat)_$(lon)")

                prec, tmin, tmax, wind = read_vic_forcing(file_src)

                prec_array[:, i] = prec
                tmin_array[:, i] = tmin
                tmax_array[:, i] = tmax
                wind_array[:, i] = wind

            end

            # Compute summary statistics

            meteo = DataFrame()

            meteo[:prec_mean] = mean(prec_array, 2)[:,1]
            meteo[:tmin_mean] = mean(tmin_array, 2)[:,1]
            meteo[:tmax_mean] = mean(tmax_array, 2)[:,1]
            meteo[:wind_mean] = mean(wind_array, 2)[:,1]

            meteo[:prec_std] = std(prec_array, 2)[:,1]
            meteo[:tmin_std] = std(tmin_array, 2)[:,1]
            meteo[:tmax_std] = std(tmax_array, 2)[:,1]
            meteo[:wind_std] = std(wind_array, 2)[:,1]

            meteo[:prec_range] = maximum(prec_array, 2)[:,1] - minimum(prec_array, 2)[:,1]
            meteo[:tmin_range] = maximum(tmin_array, 2)[:,1] - minimum(tmin_array, 2)[:,1]
            meteo[:tmax_range] = maximum(tmax_array, 2)[:,1] - minimum(tmax_array, 2)[:,1]
            meteo[:wind_range] = maximum(wind_array, 2)[:,1] - minimum(wind_array, 2)[:,1]

            meteo_res[resolution] = meteo

        end

        meteo_wsh[watershed] = meteo_res

    end

    return meteo_wsh

end

"""
Plot input data characteristics against different resolutions.
"""
function plot_elevation_meteo(meteo_wsh, watersheds, resolutions, var_name, ylabel_name, file_name)

    ioff()

    fig = plt[:figure](figsize = (8, 6))

    for watershed in watersheds

        var = Float64[]

        for resolution in resolutions

            push!(var, mean(meteo_wsh[watershed][resolution][var_name]))

        end

        var[isnan.(var)] = 0

        plot(var, label = wsh_name[watershed])

    end

    legend()

    plt[:xticks](collect(0:length(resolutions)-1), resolutions)
    plt[:xlabel]("Spatial resolution")
    plt[:ylabel](ylabel_name)

    file_name = joinpath(Pkg.dir("IcanProj", "plots", "$(file_name).png"))

    rm(file_name, force=true)
    savefig(file_name, dpi = 300)

    close(fig)

end

# Get settings

opt = get_options()

wsh_info = get_wsh_info()

wsh_name = get_wsh_name(wsh_info)

resolutions = opt["resolutions"]

watersheds = opt["stat_sel"]

# Load inputs from forcing files

meteo_wsh = load_from_forcings(watersheds, resolutions, opt)

# Plot the different variables

colnames = names(meteo_wsh[watershed][resolution])

for name in colnames

    var_name = name
    ylabel_name = String(name)
    file_name = String(name)

    plot_elevation_meteo(meteo_wsh, watersheds, resolutions, var_name, ylabel_name, file_name)

end


text:
# Load Julia packages (libraries)

using Pkg, DrWatson

#@quickactivate "StatisticalRethinkingStan"
using StanSample
using StatisticalRethinking

# Define the Stan language model

stan8_2s = "
data{
    int N;
    vector[N] y;
}
parameters{
    real mu;
    real<lower=0> sigma;
}
model{
    y ~ normal( mu , sigma );
}
";

# Define the Stanmodel.

m8_2s = SampleModel("m8.2s", stan8_2s);

# Input data for cmdstan

m8_2_data = Dict("N" => 2, "y" => [-1, 1]);
m8_2_init = Dict("mu" => 0.0, "sigma" => 1.0);

# Sample using cmdstan

rc8_2s = stan_sample(m8_2s; data=m8_2_data, init=m8_2_init);

# Describe the draws
if success(rc8_2s)
  part8_2s = read_samples(m8_2s, :particles)
  part8_2s |> display
end

# End of m8_2s


text:
# Version

version() = qiskit.__version__

qiskit_version() = qiskit.__qiskit_version__


text:
module ROyFit

using KM3NeT

const n = 1.3797
const c = 299792458.0 / 1e9

function prefit(hits::Vector{Hit}, calib::Calibration; min_n_doms = 3, tmax = 15, n_hits = 4)
    shits = nfoldhits(hits, tmax, n_hits)
    n_doms = length(Set(sh.dom_id for sh in shits))
    if n_doms < min_n_doms
        return NoRecoTrack()
    end
    chits = calibrate(shits, calib)
    sort!(chits)
    t0 = chits[div(length(shits), 2)].t
    pos, dir = svdfit(matrix([h.pos for h in chits]))
    if (last(chits).pos - first(chits).pos) ⋅ dir  < 0.0
        dir *= -1
    end
    return RecoTrack(dir, pos, t0)
end

struct SingleDUParameters
    t::Vector{HitTime}
    z::Vector{Float32}
    tot::Vector{ToT}
    σ_t::Float32
    mean_tot::Float32

    function SingleDUParameters(t, z, tot, σ_t)
        mean_tot = mean(tot)
        return new(t, z, tot, σ_t, mean_tot)
    end
end

function γ_d(p::SingleDUParameters, uz, zc, dc)
    return n/sqrt(n^2-1) * sqrt(dc^2 + ((p.z-zc)^2) * (1-uz^2))
end

function γ_t(p::SingleDUParameters, uz, zc, dc, tc)
    return tc + ((p.z-zc)*uz+(n^2 - 1) * γ_d(p, uz, zc, dc) / n ) / c
end

function cosθ(p::SingleDUParameters, uz, zc, dc)
    return (1 - uz^2) * (p.z - zc) / γ_D(p, uz, zc, dc) + uz/n
end

function make_quality_function(p::SingleDUParameters)
    function quality_function(uz, zc, dc, tc)
        tot_weights = 2p.tot / cosθ(p, uz, zc, dc) + 1.0
        mean_tot = mean(tot_weights)
        d = sqrt(p.d1^2 + γ_d(p, uz, zc, dc)^2)
        #return sum(γ_t(p, uz, zc, dc, tc) - p.t)^2 / p.σ_t^2 +
    end
    return quality_function
end

"""
    function write_prefit_summary(fname::AbstractString, detx_fname::AbstractString)

Perform the prefit on a given file and create a CSV summary file.
"""
function write_prefit_summary(fname::AbstractString, detx_fname::AbstractString)
    run_id = mc_run_id(fname)
    calib = read_calibration(detx_fname)
    reader = KM3NeT.EventReader(fname, load_tracks=true)
    output_fname = basename(fname) * ".royprefit.csv"
    out_fobj = open(output_fname, "w")
    write(out_fobj, "run_id,event_id,n_muons,muon_energy,dir_x,dir_y,dir_z,alpha\n")
    for (event_id, event) in enumerate(reader)
        n_muons = length([t for t in event.mc_tracks if t.particle_type == 5])
        muon_energy = sum([t.E for t in event.mc_tracks if t.particle_type == 5])
        write(out_fobj, "$run_id,$event_id,$n_muons,$muon_energy,")
        reco = ROyFit.prefit(event.hits, calib)
        if isa(reco, NoRecoTrack)
            write(out_fobj, "nan,nan,nan,nan\n")
            continue
        end
        muon = event.mc_tracks[2]
        alpha = angle(Direction(reco.dir), muon.dir)
        write(out_fobj, "$(reco.dir.x),$(reco.dir.y),t$(reco.dir.z),$alpha\n")
    end
    close(out_fobj)
end

end # module


text:
let
    #= line 1 =#
    var"|Println|b_unit"::Symbol = Symbol("unit@3_|Println|c")::Symbol
    var"|Println|d_println"::Function = println::Function
    (var"|Println|d_println"::Function)("hello!"::String)::Any
    var"|Println|e"::Function = function (_,)
            #= ../examples/println.mlfsf:1 =#
            var"|Println|d_println"
        end
    var"|Println|f"::Function = function (_,)
            #= ../examples/println.mlfsf:1 =#
            Symbol("unit@3_|Println|c")
        end
    var"|Println|"::Any = Symbol("Module Println")
    Println::Any = Symbol("namespace {Module Println, {{unit, |Println|b_unit, Type unit@3_|Println|c}, {println, |Println|d_println, forall a@10546255710782205774.a@10546255710782205774->unit@3_|Println|c}}}")
    0
end


text:
using Plots

include("models/Dicke.jl")
include("modules/ClassicalDynamics.jl")

dickeParameters = [2,0.5,1,1]
energy = 50
numTrajectories = 5
PSPoints = 1000

savePath = "d:\\results\\Dicke"

pyplot(size = (1800,1000))
fig = PoincareSection(energy, dickeParameters, numTrajectories; sectionPlane=1, maxPSPoints=PSPoints, verbose=true, maxICNumber=10000, min=-sqrt(2.0*energy), max=sqrt(2.0*energy))
savefig(fig, savePath * "_PS_E=$(energy)_$dickeParameters.png")

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
using Documenter, Caesar
using RoME

import IncrementalInference: fmcmc!, localProduct, prodmultiplefullpartials, prodmultipleonefullpartials, setfreeze!
import IncrementalInference: cliqGibbs, packFromLocalPotentials!, treeProductDwn, updateFGBT!, upGibbsCliqueDensity
import IncrementalInference: initfg, downGibbsCliqueDensity
import IncrementalInference: solveGraphParametric, solveGraphParametric!

using KernelDensityEstimatePlotting
# import KernelDensityEstimatePlotting: plotKDE
using RoMEPlotting

# until namespaces are properly figured out
using DistributedFactorGraphs
import DistributedFactorGraphs: showFactor, showVariable
import DistributedFactorGraphs: deleteVariable!

makedocs(
    modules = [Caesar, RoME, IncrementalInference, RoMEPlotting, KernelDensityEstimatePlotting, DistributedFactorGraphs],
    format = Documenter.HTML(),
    sitename = "Caesar.jl",
    pages = Any[
        "Welcome" => "index.md",
        "Introduction" => [
            "Introduction" => "introduction.md",
            "Gaussian vs. Non-Gaussian" => "concepts/why_nongaussian.md",
            "Installation" => "installation_environment.md",
            "Using Julia" => "concepts/using_julia.md",
            "FAQ" => "faq.md",
        ],
        "Getting Started" => [
            "Initial Concepts" => "concepts/concepts.md",
            "Building Graphs" => "concepts/building_graphs.md",
            "Solving Graphs" => "concepts/solving_graphs.md",
            "Interact w Graphs" => "concepts/interacting_fgs.md",
            "Multi-Modal/Hypothesis" => "concepts/dataassociation.md",
            "Parallel Processing" => "concepts/parallel_processing.md",
            "[DEV] Parametric Solve" => "examples/parametric_solve.md",
        ],
        "Examples" => [
            "Caesar Examples" => "examples/examples.md",
            "Canonical 1D Example" => "examples/basic_continuousscalar.md",
            "Underconstrained Range-only" => "examples/basic_slamedonut.md",
            "Hexagonal 2D SLAM" => "examples/basic_hexagonal2d.md",
            "Fixed-Lag Solving 2D" => "examples/interm_fixedlag_hexagonal.md",
            "Dead Reckon Tether" => "examples/deadreckontether.md",
        ],
        "Graph Library" => [
            "Canonical Generators" => "examples/canonical_graphs.md",
            "Entry=>Data Blob" => "concepts/entry_data.md",
            "Variables/Factors" => "concepts/available_varfacs.md",
            "Flux (NN) Factors" => "concepts/flux_factors.md",
            "Images and AprilTags" => "examples/using_images.md",
        ],
        "Visualization" => [
            "Installing Viz" => "install_viz.md",
            "Plotting (2D)" => "concepts/2d_plotting.md",
            "Visualization (3D)" => "concepts/arena_visualizations.md",
        ],
        "Middlewares" => [
            "ROS Middleware" => "examples/using_ros.md",
            "Compile Binaries" => "concepts/compile_binary.md",
            "Zero Install Solution" => "concepts/zero_install.md",
            "Multi-session/agent Solving" => "concepts/multisession.md",
            "Multi-Language Support" => "concepts/multilang.md",
        ],
        "How to Expand?" => [
            "Pkg Framework" => "caesar_framework.md",
            "Custom Variables" => "examples/custom_variables.md",
            "Custom Prior Factor" => "examples/basic_definingfactors.md",
            "Custom Relative Factor" => "examples/custom_relative_factors.md",
            "Important Factor Features" => "examples/custom_factor_features.md",
            "Variable/Factor Considerations" => "examples/adding_variables_factors.md",
            "More Functions" => "func_ref.md",
        ],
        "Principles" => [
            "Filters vs. Graphs" => "principles/filterCorrespondence.md",
            "Generic Convolutions" => "principles/approxConvDensities.md",
            "Multiplying Functions (.py)" => "principles/multiplyingDensities.md",
            "Bayes (Junction) tree" => "principles/bayestreePrinciples.md",
            "Advanced Bayes Tree Topics" => "principles/initializingOnBayesTree.md",
            "Non-Gaussian Algorithm" => "concepts/mmisam_alg.md",
        ],
        "Developer Zone" => [
            "Wiki Pointers" => "dev/wiki.md",
            "Legacy Factors" => "examples/legacy_deffactors.md",
            "Creating DynPose Factor" => "principles/interm_dynpose.md",
            "Known Issue List" => "dev/known_issues.md",
            "Internal Functions" => "dev/internal_fncs.md",
            ],
        "Literature" => [
            "References" => "refs/literature.md"
        ],
    ]
    # html_prettyurls = !("local" in ARGS),
    )

deploydocs(
    repo   = "github.com/JuliaRobotics/Caesar.jl.git",
    target = "build"
)


text:
using HttpServer
using WebSockets
using JSON

decodeMessage(msg) = bytestring(msg)

global jsonwebsocketconnections = Dict{Int,WebSocket}()

function send(msg)
  global jsonwebsocketconnections
  jsonmsg = JSON.json(msg)
  if length(jsonwebsocketconnections) == 0
    println("No one to send to! ", jsonmsg)
  end
  for (id, client) in jsonwebsocketconnections
    println("Sending to ", id, ": ", jsonmsg)
    write(client, jsonmsg)
  end
end

function start_websocket_handler_while_processing(processingFunc; mainpage = "index.html", port = 8084)

  onepage = readall(mainpage)
  httph = HttpHandler() do req::Request, res::Response
    Response(onepage)
  end

  @async begin
    while true
      processingFunc()
    end
  end

  wsh = WebSocketHandler() do req, client
    global jsonwebsocketconnections
    println("Got request from client: ", client)
    jsonwebsocketconnections[client.id] = client
    while true
      msg = decodeMessage(read(client))
      println("Received from ", client.id, ": ", msg)
      try
        jsonmsg = JSON.parse(msg)
        handle_json_message(req, client, jsonmsg)
      catch error
        error("Could not parse and/or handle the msg: ", msg)
      end
    end
  end
  println("Now starting server!")
  server = Server(httph, wsh)
  run(server, 8084)
end

global currentrow = 1

# For now the processing is just to send some random data every now and then...
function processing_func()
  global currentrow
  println("Starting processing of row ", currentrow)
  numcols = 10
  @sync @parallel for col in 1:numcols
    sleep(5 + 10 * rand()) # Fake that processing takes time...
    send({"matrix" => "A", "action" => "set_value",
      "value" => rand(), "row" => currentrow, "col" => col})
  end
  currentrow += 1
end

start_websocket_handler_while_processing(processing_func)