text:
## This file is adapted from the code: https://www.mathworks.com/matlabcentral/fileexchange/53796-truncated-normal-and-student-s-t-distribution-toolbox
# Computes logarithm of tail of Z~N(0,1) mitigating numerical roundoff errors
lnPhi(x) = -0.5 * x^2 - 0.69314718055994530941723212 +
            log(erfcx(x / 1.4142135623730950488016887242))

sqrt2pi = sqrt(2*pi)

# Computes ln(P(a<Z 0) && return pa + log1p(-exp(pb - pa)) # case b > a > 0
    (b < 0) && return pb + log1p(-exp(pa - pb)) # case a < b < 0
    # case a < 0 < b
    log1p(-erfc(-a / 1.4142135623730950488016887242) / 2 -
           erfc(b / 1.4142135623730950488016887242) / 2)
end

# Sample random number from truncated normal distribution with l < x < u
trandn(rng::AbstractRNG, l, u) = rand(rng, Truncated(Normal(), l, u), 1)[1]

# Generates the proposals from the exponentially tilted sequential importance sampling pdf
# input: n is the number of samples
#        L is ?
#        l is the lower bound
#        u is the upper bound
#        mu is ?
# output:    'p', log-likelihood of sample
#             Z, random sample
function mvnrnd(rng::AbstractRNG, n, L, l, u, mu)
    d = length(l)
    p, Z = 0, zeros(d,n)
    for k = 1:d
        # compute matrix multiplication L*Z)
        col = dropdims(L[k:k, 1:k]*Z[1:k, :], dims = 1)

        # compute limits of truncation
        tl = l[k] .- mu[k] .- col
        tu = u[k] .- mu[k] .- col

        # simulate N(mu,1) conditional on [tl,tu]
        Z[k,:] .= mu[k] .+ trandn.(rng, tl, tu)

        # update likelihood ratio
        p = p .+ lnNpr.(tl, tu) .+ 0.5*mu[k]^2 .- mu[k]*Z[k,:]
    end
    p, Z
end

# Implements psi(x,mu). assumes scaled 'L' without diagonal
function psy(x, L, l, u, mu)
    d = length(u)

    # compute  ~l and ~u
    c = L * x
    nl = l - mu - c
    nu = u - mu - c
    sum(lnNpr.(nl,nu) .+ .5*mu.^2 .- x.*mu)
end

# implements gradient of psi(x) to find optimal exponential twisting
# assumes scaled 'L' with zero diagonal
function gradpsi(y, L, l, u)
    d = length(u)
    c = zeros(d,1)
    x, mu = deepcopy(c), deepcopy(c)
    x[1:(d-1)] = y[1:(d-1)]
    mu[1:(d-1)] = y[d:end]

    # compute now ~l and ~u
    c[2:d] = L[2:d,:]*x
    lt = l - mu - c
    ut = u - mu - c

    # compute gradients avoiding catastrophic cancellation
    w = lnNpr.(lt, ut)
    pl = exp.(-0.5*lt.^2 .- w)./sqrt(2*pi)
    pu = exp.(-0.5*ut.^2 .- w)./sqrt(2*pi)
    P = pl - pu

    # output the gradient
    dfdx = -mu[1:(d-1)] .+ (P'*L[:,1:(d-1)])'
    dfdm = mu - x + P
    grad = [dfdx; dfdm[1:(d-1)]]

    # Compute the jacbian
    lt[isinf.(lt)] .= 0
    ut[isinf.(ut)] .= 0
    dP = -P.^2 + lt.*pl .- ut.*pu # dPdm
    DL = repeat(dP, 1, d).*L
    mx = -Matrix{Float64}(I, d, d) + DL
    xx = L'*DL
    mx = mx[1:d-1,1:d-1]

    xx = xx[1:d-1,1:d-1]

    J = [xx mx'; mx diagm(0 => (1 .+ dP[1:d-1]))]

    grad, J
end

#  Computes permuted lower Cholesky factor L for Σ
#  by permuting integration limit vectors l and u.
#  Outputs perm, such that Σ(perm,perm)=L*L'.
#
# Reference:
#  Gibson G. J., Glasbey C. A., Elston D. A. (1994),
#  "Monte Carlo evaluation of multivariate normal integrals and
#  sensitivity to variate ordering",
#  In: Advances in Numerical Methods and Applications, pages 120--126
function cholperm(Σ_in, l, u)
    Σ = copy(Σ_in)
    d = length(l)
    perm, L, z = [1:d...], zeros(d,d), zeros(d,1)
    for j = 1:d
        pr = Inf*ones(d) # compute marginal prob.
        rd = collect(j:d) # search remaining dimensions
        ud = collect(1:j-1)
        D = diag(Σ)
        s = D[rd] .- sum(L[rd,ud].^2, dims=2)
        s[s.<0] .= 1e-16
        s = sqrt.(s)

        L_z = L[rd,ud]*z[ud]

        tl = (l[rd] .- L_z)./s
        tu = (u[rd] .- L_z)./s

        pr[rd] = lnNpr.(tl,tu)

        # find smallest marginal dimension
        k = argmin(pr)

        # flip dimensions k-->j
        jk, kj = [j,k], [k,j]

        # update rows and cols of Σ
        Σ[jk, :] .= Σ[kj, :]
        Σ[:, jk] .= Σ[:, kj]

        # update only rows of L
        L[jk,:] .= L[kj,:]

        # update integration limits
        l[jk] .= l[kj]
        u[jk] .= u[kj]
        perm[jk] .= perm[kj] # keep track of permutation

        # construct L sequentially via Cholesky computation
        s = Σ[j,j] - sum(L[j,ud].^2)

        if s<-0.01 error("Σ is not positive semi-definite") end
        if s < 0 s = 1e-16 end
        L[j,j] = sqrt(s)
        Ld = L[j,j]

        L[j+1:d,j:j] .= (Σ[j+1:d,j]-L[j+1:d,1:j-1]*L[j:j,1:j-1]')./Ld

        # find mean value, z(j), of truncated normal:
        L_z = L[j:j,ud]*z[ud]

        tl = (l[j].-L_z)./Ld
        tu = (u[j].-L_z)./Ld
        w = lnNpr.(tl,tu) # aids in computing expected value of trunc. normal

        @assert length(tl) == 1 && length(tu) == 1 && length(w) == 1
        z[j] = (exp(-.5*tl[1]^2 - w[1]) - exp(-5*tu[1]^2. - w[1]))/sqrt2pi
    end
    L, l, u, perm
end

# truncated multivariate normal generator
# simulates 'n' random vectors exactly/perfectly distributed
# from the d-dimensional N(0,Sig) distribution (zero-mean normal
# with covariance 'Sig') conditional on l<X<u;
# infinite values for 'l' and 'u' are accepted;
# output:   'd' times 'n' array 'rv' storing random vectors;
#
# * Example:
#  d=60;n=10^3;Sig=0.9*ones(d,d)+0.1*eye(d);l=(1:d)/d*4;u=l+2;
#  X=mvrandn(l,u,Sig,n);boxplot(X','plotstyle','compact')
# plot marginals
#
# * Notes: Algorithm may not work if 'Sig' is close to being rank deficient.
# See also: mvNcdf, mvNqmc, mvrorth
# For more help, see Truncated Multivariate Student & Normal documentation at the bottom.
function mvrandn(rng::AbstractRNG, l, u, Σ, n)
    # basic input check
    d = length(l)
    if  length(u) != d || size(Σ) != (d,d) || any(l > u)
        error("l, u, and Σ have to match in dimension with u>l")
    end

    # Cholesky decomposition of matrix with permuation
    Lfull, l, u, perm = cholperm(Σ,l,u) # outputs the permutation
    D = diag(Lfull)
    if any(D .< 1e-16)
        warning("Method may fail as covariance matrix is singular!")
    end
    L = Lfull ./ repeat(D,1,d)

    u = u ./ D
    l = l ./ D #rescale
    L = L .- Matrix{Float64}(I, d, d) # remove diagonal

    # find optimal tilting parameter non-linear equation solver
    soln = nlsolve(only_fj((x) -> gradpsi(x, L, l, u)), zeros(2*(d-1),1))

    x = [soln.zero[1:(d-1)]..., 0]
    mu = [soln.zero[d:(2*d-2)]..., 0]

    # compute psi star
    psistar = psy(x, L, l, u, mu)

    # start acceptance rejection sampling
    rv, accept, iter = zeros(d, 0), 0, 0

    while accept < n #  while # of accepted is less than n
        logpr, Z = mvnrnd(rng, n, L, l, u, mu) # simulate n proposals
        idx = -log.(rand(n)) .> (psistar .- logpr) #acceptance tests
        rv = hcat(rv, Z[:, idx]) #accumulate accepted
        accept = size(rv,2) #keep track of # of accepted
        iter = iter+1 #keep track of while loop iterations
        if iter > 1e3 # if iterations too large, seek approximation only
            accept = n
            rv = hcat(rv, Z) # add the approximate samples
            # println("WARNING: Sample is only approximately distributed.")
        end
    end
    # finish sampling postprocessing
    order = sortperm(perm)

    rv = rv[:, 1:n] # cut-down the array to desired n samples
    rv = Lfull*rv # reverse scaling of L
    rv = rv[order,:] # reverse the Cholesky permutation
end

# Samples a constrained mvrandn with a nonzero mean
function mvrandn_μ(rng::AbstractRNG, μ, l, u, Σ, n)
    Y = mvrandn(rng, l - μ, u - μ, Σ, n)
    return Y .+ μ
end


text:
# Suppress warnings during testing.
setlevel!(getlogger(PowerModels), "error")

# is this the best way to find a file in a package?
base_dir = dirname(dirname(pathof(PowerSystems)))
case5_data = PM.parse_file(joinpath(base_dir,"data/matpower/case5.m"))
case5_data = PM.replicate(case5_data, 2)

case5_dc_data = PM.parse_file(joinpath(base_dir,"data/matpower/case5_dc.m"))
case5_dc_data = PM.replicate(case5_dc_data, 2)

# TODO: currently JuMP.num_variables is the best we can do to introspect the JuMP model.
#  Ideally this would also test the number of constraints generated

@testset "PowerModels Model Build" begin
    pm = PowerSimulations.build_nip_model(case5_data, PM.DCPPowerModel)
    @test JuMP.num_variables(pm.model) == 34
    pm = PowerSimulations.build_nip_model(case5_data, PM.ACPPowerModel)
    @test JuMP.num_variables(pm.model) == 96
    pm = PowerSimulations.build_nip_model(case5_data, PM.SOCWRPowerModel)
    @test JuMP.num_variables(pm.model) == 110
end

# test PowerSimulations type extentions
DCAngleModel = (data::Dict{String,Any}; kwargs...) -> PM.GenericPowerModel(data, PM.DCPlosslessForm; kwargs...)

StandardACModel = (data::Dict{String,Any}; kwargs...) -> PM.GenericPowerModel(data, PM.StandardACPForm; kwargs...)

@testset "PM with type extensions" begin
    pm = PowerSimulations.build_nip_model(case5_data, DCAngleModel)
    JuMP.num_variables(pm.model) == 34
    pm = PowerSimulations.build_nip_model(case5_data, StandardACModel)
    JuMP.num_variables(pm.model) == 96
    pm = PowerSimulations.build_nip_model(case5_dc_data, PM.DCPPowerModel)
    JuMP.num_variables(pm.model) == 36
    pm = PowerSimulations.build_nip_model(case5_dc_data, DCAngleModel)
    JuMP.num_variables(pm.model) == 48
end

@testset "PM integration into PS" begin
    base_dir = dirname(dirname(pathof(PowerSystems)))
    include(joinpath(base_dir,"data/data_5bus_pu.jl"))
    PS_struct = PSY.System(nodes5, generators5, loads5_DA, branches5, nothing,  100.0);
    PM_dict = PowerSimulations.pass_to_pm(PS_struct)
    PM_object = PowerSimulations.build_nip_model(PM_dict, DCAngleModel);
    @test JuMP.num_variables(PM_object.model) == 384
end

ext: jl


text:
abstract type LogLevel end

abstract type None <: LogLevel end
abstract type Value <: LogLevel end
abstract type Decision <: LogLevel end
abstract type Gradient <: LogLevel end
abstract type Full <: LogLevel end

include("simplelogger.jl")
include("progresslogger.jl")
include("historylogger.jl")

History() = HistoryLogger.History()
History(::Type{L}) where L <: LogLevel = HistoryLogger.History(L)
History(logger::AbstractLogger) = HistoryLogger.History(logger)


text:
import JSON

# JSON.jl doesn't define a serialization method for MIME objects, so we add one ourselves:
import JSON: lower
JSON.lower(m::MIME) = string(m)
JSON.lower(u::UUID) = string(u)

function serialize_message_to_stream(io::IO, message::UpdateMessage)
    to_send = Dict(:type => message.type, :message => message.message)
    if message.notebook !== nothing
        to_send[:notebookID] = message.notebook.uuid
    end
    if message.cell !== nothing
        to_send[:cellID] = message.cell.uuid
    end
    if message.initiator !== missing
        to_send[:initiatorID] = message.initiator.clientID
        to_send[:requestID] = message.initiator.requestID
    end

    JSON.print(io, to_send)
end

function serialize_message(message::UpdateMessage)
    sprint(serialize_message_to_stream, message)
end

"To be used in `make_paths_distinct!`"
mutable struct NotebookPath
    uuid
    path_split
    current_path
    current_depth
end

function count_occurances(vals)
    counts = Dict()
    for v in vals
        old = get(counts, v, 0)
        counts[v] = old + 1
    end
    counts
end

"""For internal use. Takes a Set of `NotebookPath`s and gives each a short path (e.g. `to/file.jl` from `/path/to/file.jl`), with the guarantee that all final short paths will be distinct.

For example, the set

`/a/b/c.jl`, `/a/P/c.jl`, `/Q/b/c.jl`, '/a/b/R.jl'

will become

`/a/b/c.jl`, `P/c.jl`, `/Q/b/c.jl`, 'R.jl'"""
function make_paths_distinct!(notebookpaths::Set{NotebookPath})
    counts = count_occurances(np.current_path for np in notebookpaths)
    for (current_path, count) in counts
        if count == 1 && !isempty(current_path)
            # done!
        else
            # these need to be made distinct by extending their paths

            not_yet_distinct = filter(notebookpaths) do np
                np.current_path == current_path
            end

            for np in not_yet_distinct
                np.current_depth += 1
                np.current_path = joinpath(np.path_split[end-np.current_depth : end]...)
                if np.current_depth == length(np.path_split) - 1
                    delete!(not_yet_distinct, np)
                    if !Sys.iswindows()
                        np.current_path = '/' * np.current_path
                    end
                end
            end

            make_paths_distinct!(not_yet_distinct)
        end
    end
end

function handle_changecell(notebook, cell, newcode; initiator::Union{Initiator, Missing}=missing)::Task
    # i.e. Ctrl+Enter was pressed on this cell
    # we update our `Notebook` and start execution

    # don't reparse when code is identical (?)
    if cell.code != newcode
        cell.code = newcode
        cell.parsedcode = nothing
    end

    # TODO: feedback to user about File IO
    save_notebook(notebook)

    putnotebookupdates!(notebook, clientupdate_cell_input(notebook, cell, initiator=initiator))

    run_reactive_async!(notebook, cell)
end

responses[:connect] = (body, notebook=nothing; initiator::Union{Initiator, Missing}=missing) -> begin
    putclientupdates!(initiator, UpdateMessage(:👋, Dict(:notebookExists => (notebook != nothing)), nothing, nothing, initiator))
end

responses[:getversion] = (body, notebook=nothing; initiator::Union{Initiator, Missing}=missing) -> begin
    putclientupdates!(initiator, UpdateMessage(:versioninfo, Dict(
        :pluto => PLUTO_VERSION_STR,
        :julia => JULIA_VERSION_STR,
    ), nothing, nothing, initiator))
end

# TODO: actions on the notebook are not thread safe
responses[:addcell] = (body, notebook::Notebook; initiator::Union{Initiator, Missing}=missing) -> begin
    new_index = body["index"] + 1 # 0-based index (js) to 1-based index (julia)

    new_cell = Cell("")

    insert!(notebook.cells, new_index, new_cell)

    putnotebookupdates!(notebook, clientupdate_cell_added(notebook, new_cell, new_index, initiator=initiator))
end

responses[:deletecell] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator, Missing}=missing) -> begin
    to_delete = cell

    # replace the cell's code with "" and do a reactive run
    runtask = handle_changecell(notebook, to_delete, "", initiator=initiator)

    # wait for the reactive run to finish, then delete the cells
    # we wait async, to make sure that the web server remains responsive
    @async begin
        wait(runtask)

        filter!(c->c.uuid ≠ to_delete.uuid, notebook.cells)
        putnotebookupdates!(notebook, clientupdate_cell_deleted(notebook, to_delete, initiator=initiator))
        save_notebook(notebook) # this might be "too late", but it will save the latest version of `notebook` anyways
    end
end

responses[:movecell] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator, Missing}=missing) -> begin
    to_move = cell

    # Indexing works as if a new cell is added.
    # e.g. if the third cell (at julia-index 3) of [0, 1, 2, 3, 4]
    # is moved to the end, that would be new julia-index 6

    new_index = body["index"] + 1 # 0-based index (js) to 1-based index (julia)
    old_index = findfirst(isequal(to_move), notebook.cells)

    # Because our cells run in _topological_ order, we don't need to reevaluate anything.
    if new_index < old_index
        deleteat!(notebook.cells, old_index)
        insert!(notebook.cells, new_index, to_move)
    elseif new_index > old_index + 1
        insert!(notebook.cells, new_index, to_move)
        deleteat!(notebook.cells, old_index)
    end

    putnotebookupdates!(notebook, clientupdate_cell_moved(notebook, to_move, new_index, initiator=initiator))
end

responses[:changecell] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator, Missing}=missing) -> begin
    newcode = body["code"]

    handle_changecell(notebook, cell, newcode, initiator=initiator)
end

responses[:runall] = (body, notebook::Notebook; initiator::Union{Initiator, Missing}=missing) -> begin
    to_update = run_reactive_async!(notebook, notebook.cells)
end

responses[:getinput] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator, Missing}=missing) -> begin
    putclientupdates!(initiator, clientupdate_cell_input(notebook, cell, initiator=initiator))
end

responses[:getoutput] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator, Missing}=missing) -> begin
    putclientupdates!(initiator, clientupdate_cell_output(notebook, cell, initiator=initiator))
end

responses[:getallcells] = (body, notebook::Notebook; initiator::Union{Initiator, Missing}=missing) -> begin
    # TODO: the client's update channel might get full
    update = UpdateMessage(:cell_list,
        Dict(:cells => [Dict(
                :uuid => string(cell.uuid),
                ) for cell in notebook.cells]), nothing, nothing, initiator)

    putclientupdates!(initiator, update)
end

responses[:getallnotebooks] = (body, notebook=nothing; initiator::Union{Initiator, Missing}=missing) -> begin
    putplutoupdates!(clientupdate_notebook_list(notebooks, initiator=initiator))
end

responses[:movenotebookfile] = (body, notebook::Notebook; initiator::Union{Initiator, Missing}=missing) -> begin
    newpath = body["path"]
    result = try
        if isfile(newpath)
            (success=false,message="File already exists")
        else
            move_notebook(notebook, newpath)
            putplutoupdates!(clientupdate_notebook_list(notebooks, initiator=initiator))
            (success=true, message="")
        end
    catch ex
        showerror(stderr, stacktrace(backtrace()))
        (success=false, message=sprint(showerror, ex))
    end

    update = UpdateMessage(:move_notebook_result, result)
    putclientupdates!(initiator, update)
end

responses[:interruptall] = (body, notebook::Notebook; initiator::Union{Initiator, Missing}=missing) -> begin
    success = WorkspaceManager.kill_workspace(notebook)
    # TODO: notify user whether interrupt was successful (i.e. whether they are using a `ProcessWorkspace`)
end


text:
# Function that reduces the forces vector and energy after parallel computation

function reduceForces!(f :: Forces)
  u = 0.
  for i in 1:Threads.nthreads()
    u = u + f.upartial[i]
    @. f.f = f.f + f.fpartial[i]
  end
  return u
end


text:
export tree_greedy, MinSpaceOut, MinSpaceDiff

struct MinSpaceOut end
struct MinSpaceDiff end

struct LegInfo{ET}
    l1::Vector{ET}
    l2::Vector{ET}
    l12::Vector{ET}
    l01::Vector{ET}
    l02::Vector{ET}
    l012::Vector{ET}
end

"""
    tree_greedy(incidence_list, log2_sizes; method=MinSpaceOut())

Compute greedy order, and the time and space complexities, the rows of the `incidence_list` are vertices and columns are edges.
`log2_sizes` are defined on edges.

```julia
julia> code = ein"(abc,cde),(ce,sf,j),ak->ael"
aec, ec, ak -> ael
├─ ce, sf, j -> ec
│  ├─ sf
│  ├─ j
│  └─ ce
├─ ak
└─ abc, cde -> aec
   ├─ cde
   └─ abc

julia> optimize_greedy(code, Dict([c=>2 for c in "abcdefjkls"]))
ae, ak -> ea
├─ ak
└─ aec, ec -> ae
   ├─ ce,  -> ce
   │  ├─ sf, j ->
   │  │  ├─ j
   │  │  └─ sf
   │  └─ ce
   └─ abc, cde -> aec
      ├─ cde
      └─ abc
```
"""
function tree_greedy(incidence_list::IncidenceList{VT,ET}, log2_edge_sizes; method=MinSpaceOut(), nrepeat=10) where {VT,ET}
    @assert nrepeat >= 1
    best_tree, best_tcs, best_scs = _tree_greedy(incidence_list, log2_edge_sizes; method=method)
    best_tc, best_sc = log2sumexp2(best_tcs), maximum(best_scs)
    for _ = 1:nrepeat-1
        tree, tcs, scs = _tree_greedy(incidence_list, log2_edge_sizes; method=method)
        tc, sc = log2sumexp2(tcs), maximum(scs)
        if sc < best_sc || (sc <= best_sc && tc < best_tc)
            best_tcs, best_scs, best_tc, best_sc, best_tree = tcs, scs, tc, sc, tree
        end
    end
    return best_tree, best_tcs, best_scs
end

function _tree_greedy(incidence_list::IncidenceList{VT,ET}, log2_edge_sizes; method=MinSpaceOut()) where {VT,ET}
    incidence_list = copy(incidence_list)
    n = nv(incidence_list)
    if n == 0
        return nothing
    elseif n == 1
        return collect(vertices(incidence_list))[1]
    end
    log2_tcs = Float64[] # time complexity
    log2_scs = Float64[]

    tree = Dict{VT,Any}([v=>v for v in vertices(incidence_list)])
    cost_values = evaluate_costs(method, incidence_list, log2_edge_sizes)
    while true
        if length(cost_values) == 0
            vpool = collect(vertices(incidence_list))
            pair = minmax(vpool[1], vpool[2])  # to prevent empty intersect
        else
            pair = find_best_cost(cost_values)
        end
        log2_tc_step, sc, code = contract_pair!(incidence_list, pair..., log2_edge_sizes)
        push!(log2_tcs, log2_tc_step)
        push!(log2_scs, space_complexity(incidence_list, log2_edge_sizes))
        if nv(incidence_list) > 1
            tree[pair[1]] = ContractionTree(tree[pair[1]], tree[pair[2]])
        else
            return ContractionTree(tree[pair[1]], tree[pair[2]]), log2_tcs, log2_scs
        end
        update_costs!(cost_values, pair..., method, incidence_list, log2_edge_sizes)
    end
end

function contract_pair!(incidence_list, vi, vj, log2_edge_sizes)
    log2dim(legs) = isempty(legs) ? 0 : sum(l->log2_edge_sizes[l], legs)  # for 1.5, you need this patch because `init` kw is not allowed.
    @assert vj > vi
    # compute time complexity and output tensor
    legsets = analyze_contraction(incidence_list, vi, vj)
    D12,D01,D02,D012 = log2dim.(getfield.(Ref(legsets),3:6))
    tc = D12+D01+D02+D012  # dangling legs D1 and D2 do not contribute

    # einsum code
    eout = legsets.l01 ∪ legsets.l02 ∪ legsets.l012
    code = (edges(incidence_list, vi), edges(incidence_list, vj)) => eout
    sc = log2dim(eout)

    # change incidence_list
    delete_vertex!(incidence_list, vj)
    change_edges!(incidence_list, vi, eout)
    for e in eout
        replace_vertex!(incidence_list, e, vj=>vi)
    end
    remove_edges!(incidence_list, legsets.l1 ∪ legsets.l2 ∪ legsets.l12)
    return tc, sc, code
end

function evaluate_costs(method, incidence_list::IncidenceList{VT,ET}, log2_edge_sizes) where {VT,ET}
    # initialize cost values
    cost_values = Dict{Tuple{VT,VT},Float64}()
    for vi = vertices(incidence_list)
        for vj in neighbors(incidence_list, vi)
            if vj > vi
                cost_values[(vi,vj)] = greedy_loss(method, incidence_list, log2_edge_sizes, vi, vj)
            end
        end
    end
    return cost_values
end

function update_costs!(cost_values, va, vb, method, incidence_list::IncidenceList{VT,ET}, log2_edge_sizes) where {VT,ET}
    for vj in neighbors(incidence_list, va)
        vx, vy = minmax(vj, va)
        cost_values[(vx,vy)] = greedy_loss(method, incidence_list, log2_edge_sizes, vx, vy)
    end
    for k in keys(cost_values)
        if vb ∈ k
            delete!(cost_values, k)
        end
    end
end

function find_best_cost(cost_values::Dict{PT}) where PT
    length(cost_values) < 1 && error("cost value information missing")
    minval = minimum(Base.values(cost_values))
    pairs = PT[]
    for (k, v) in cost_values
        if v == minval
            push!(pairs, k)
        end
    end
    return rand(pairs)
end

function analyze_contraction(incidence_list::IncidenceList{VT,ET}, vi::VT, vj::VT) where {VT,ET}
    ei = edges(incidence_list, vi)
    ej = edges(incidence_list, vj)
    leg012,leg12,leg1,leg2,leg01,leg02 = ET[], ET[], ET[], ET[], ET[], ET[]
    # external legs
    for leg in ei ∪ ej
        isext = leg ∈ incidence_list.openedges || !all(x->x==vi || x==vj, vertices(incidence_list, leg))
        if isext
            if leg ∈ ei
                if leg ∈ ej
                    push!(leg012, leg)
                else
                    push!(leg01, leg)
                end
            else
                push!(leg02, leg)
            end
        else
            if leg ∈ ei
                if leg ∈ ej
                    push!(leg12, leg)
                else
                    push!(leg1, leg)
                end
            else
                push!(leg2, leg)
            end
        end
    end
    return LegInfo(leg1, leg2, leg12, leg01, leg02, leg012)
end

function greedy_loss(::MinSpaceOut, incidence_list, log2_edge_sizes, vi, vj)
    log2dim(legs) = isempty(legs) ? 0 : sum(l->log2_edge_sizes[l], legs)  # for 1.5, you need this patch because `init` kw is not allowed.
    legs = analyze_contraction(incidence_list, vi, vj)
    log2dim(legs.l01)+log2dim(legs.l02)+log2dim(legs.l012)
end

function greedy_loss(::MinSpaceDiff, incidence_list, log2_edge_sizes, vi, vj)
    log2dim(legs) = isempty(legs) ? 0 : sum(l->log2_edge_sizes[l], legs)  # for 1.5, you need this patch because `init` kw is not allowed.
    legs = analyze_contraction(incidence_list, vi, vj)
    D1,D2,D12,D01,D02,D012 = log2dim.(getfield.(Ref(legs), 1:6))
    exp2(D01+D02+D012) - exp2(D01+D12+D012) - exp2(D02+D12+D012)  # out - in
end

function space_complexity(incidence_list, log2_sizes)
    sc = 0.0
    for v in vertices(incidence_list)
        for e in edges(incidence_list, v)
            sc += log2_sizes[e]
        end
    end
    return sc
end

function contract_tree!(incidence_list::IncidenceList, tree::ContractionTree, log2_edge_sizes, tcs, scs)
    vi = tree.left isa ContractionTree ? contract_tree!(incidence_list, tree.left, log2_edge_sizes, tcs, scs) : tree.left
    vj = tree.right isa ContractionTree ? contract_tree!(incidence_list, tree.right, log2_edge_sizes, tcs, scs) : tree.right
    tc, sc, code = contract_pair!(incidence_list, vi, vj, log2_edge_sizes)
    push!(tcs, tc)
    push!(scs, sc)
    return vi
end


text:
@testset "Minimal models of elliptic curves" begin

  @testset "Minimal model (Laska-Kraus-Connell" begin
    E = EllipticCurve([1,2,3,4,5])
    EE = @inferred laska_kraus_connell(E)
    @test a_invars(EE) == (1, -1, 0, 4, 3)

    E = EllipticCurve([625, -15625, 19531250, -2929687500, -34332275390625])
    EE = @inferred laska_kraus_connell(E)
    @test a_invars(EE) == (1, -1, 0, 4, 3)
  end

  @testset "Tates algorithm" begin
    E = EllipticCurve([625, -15625, 19531250, -2929687500, -34332275390625])
    EE = @inferred tates_algorithm_global(E)
    @test a_invars(EE) == (1, -1, 0, 4, 3)

    E = EllipticCurve([1,2,3,4,5])
    EE = @inferred tates_algorithm_global(E)
    @test a_invars(EE) == (1, -1, 0, 4, 3)

    #  25350.a1
    E = EllipticCurve([1, 1, 0, 40050, 7557750])
    Ep, K, f, c = tates_algorithm_local(E, 2)
    @test a_invars(tidy_model(Ep)) == a_invars(E)
    @test K == "I1"
    @test f == 1
    @test c == 1

    Ep, K, f, c = tates_algorithm_local(E, 3)
    @test a_invars(tidy_model(Ep)) == a_invars(E)
    @test K == "I2"
    @test f == 1
    @test c == 2

    Ep, K, f, c = tates_algorithm_local(E, 5)
    @test a_invars(tidy_model(Ep)) == a_invars(E)
    @test K == "III*"
    @test f == 2
    @test c == 2

    Ep, K, f, c = tates_algorithm_local(E, 13)
    @test a_invars(tidy_model(Ep)) == a_invars(E)
    @test K == "IV*"
    @test f == 2
    @test c == 1

    # 150.a1
    E = EllipticCurve([1, 1, 0, -20700, 1134000])
    Ep, K, f, c = tates_algorithm_local(E, 2)
    @test a_invars(tidy_model(Ep)) == a_invars(E)
    @test K == "I5"
    @test f == 1
    @test c == 1

    Ep, K, f, c = tates_algorithm_local(E, 3)
    @test a_invars(tidy_model(Ep)) == a_invars(E)
    @test K == "I10"
    @test f == 1
    @test c == 2

    Ep, K, f, c = tates_algorithm_local(E, 5)
    @test a_invars(tidy_model(Ep)) == a_invars(E)
    @test K == "III*"
    @test f == 2
    @test c == 2
  end

  @testset "Tates algorithm over number fields" begin
    #100.1-b2
    Rx, x = PolynomialRing(QQ, "x")
    L, a = number_field(x^2-x-1)
    E = EllipticCurve(L, [1, 1, 1, -3, 1])
    F, phi = transform_rstu(E,[12, -1, 1+2*a, 3+a])
    F, phi = integral_model(F)
    OL = ring_of_integers(L)
    P = 2*OL
    Ep, K, f, c, s = tates_algorithm_local(F, P)
    @test K == "I5"
    @test f == 1
    @test c == 5
    @test s == true
    @test valuation(discriminant(E),P) == valuation(discriminant(Ep),P)

    P = numerator(ideal(OL, -2*a+1))

    Ep, K, f, c, s = tates_algorithm_local(Ep, P)
    @test K == "IV"
    @test f == 2
    @test c == 3
    @test s == true

    @test valuation(discriminant(E),P) == valuation(discriminant(Ep),P)

    #49.1-CMa1
    L, a = number_field(x^2-x+1)
    E = EllipticCurve(L, [0, 1+a, a, a, 0])
    F, phi = transform_rstu(E,[4, -1+a, 7, a-4])
    F, phi = integral_model(F)
    OL = ring_of_integers(L)

    P = numerator(ideal(OL, -3*a+1))

    Ep, K, f, c, s = tates_algorithm_local(F, P)
    @test K == "II"
    @test f == 2
    @test c == 1
    @test s == true

    @test valuation(discriminant(E),P) == valuation(discriminant(Ep),P)
    #673.1-a1
    E = EllipticCurve(L, [0, a, a, -1+a, 0])
    F, phi = transform_rstu(E,[a, 0, -3+a, 7])
    F, phi = integral_model(F)
    P = numerator(ideal(OL, 29*a-8))
    Ep, K, f, c, s = tates_algorithm_local(F, P)
    @test K == "I1"
    @test f == 1
    @test c == 1
    @test s == false

    @test valuation(discriminant(E),P) == valuation(discriminant(Ep),P)
    #121.1-c3
    L, a = number_field(x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1)
    E = EllipticCurve(L, [0, a-1, a^3+a^2-2*a-1, -2*a-1, -a^4-a^3+a^2-a-2])
    F, phi = transform_rstu(E,[5, -1+a^2, a^3, 2*a-a^4])
    F, phi = integral_model(F)
    OL = ring_of_integers(L)
    P = numerator(ideal(OL, a^2+a-2))
    Ep, K, f, c, s = tates_algorithm_local(F, P)
    @test K == "I5*"
    @test f == 2
    @test c == 2
    @test s == true

    @test valuation(discriminant(E),P) == valuation(discriminant(Ep),P)
  end
end


text:
using HTTP
using CSV
using DataFrames
using GLMakie
using LinearAlgebra
GLMakie.inline!(false)
url = "https://raw.githubusercontent.com/plotly/datasets/master/vortex.csv"
df = CSV.File(HTTP.get(url).body)|> DataFrame;

uvw = Vec3f.(df.u, df.v, df.w)
len = norm.(uvw)

arrows(df.x, df.y, df.z, df.u, df.v, df.w, linewidth=0, color=len, arrowsize = Vec3f(0.5, 0.5, 1)) |> display

# surface parameterization  (u, v)----> (f(u,v), g(u,v), h(u,v)):
f = (u, v) -> cos(v)*(6 - (5/4 +sin(3*u))*sin(u-3*v))
g = (u, v) -> sin(v)*(6 - (5/4 +sin(3*u))*sin(u-3*v))
h = (u, v) -> -cos(u-3*v) * (5/4 +sin(3*u));

u = range(0, stop=2π, length=150)
v = range(0, stop=2π, length=150)

tr = surface(f.(u,v'),
             g.(u,v'),
             h.(u,v'),
             ambient=Vec3f(0.5),
             diffuse=Vec3f(1),
             specular=0.5,
             colormap=:balance) |> display


text:
using Random, Distributions, Plots; pyplot()
Random.seed!(1)

Z() = sqrt(-2*log(rand()))*cos(2*pi*rand())
xGrid = -4:0.01:4

histogram([Z() for _ in 1:10^6], bins=50,
		normed=true, label="MC estimate")
plot!(xGrid, pdf.(Normal(),xGrid),
	 c=:red, lw=4, label="PDF",
	 xlims=(-4,4), ylims=(0,0.5), xlabel="x", ylabel="f(x)")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
export MADDPGManager

"""
    MADDPGManager(; agents::Dict{<:Any, <:Agent}, args...)
Multi-agent Deep Deterministic Policy Gradient(MADDPG) implemented in Julia. By default, `MADDPGManager` uses for simultaneous environments with [continuous action space](https://spinningup.openai.com/en/latest/spinningup/rl_intro.html#stochastic-policies).
See the paper https://arxiv.org/abs/1706.02275 for more details.

# Keyword arguments
- `agents::Dict{<:Any, <:Agent}`, here each agent collects its own information. While updating the policy, each **critic** will assemble all agents'
  trajectory to update its own network. **Note that** here the policy of the `Agent` should be `DDPGPolicy` wrapped by `NamedPolicy`, see the relative
  experiment([`MADDPG_KuhnPoker`](https://juliareinforcementlearning.org/docs/experiments/experiments/Policy%20Gradient/JuliaRL_MADDPG_KuhnPoker/#JuliaRL\\_MADDPG\\_KuhnPoker) or [`MADDPG_SpeakerListener`](https://juliareinforcementlearning.org/docs/experiments/experiments/Policy%20Gradient/JuliaRL_MADDPG_SpeakerListener/#JuliaRL\\_MADDPG\\_SpeakerListener)) for references.
- `traces`, set to `SARTS` if you are apply to an environment of `MINIMAL_ACTION_SET`, or `SLARTSL` if you are to apply to an environment of `FULL_ACTION_SET`.
- `batch_size::Int`
- `update_freq::Int`
- `update_step::Int`, count the step.
- `rng::AbstractRNG`.
"""
mutable struct MADDPGManager <: AbstractPolicy
    agents::Dict{<:Any, <:Agent}
    traces
    batch_size::Int
    update_freq::Int
    update_step::Int
    rng::AbstractRNG
end

# used for simultaneous environments.
function (π::MADDPGManager)(env::AbstractEnv)
    while current_player(env) == chance_player(env)
        env |> legal_action_space |> rand |> env
    end
    Dict(
        player => agent.policy(env)
    for (player, agent) in π.agents)
end

function (π::MADDPGManager)(stage::Union{PreEpisodeStage, PostActStage}, env::AbstractEnv)
    # only need to update trajectory.
    for (_, agent) in π.agents
        update!(agent.trajectory, agent.policy, env, stage)
    end
end

function (π::MADDPGManager)(stage::PreActStage, env::AbstractEnv, actions)
    # update each agent's trajectory.
    for (player, agent) in π.agents
        update!(agent.trajectory, agent.policy, env, stage, actions[player])
    end

    # update policy
    update!(π, env)
end

function (π::MADDPGManager)(stage::PostEpisodeStage, env::AbstractEnv)
    # collect state and a dummy action to each agent's trajectory here.
    for (_, agent) in π.agents
        update!(agent.trajectory, agent.policy, env, stage)
    end

    # update policy
    update!(π, env)
end

# update policy
function RLBase.update!(π::MADDPGManager, env::AbstractEnv)
    π.update_step += 1
    π.update_step % π.update_freq == 0 || return

    for (_, agent) in π.agents
        length(agent.trajectory) > agent.policy.policy.update_after || return
        length(agent.trajectory) > π.batch_size || return
    end

    # get training data
    temp_player = collect(keys(π.agents))[1]
    t = π.agents[temp_player].trajectory
    inds = rand(π.rng, 1:length(t), π.batch_size)
    batches = Dict((player, RLCore.fetch!(BatchSampler{π.traces}(π.batch_size), agent.trajectory, inds))
                for (player, agent) in π.agents)

    # get s, a, s′ for critic
    s = vcat((batches[player][:state] for (player, _) in π.agents)...)
    a = vcat((batches[player][:action] for (player, _) in π.agents)...)
    s′ = vcat((batches[player][:next_state] for (player, _) in π.agents)...)

    for (player, agent) in π.agents
        p = agent.policy.policy # get agent's concrete DDPGPolicy.

        A = p.behavior_actor
        C = p.behavior_critic
        Aₜ = p.target_actor
        Cₜ = p.target_critic

        γ = p.γ
        ρ = p.ρ

        # by default A, C, Aₜ, Cₜ on the same device.
        batches, s, a, s′ = send_to_device(device(A), (batches, s, a, s′))
        r = batches[player][:reward]
        t = batches[player][:terminal]
        # for training behavior_actor.
        mu_actions = vcat(
            ((
                batches[p][:next_state] |>
                a.policy.policy.behavior_actor
            ) for (p, a) in π.agents)...
        )
        # for training behavior_critic.
        new_actions = vcat(
            ((
                batches[p][:next_state] |>
                a.policy.policy.target_actor
            ) for (p, a) in π.agents)...
        )

        if π.traces == SLARTSL
            # Note that by default **MADDPG** is used for the environments with continuous action space, and `legal_action_space_mask` is
            # defined in the environments with discrete action space. So we need `env.action_mapping` to transform the actions
            # getting from the trajectory.
            @assert env isa ActionTransformedEnv

            mask = batches[player][:next_legal_actions_mask]
            mu_l′ = Flux.batch(
                (begin
                    actions = env.action_mapping(mu_actions[:, i])
                    mask[actions[player]]
                end for i = 1:π.batch_size)
            )
            new_l′ = Flux.batch(
                (begin
                    actions = env.action_mapping(new_actions[:, i])
                    mask[actions[player]]
                end for i = 1:π.batch_size)
            )
        end

        qₜ = Cₜ(vcat(s′, new_actions)) |> vec
        if π.traces == SLARTSL
            qₜ .+= ifelse.(new_l′, 0.0f0, typemin(Float32))
        end
        y = r .+ γ .* (1 .- t) .* qₜ

        gs1 = gradient(Flux.params(C)) do
            q = C(vcat(s, a)) |> vec
            loss = mean((y .- q) .^ 2)
            ignore() do
                p.critic_loss = loss
            end
            loss
        end

        update!(C, gs1)

        gs2 = gradient(Flux.params(A)) do
            v = C(vcat(s, mu_actions)) |> vec
            if π.traces == SLARTSL
                v .+= ifelse.(mu_l′, 0.0f0, typemin(Float32))
            end
            reg = mean(A(batches[player][:state]) .^ 2)
            loss = -mean(v) +  reg * 1e-3
            ignore() do
                p.actor_loss = loss
            end
            loss
        end

        update!(A, gs2)

        # polyak averaging
        for (dest, src) in zip(Flux.params([Aₜ, Cₜ]), Flux.params([A, C]))
            dest .= ρ .* dest .+ (1 - ρ) .* src
        end
    end
end


text:
@testset "Trainable Potential Default Implmentation Unit Tests" begin
    p = MockTrainablePotential(1.0, 2.0, 3.0)

    @test get_parameters(p) == (; A=1.0, B=2.0)
    @test set_parameters(p, (; A=2.0, B=4.0)) == MockTrainablePotential(2.0, 4.0, 3.0)
    @test serialize_parameters(p) == [1.0, 2.0]
    @test deserialize_parameters(p, [2.0, 4.0]) == MockTrainablePotential(2.0, 4.0, 3.0)

    @test get_hyperparameters(p) == (; C=3.0)
    @test set_hyperparameters(p, (; C=6.0)) == MockTrainablePotential(1.0, 2.0, 6.0)
    @test serialize_hyperparameters(p) == [3.0]
    @test deserialize_hyperparameters(p, [6.0]) == MockTrainablePotential(1.0, 2.0, 6.0)
end


text:
struct LLVMOrdering{name} <: Ordering end

const unordered = LLVMOrdering{:unordered}()
const monotonic = LLVMOrdering{:monotonic}()
const acquire = LLVMOrdering{:acquire}()
const release = LLVMOrdering{:release}()
const acq_rel = LLVMOrdering{:acq_rel}()
const seq_cst = LLVMOrdering{:seq_cst}()

const orderings = (unordered, monotonic, acquire, release, acq_rel, seq_cst)

const ConcreteOrdering = Union{map(typeof, orderings)...}

llvm_ordering(::LLVMOrdering{name}) where {name} = name
Base.string(o::LLVMOrdering) = String(llvm_ordering(o))
Base.print(io::IO, o::LLVMOrdering) = print(io, string(o))

Base.show(io::IO, o::ConcreteOrdering) = print(io, UnsafeAtomics, '.', llvm_ordering(o))


text:
compose(f::Function, g::Function) = (x) -> g(f(x))
@show compose(sin, asin)(0.5)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
using Printf, Statistics, DiscreteEvents, BenchmarkTools
import Dates.now

struct RES
    min::Float64
    mean::Float64
    nevents::Int64
    nsample::Int64
    nactions::Int64
    minact::Float64
    meanact::Float64
end
Base.Array(r::RES) = [r.min, r.mean, r.nevents, r.nsample, r.nactions, r.minact, r.meanact]
function benchresult(trial::BenchmarkTools.Trial, clk::Clock)
    _min = minimum(t.times)*1e-6
    _mean = mean(t.times)*1e-6
    nevents = clk.evcount
    nsample = clk.scount
    nactions = clk.evcount+clk.scount
    minact = _min/nactions*1e3
    meanact = _mean/nactions*1e3
    return RES(_min,_mean,nevents,nsample,nactions,minact,meanact)
end
res = Dict()

btime = now()

println("... channel benchmark ...")
include("bench_channel.jl")
res["channel"] = benchresult(t, 𝐶)

println("... dice benchmark ...")
include("bench_dice.jl")
res["dice"] = benchresult(t, clk)

println("... heating benchmark ...")
include("bench_heating.jl")
res["heating"]  = benchresult(t, 𝐶)

arr = sum(hcat(Array(res["channel"]), Array(res["dice"]), Array(res["heating"])), dims=2)
arr[6] = arr[1]/arr[5]*1e3
arr[7] = arr[2]/arr[5]*1e3
res["sum"] = RES(arr...)

println()
println("Benchmark results:")
println("==================")
@printf "time    \tdatetime               \tchannel [ms] \tdice [ms] \theating [ms] \t    sum [ms]\n"
@printf "minimum \t%s\t%12.3f \t%9.3f \t%12.3f \t%12.3f\n" btime res["channel"].min res["dice"].min res["heating"].min res["sum"].min
@printf "mean    \t%s\t%12.3f \t%9.3f \t%12.3f \t%12.3f\n" btime res["channel"].mean res["dice"].mean res["heating"].mean res["sum"].mean
println()
println("action times:")
println("-------------")
@printf "time    \tdatetime               \tchannel [μs] \tdice [μs] \theating [μs] \toverall [μs]\n"
@printf "minimum \t%s\t%12.3f \t%9.3f \t%12.3f \t%12.3f\n" btime res["channel"].minact res["dice"].minact res["heating"].minact res["sum"].minact
@printf "mean    \t%s\t%12.3f \t%9.3f \t%12.3f \t%12.3f\n" btime res["channel"].meanact res["dice"].meanact res["heating"].meanact res["sum"].meanact


text:
using Documenter, NeXLSpectrum, Gadfly

weavedocs = ( "continuummodel", "errorbars", "K412fit", "K412quick", "XRFSpectra" )

rebuildweave = !all(map(wd->isfile(joinpath("src", "$wd.md")), weavedocs))

if rebuildweave
    map(name->rm(joinpath("src","$name.md")), filter(wd->isfile(joinpath("src","$name.md")),weavedocs))
    let curr=pwd()
        try
            include(joinpath("..","weave","buildweave.jl"))
        finally
            cd(curr)
        end
    end
end

pages = [
    "Home" => "index.md",
    "Spectrum Methods" => "spectrum.md",
    "Fitting K412 (simple API)" => "k412refs.md",
    "Fitting K412 (flexible API)" => "K412fit.md",
    "Fitting K412 (quick fit)" => "K412quick.md",
    # "Quantifying AMD-6005a glass" => "quantAMDglass.md",
    "Fitting XRF Spectra" => "XRFspectra.md",
    "Lovely Error Bars" => "errorbars.md",
    "Modeling the Continuum" => "continuummodel.md",
    "Methods" => "methods.md"
 ]

makedocs(
    modules = [NeXLSpectrum],
    sitename = "NeXLSpectrum.jl",
    pages = pages,
    checkdocs=:exports
)

function addNISTHeaders(htmlfile::String)
    # read HTML
    html = transcode(String,read(htmlfile))
    # Find </head>
    i = findfirst(r"</[Hh][Ee][Aa][Dd]>", html)
    # Already added???
    j = findfirst("nist-header-footer", html)
    if isnothing(j) && (!isnothing(i))
        # Insert nist-pages links right before </head>
        res = html[1:i.start-1]*
            "<link rel=\"stylesheet\" href=\"https://pages.nist.gov/nist-header-footer/css/nist-combined.css\">\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/jquery-1.9.0.min.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/nist-header-footer.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            html[i.start:end]
        write(htmlfile, res)
        println("Inserting NIST header/footer into $htmlfile")
    end
    return htmlfile
end

addNISTHeaders(joinpath(@__DIR__, "build","index.html"))
addNISTHeaders.(map(name->joinpath(@__DIR__, "build", splitext(name)[1], "index.html"), map(p->p.second, pages[2:end])))

#map(name->rm(joinpath("src","$(splitext(name)[1]).md")), weavedocs)

# deploydocs(repo = "github.com/NeXLSpectrum.jl.git")


text:
# Provides extra verbose output for testing
function ev_error(client::Client, args...)
    id = "events.ev_error\t"
    clog(id, :yellow, "client $client ", :normal, " ", args..., "\n")
end
function ev_listen(port, args...)
    id = "events.ev_listen\t"
    nothing
end
function ev_connect(client::Client, args...)
    id = "events.ev_connect\t"
    clog(id, :yellow, " from client ", :bold, client , :normal, " ", args..., "\n")
end
function ev_close(client::Client, args...)
    id = "events.ev_close\t"
    clog(id, :yellow, " from client ", :bold, client , :normal, " ", args..., "\n")
end
function ev_write(client::Client, response::Response)
    id = "events.ev_write\t"
    clog(id, "to ", :yellow, :bold, client , :normal, "\n", response, "\n")
end
function ev_reset(client::Client, response::Response)
    id = "events.ev_reset\t"
    clog(id, :yellow, " client ", :bold, client , :normal, "\n", response, "\n")
end


text:
"""
    SeisMP(in,operator,param,iter,lp,imax)

Robust Matching Pursuit algorithm via lp-norm inner product

# Arguments
- `in`: Noisy data
- `operator`: Any operator passed the dot test
- `param`: parameters for the operator
- `iter`: Max iterations for the matching pursuit
- `lp`: lp-norm space inner product
- `imax`: Max interation for the lp-norm inner product

"""
function SeisMP(in,operator,param,iter,lp,imax)
    Madj=operator(in,param,-1,lp,imax);
    x=zeros(size(Madj));
    dr=zeros(size(in));
    r=in
    for i=1:iter
        Madj=operator(r,param,-1,lp,imax);
        ind=findmax(abs.(Madj))[2];
        M=zeros(size(Madj));
        M[ind[1],ind[2]]=Madj[ind[1],ind[2]];
        g=operator(M,param,1,lp,imax);
        alpha=LP_norm(vec(g),vec(r),lp,imax,1);
        dr=dr+alpha*g
        x=x+alpha*M;
        r=r-alpha*g;
    end
    return dr
end

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
export get_fields, get_records, build_result

function get_fields
end

function get_records
end

function build_result
end

try
    # interface for `DataFrame`s
    get_fields(df :: DataFrame) = names(df)
    get_records(df :: DataFrame) = zip(DataFrames.columns(df)...)
    function build_result(::Type{DataFrame}, fields, typs, source :: Union{Base.Generator, Vector{Tp}}) where Tp <: Tuple
        res = Tuple(typ[] for typ in typs)
        for each in source
            push!.(res, each)
        end
        DataFrame(collect(res), fields)
    end
catch
    @warn """
    DataFrames.jl not installed. If you want to peform queries on DataFrames with Wanderer.jl,
    build Wanderer.jl once again after installed DataFrames.jl.
    """
end

# interface for `Dict`s
get_fields(df :: D) where D <: Dict{Symbol} = collect(keys(df))
get_records(df :: D) where D <: Dict{Symbol} = zip((df[key] for key in keys(df))...)
function build_result(::Type{D}, fields, typs, source :: Union{Base.Generator, Vector{Tp}}) where {D <: Dict{Symbol}, Tp <: Tuple}
    res = Tuple(typ[] for typ in typs)
    for each in source
        push!.(res, each)
    end
    Dict(zip(fields, res))
end


text:
struct DataMH{vType<:Real}
	h::Vector{vType}
	m::Vector{vType}
end

struct SwParticle{psiType<:Real}
	ψ::psiType
end

function gr!(storage, x, h::Float64, p::MVector)
    storage[1] = -0.5*sin(psi-x[1])+h*sin(x[1])
end

function apply_field(ψ, ϕ, h)
    energy(x) = 0.5*sin(ψ-x[1])^2- h*cos.(x[1])
    gradient(storage, x) = gr!(storage, x, h, p)
    #optimization_result = optimize(energy, [ϕ], Newton())
    optimization_result = optimize(energy,[ϕ], BFGS())
    return Optim.minimizer(optimization_result)[1]
end

function obtain_hysteresis_loop(p::SwParticle)
	hstep = 0.01
    hmax = 2.0
    h = [collect(hmax:-hstep:-hmax); collect(-hmax:hstep:hmax)]
    len = length(h)
    m = zeros(len)

    ϕ = p.ψ
    for i=1:len
        ϕ = apply_field(p.ψ, ϕ, h[i])
        m[i] = cos(ϕ)
    end
	return DataMH(h, m)
end

function calculate_particles(psis, h)
    n = length(psis)
    len = length(h)
    phis = copy(psis)

    m = zeros(len, n)
    for p = 1:n
        for i = 1:len
            phis[p] = apply_field(psis[p], phis[p], h[i])
            m[i,p] = cos(phis[p])
        end
    end
    return m
end

@recipe function plot_one_sw_particle(data::DataMH)
	legend --> :none
	xlabel --> "Field"
	ylabel --> "Magnetization"
	xlim --> (-2.2, 2.2)
	ylim --> (-1.1,1.1)
    aspect_ratio --> [1 2]

    @series begin
        x = vcat((-2.2:0.01:2.2))
        y = vcat(zeros(441))
        linecolor --> :black
        linestyle --> :dash
        (x, y)
    end

    @series begin
        x = vcat(zeros(221))
        y = vcat((-1.1:0.01:1.1))
        linecolor --> :black
        linestyle --> :dash
        (x, y)
	end

	@series begin
        x = vcat((-1:0.01:1), ones(201), (1:-0.01:-1), -1*ones(201))
        y = vcat(ones(201), (1:-0.01:-1), -1*ones(201), (-1:0.01:1))
        linecolor --> :red
        linestyle --> :dash
        (x, y)
	end

    @series begin
        linecolor --> :blue
        linewidth --> 2
        (data.h, data.m)
	end
end


text:
module MultiQuad
using Hecke

function number_of_subgroups(p::Int, n::Int)
  @assert isprime(p)
  G = fmpz[1,2]
  pn = fmpz(p)
  for i=2:n
    push!(G, 2*G[i] + (pn -1)*G[i-1])
    pn *= p
  end
  return G[end]
end

function _combine(f::fmpq_poly, g::fmpq_poly, Qxy)
  Qx = parent(f)
  x = gen(Qx)
  y = gen(Qxy)
  f1 = f(x+y)
  f2 = g(y)
  return resultant(f1, f2)
end

function multi_quad_with_aut(d::Vector{fmpz})
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  Qxy, y = PolynomialRing(Qx, "y", cached = false)
  lp = [ NumberField(x^2-a)[1] for a = d]
  aut = [ [gen(x), -gen(x)] for x = lp]
  while length(lp) > 1
    ld = []
    lau = []
    for i=1:div(length(lp), 2)
      K, m1, m2 = compositum(lp[2*i-1], lp[2*i])
      @assert m1(gen(lp[2*i-1])) + m2(gen(lp[2*i])) == gen(K)
      au = [m1(x) + m2(y) for x = aut[2*i-1] for y = aut[2*i]]
      push!(ld, K)
      push!(lau, au)
    end
    if isodd(length(lp))
      push!(ld, lp[end])
      push!(lau, aut[end])
    end
    lp = ld
    aut = lau
  end
  return lp[1], aut[1]
end

function multi_quad_with_emb(d::Vector{fmpz})
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  Qxy, y = PolynomialRing(Qx, "y", cached = false)
  lp = [ NumberField(x^2-a)[1] for a = d]
  aut = [ [gen(x)] for x = lp]
  while length(lp) > 1
    ld = []
    lau = []
    for i=1:div(length(lp), 2)
      K, m1, m2 = compositum(lp[2*i-1], lp[2*i])
      push!(ld, K)
      push!(lau, vcat([ m1(x) for x = aut[2*i-1]], [ m2(x) for x = aut[2*i]]))
    end
    if isodd(length(lp))
      push!(ld, lp[end])
      push!(lau, aut[end])
    end
    lp = ld
    aut = lau
  end
  return lp[1], aut[1]
end

function multi_quad(d::Vector{fmpz}, B::Int)
  K, rt = multi_quad_with_emb(d)

  b = [K(1)]
  bb = [K(1)]
  for i = 1:length(d)
    if (d[i] < 0 && d[i] % 4 == -3) || (d[i] > 0 && d[i] % 4 == 1)
      o = (1+rt[i])//2
    else
      o = rt[i]
    end
    append!(b, Ref(o) .* b)
    append!(bb, Ref(rt[i]) .* bb)
  end

  all_d = fmpz[1]
  for i = d
    append!(all_d, Ref(i) .* all_d)
  end

  # @show all_d

  ZK = Order(K, b)
  ZK = pmaximal_overorder(ZK, fmpz(2))
  ZK.ismaximal = 1
  Hecke._set_maximal_order_of_nf(K, ZK)

  c = Hecke.class_group_init(ZK, B, complete = false, add_rels = false, min_size = 0)
  cp = Set(minimum(x) for x = c.FB.ideals)
  t_ord = 0
  local t_u

  Zx, x = FlintZZ["x"]

  for i = 2:length(all_d)
    k, a = number_field(x^2-all_d[i], cached = false)
    zk = maximal_order(k)
    class_group(zk)
    lp = prime_ideals_up_to(zk, Int(B), complete = false, degree_limit = 1)
    #only need split primes ...
    lp = [ x for x = lp if minimum(x) in cp]
    @assert all(x->minimum(x) == norm(x), lp)
    if length(lp) > 0
      S, mS = Hecke.sunit_group_fac_elem(lp)
    else
      S, mS = Hecke.unit_group_fac_elem(zk)
    end
    h = Hecke.NfToNfMor(k, K, bb[i])
    @assert bb[i]^2 == all_d[i]

    for i=2:ngens(S) # don't need torsion here - it's the "same" everywhere
      u = mS(S[i])  #do compact rep here???
      u = Hecke.compact_presentation(u, 2, decom = Dict((P, valuation(u, P)) for P = lp))
      Hecke.class_group_add_relation(c, FacElem(Dict((h(x), v) for (x,v) = u.fac)))
    end
    if t_ord < order(S[1])
      t_ord = order(S[1])
      t_u = FacElem(Dict((h(x), v) for (x,v) = mS(S[1]).fac))
    end
  end
  zeta = evaluate(t_u)
  z_all = [K(1)]
  for i=1:t_ord-1
    push!(z_all, z_all[end]*zeta)
  end
  Hecke._set_nf_torsion_units(K, (z_all, zeta))

  return c
end

function dlog(dl::Dict, x, p::Int)
  if iszero(x)
    throw(Hecke.BadPrime(1))
  end
  if haskey(dl, x)
    return dl[x]
  end
#  println("difficult for ", parent(x))
  i = 2
  y = x*x
  while !haskey(dl, y)
    y *= x
    i += 1
    @assert i <= p
  end
  #OK: we know x^i = g^dl[y] (we don't know g)
  v = dl[y]
  g = gcd(p, i)
  r = div(p, g)
  @assert v % g == 0
  e = invmod(div(i, g), r)*div(v, g) % r
  if e == 0
    e = r
  end
  dl[x] = e
  y = x*x
  f = (e*2) % p
  while !isone(y)
    if haskey(dl, y)
      @assert dl[y] == f
    end
    dl[y] = f
    y *= x
    f = (f+e) % p
  end
  g = [ a for (a,b) = dl if b == 1]
  @assert length(g) == 1
  @assert g[1]^dl[x] == x
  return dl[x]
end

function Hecke.matrix(R::Hecke.Ring, M::MatElem)
  return matrix(R, nrows(M), ncols(M), elem_type(R)[R(M[i,j]) for i=1:nrows(M) for j=1:ncols(M)])
end

function _nullspace(A::nmod_mat)
  A_orig = A
  p = fmpz(A.n)
  if isprime(p)
    return nullspace(A)
  end
  A = A'
  R = base_ring(A)
  r = nrows(A)
  c = ncols(A)
  A = hcat(A, identity_matrix(R, nrows(A)))
  A = vcat(A, zero_matrix(R, ncols(A) - nrows(A), ncols(A)))

  howell_form!(A)
  i = nrows(A)
  while iszero(A[i, :])
    i -= 1
  end
  r = i
  while i>0 && iszero(A[i:i, 1:c])
    i-= 1
  end
  if i < nrows(A)
    if i<r
      A = sub(A, i+1:r, c+1:ncols(A))
    else
      A = zero_matrix(base_ring(A), 0, ncols(A)-c)
    end
  else
    A = sub(A, i:r, c+1:ncols(A))
  end
  A = A'
  @assert iszero(A_orig * A)
  for i = keys(factor(p).fac)
    if valuation(p, i) > 1
      continue
    end
    b = matrix(ResidueRing(FlintZZ, Int(i)), A_orig)
    b = nullspace(b)
    b = rref(b[1]')
    c = matrix(base_ring(b[2]), A)'
    c = rref(c)
    if c[1] != b[1]
      global bla
      bla = A_orig, A, c, b
    end
    @assert c[1] == b[1]
  end
  return ncols(A), A
end

function mod_p(R, Q::NfOrdIdl, p::Int)
  F, mF = Hecke.ResidueFieldSmall(order(Q), Q)
  mF = Hecke.extend_easy(mF, nf(order(Q)))
  @assert size(F) % p == 1
  pp,e = Hecke.ppio(Int(size(F)-1), p)
#  @show pp, e, p
  dl = Dict{elem_type(F), Int}()
  dl[F(1)] = 0
#  #=
  lp = factor(p)
  while true
    x = rand(F)
    if iszero(x)
      continue
    end
    x = x^e
    if any(i-> x^div(pp, Int(i)) == 1, keys(lp.fac))
      continue
    else
      dlog(dl, x, pp)
      @assert length(dl) == pp
      break
    end
  end
#  =#
  return matrix(ResidueRing(FlintZZ, p), 1, length(R), [dlog(dl, mF(x)^e, pp) % p for x = R])
end

Hecke.lift(A::fmpz_mat) = A
#Lorenz: does not work for 8|n in general...
function saturate_exp(c::Hecke.ClassGrpCtx, p::Int, stable = 1.5)
  # p does NOT have to be a prime!!!

  ZK = order(c.FB.ideals[1])
  T = torsion_unit_group(ZK)[1]
  sT = Int(length(T))

  R = vcat(c.R_gen, c.R_rel)
  K = nf(ZK)
  _, zeta = Hecke._get_nf_torsion_units(K)
  if !(hash(zeta) in c.RS)
    #println("adding zeta = ", zeta)
    push!(R, zeta)
  else
    #println("NOT doint zeta")
  end
  T = ResidueRing(FlintZZ, p)
  A = identity_matrix(T, length(R))
  i = 1
  for (up, k) = factor(p).fac
    if sT % Int(up) == 0
      all_p = [up^i for i=1:k]
    else
      all_p = [up^k]
    end
    #@show all_p
    AA = identity_matrix(FlintZZ, ncols(A))
    for pp = all_p
      #println("doin' $pp")
      AA = matrix(ResidueRing(FlintZZ, Int(pp)), lift(AA))
      Ap = matrix(base_ring(AA), A)
      i = 1
      S = Hecke.PrimesSet(Hecke.p_start, -1, Int(pp), 1)
      cAA = ncols(AA)
      for q in S
        if isindex_divisor(ZK, q)
          continue
        end
        if discriminant(ZK) % q == 0
          continue
        end
        #if gcd(div(q-1, Int(pp)), pp) > 1 # not possible if cond(k) is involved
        #  continue
        #end
        lq = prime_decomposition(ZK, q, 1)
        if length(lq) == 0
          continue
        end
        for Q in lq
          try
            z = mod_p(R, Q[1], Int(pp))
            z = z*Ap
            z = _nullspace(z)
            B = hcat(AA, sub(z[2], 1:nrows(z[2]), 1:z[1]))
            B = _nullspace(B)
            AA = AA*sub(B[2], 1:ncols(AA), 1:B[1])
            if !isprime(p)
              AA = AA'
              if nrows(AA) < ncols(AA)
                AA = vcat(AA, zero_matrix(base_ring(AA), ncols(AA) - nrows(AA), ncols(AA)))
              end
              howell_form!(AA)
              local i = nrows(AA)
              while i>0 && iszero(AA[i, :])
                i -= 1
              end
              AA = sub(AA, 1:i, 1:ncols(AA))'
            else
              @assert rank(AA') == ncols(AA)
            end
#            @show cAA, pp, q, size(AA)
            if cAA == ncols(AA)
              break #the other ideals are going to give the same info
                    #for multi-quad as the field is normal
            end
          catch e
            @show "BadPrime"
            if !isa(e, Hecke.BadPrime)
              rethrow(e)
            end
          end
        end
        if cAA == ncols(AA)
          #println("good $i")
          i += 1
        else
          #println("bad")
          i = 0
        end
        cAA = ncols(AA)
        if i> stable*ncols(AA)
          break
        end
      end
    end
    pp = Int(modulus(base_ring(AA)))
    #@show "done $pp"
    # A is given mod p, AA mod pp
    #we need AA mod p where the lift is any lift, modulo powers of pp
    #                                   identity modulo coprime (CRT)
    q, w = Hecke.ppio(p, pp) # q is a "power" of pp and w is coprime
    g, e, f = gcdx(q, w)
    AA = AA'
    AA = vcat(AA, zero_matrix(base_ring(AA), ncols(AA) - nrows(AA), ncols(AA)))
    strong_echelon_form!(AA)

    X  = similar(AA)
    for j=1:min(nrows(X), ncols(X))
      X[j,j] = 1
    end
    _A = matrix(base_ring(A), e*q*lift(X) + f*w*lift(AA))
    A = _A*A'
    howell_form!(A)
    i = nrows(A)
    while iszero(A[i, :])
      i -= 1
    end
    A = sub(A, 1:i, 1:ncols(A))'
    #@show size(A)
  end
  return A
end

fe(a::FacElem) = a
fe(a::nf_elem) = FacElem(a)

function elems_from_sat(c::Hecke.ClassGrpCtx, z)
  res = []
  fac = []
  for i=1:ncols(z)
    a = fe(c.R_gen[1])^FlintZZ(z[1, i])
    b = FlintZZ(z[1, i]) * c.M.bas_gens[1]
    for j=2:length(c.R_gen)
      a *= fe(c.R_gen[j])^FlintZZ(z[j, i])
      b += FlintZZ(z[j, i]) * c.M.bas_gens[j]
    end
    for j=1:length(c.R_rel)
      a *= fe(c.R_rel[j])^FlintZZ(z[j + length(c.R_gen), i])
      b += FlintZZ(z[j + length(c.R_gen), i]) * c.M.rel_gens[j]
    end

    push!(res, (a, b))
  end
  return res
end

function saturate(c::Hecke.ClassGrpCtx, n::Int, stable = 3.5)
  e = matrix(FlintZZ, saturate_exp(c, n%8 == 0 ? 2*n : n, stable))
  se = sparse_matrix(e)'

  A = sparse_matrix(FlintZZ)
  K = nf(c)
  _, zeta = Hecke._get_nf_torsion_units(K)

  #println("Enlarging by $(ncols(e)) elements")
  n_gen = []
  for i=1:ncols(e)
    r = e[:, i]
    g = content(r)
    g = gcd(g, n)
    divexact!(r, r, g)
#    g == 1 || println("non triv. content $g in ", e[:, i])
    a = fe(c.R_gen[1])^r[1, 1]
    fac_a = r[1, 1] * c.M.bas_gens[1]
    for j = 2:length(c.R_gen)
      a *= fe(c.R_gen[j])^r[j, 1]
      fac_a += r[j, 1] * c.M.bas_gens[j]
    end
    for j=1:length(c.R_rel)
      a *= fe(c.R_rel[j])^r[j + length(c.R_gen), 1]
      fac_a += r[j + length(c.R_gen), 1] * c.M.rel_gens[j]
    end
    if nrows(e) > length(c.R_gen) + length(c.R_rel)
      @assert length(c.R_gen) + length(c.R_rel) + 1 == nrows(e)
      a *= fe(zeta)^r[nrows(e), 1]
    end

    decom = Dict((c.FB.ideals[k], v) for (k,v) = fac_a)
    if n == g
      fl = true
      x = a
    else
      fl, x = ispower(a, div(n, Int(g)), decom = decom)
      if !fl
        fl, x = ispower(nf(c)(-1)*a, div(n, Int(g)), decom = decom)
      end
    end
    if fl
      push!(n_gen, FacElem(x))
      r = divexact(se.rows[i], g)
      push!(r.pos, nrows(e) + length(n_gen))
      push!(r.values, -div(n, Int(g)))
      push!(A, r)
    else
      global bad = (a, div(n, Int(g)))
      error("not a power")
    end
  end

  #= Idea - before I forget:
  we have generators g_1, ..., g_n on input and enlarge by
                     h_1, ..., h_r
  And we have relations: n*h_i = some word in g
  in matrices:
  A = (words in g | n*I)
  now, using the column HNF
  A T = H = (I|0) - if the relations were primitive
  Thus
  A * (g | h)^t = 0 (using the relations) (possibly missing a sign)
  T^-1 = (R|S)^t
  then
  A T T^-1 (g|h)^t = (I|0) T^-1 (g|h)^t
  => R^t (g|h)^t = 0
  => S^t (g|h) is the new basis (by dimensions)

  now: since Hecke is row based, we have to transpose..
  =#
  A = A'
#    @show fmpz_mat(A)
  H, T = hnf_with_trafo(fmpz_mat(A))
  @assert isone(sub(H, 1:ncols(A), 1:ncols(A))) #otherwise, relations sucked.
  Ti = inv(T')
  Ti = sub(Ti, length(n_gen)+1:nrows(Ti), 1:ncols(Ti))

  R = vcat(c.R_gen, c.R_rel)
  if !(hash(zeta) in c.RS)
    push!(R, zeta)
  end
  R = vcat(R, n_gen)
  @assert ncols(Ti) == length(R)

  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)

  for i=1:nrows(Ti)
    a = FacElem(K(1))
    for j=1:ncols(Ti)
      a *= R[j]^Ti[i, j]
    end
      #TODO remove zeta from relations!!
    Hecke.class_group_add_relation(d, a)
  end

  return d
end

function sunits_mod_units(c::Hecke.ClassGrpCtx)
  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo
  su = Vector{FacElem{nf_elem, AnticNumberField}}()
  for i=1:length(c.FB.ideals)
    x = zeros(fmpz, length(c.R_gen) + length(c.R_rel))
    x[i] = 1
    for j in length(trafos):-1:1
      Hecke.apply_right!(x, trafos[j])
    end
    y = FacElem(vcat(c.R_gen, c.R_rel), x)
    push!(su, y)
  end
  return su
end

function simplify(c::Hecke.ClassGrpCtx)
  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)
  U = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(order(d))

  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo

  for i=1:length(c.FB.ideals)
    x = zeros(fmpz, length(c.R_gen) + length(c.R_rel))
    x[i] = 1
    for j in length(trafos):-1:1
      Hecke.apply_right!(x, trafos[j])
    end
    y = FacElem(vcat(c.R_gen, c.R_rel), x)
    fl = Hecke.class_group_add_relation(d, y, deepcopy(c.M.basis.rows[i]))
    @assert fl
  end
  for i=1:nrows(c.M.rel_gens)
    if iszero(c.M.rel_gens.rows[i])
      Hecke._add_unit(U, c.R_rel[i])
    end
  end
  for i=1:length(U.units)
    Hecke.class_group_add_relation(d, U.units[i], SRow(FlintZZ))
  end
  return d
end

function units(c::Hecke.ClassGrpCtx)
  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)
  U = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(order(d))

  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo

  for i=1:nrows(c.M.rel_gens)
    if iszero(c.M.rel_gens.rows[i])
      Hecke._add_unit(U, c.R_rel[i])
    end
  end

  U.units = Hecke.reduce(U.units, U.tors_prec)
  U.tentative_regulator = Hecke.regulator(U.units, 64)

  return U
end

#TODO:
#  use the essential part only for the saturation (pointless for MultiQuad:
#    the Brauer relations force every prime block to have 2 on the
#    diagonal)
#  in MultiQuad, we get the "unit-part" for free without the expensive
#    real-part, so the S-units are the image of the rel mat, and
#    no need for the kernel.
#  keep track if the relations and update the unit group as well
#  S-units: easy: add the relation, here if n is prime, either
#   the S-class number or the regulator changes, never both
#  units: we have new^n = prod old. use this to obtain new basis
#
#  track the torsion as well
#  if n is divisible by 8, then, generically, the saturation needs to
#  be followed by a second saturation at 2:
#    Elements look like (locally) an 8-th power but are only a 4-th
#    so I can only extract a 4-th.
#    However, it might be an 8-th (or the product of 2 might be an 8-th)
#  Darn. Math is unfair.
#
#  extend to gen. mult group...
end


text:
using OrdinaryDiffEq, Test
f(u,p,t) = 1.01*u
u0=1/2
tspan = (0.0,Inf)
prob = ODEProblem(f,u0,tspan)

# Shouldn't error, but should go unstable and abort
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,adaptive=false,dt=0.1)
@test sol.retcode == :Unstable


text:
"""
    intersect(a::Line, b::Line) -> Point

Intersection of 2 line segmens `a` and `b`.
Returns intersection_found::Bool, intersection_point
"""
function intersects(a::Line{2, T1}, b::Line{2, T2}) where {T1, T2}
    T = promote_type(T1, T2)
    v1, v2 = a; v3, v4 = b;
    MT = Mat{2,2,T,4}
    p0 = zero(Point2{T})

    verticalA = v1[1] == v2[1]
    verticalB = v3[1] == v4[1]

    # if a segment is vertical the linear algebra might have trouble
    # so we will rotate the segments such that neither is vertical
    dorotation = verticalA || verticalB

    if dorotation
        θ = T(0.0)
        if verticalA && verticalB
            θ = T(π / 4)
        elseif verticalA || verticalB # obviously true, but make it clear
            θ34 = -atan(v4[2] - v3[2], v4[1] - v3[1])
            θ12 = -atan(v2[2] - v1[2], v2[1] - v1[1])
            θ = verticalA ? θ34 : θ12
            θ = abs(θ) == T(0) ? (θ12 + θ34) / 2 : θ
            θ = abs(θ) == T(pi) ? (θ12 + θ34) / 2 : θ
        end
        rotation = MT(cos(θ), sin(θ), -sin(θ), cos(θ))
        v1 = rotation * v1
        v2 = rotation * v2
        v3 = rotation * v3
        v4 = rotation * v4
    end

    a = det(MT(v1[1] - v2[1], v1[2] - v2[2],
               v3[1] - v4[1], v3[2] - v4[2]))

    (abs(a) < eps(T)) && return false, p0 # Lines are parallel

    d1 = det(MT(v1[1], v1[2], v2[1], v2[2]))
    d2 = det(MT(v3[1], v3[2], v4[1], v4[2]))
    x = det(MT(d1, v1[1] - v2[1], d2, v3[1] - v4[1])) / a;
    y = det(MT(d1, v1[2] - v2[2], d2, v3[2] - v4[2])) / a;

    (x < prevfloat(min(v1[1], v2[1])) || x > nextfloat(max(v1[1], v2[1]))) && return false, p0
    (y < prevfloat(min(v1[2], v2[2])) || y > nextfloat(max(v1[2], v2[2]))) && return false, p0
    (x < prevfloat(min(v3[1], v4[1])) || x > nextfloat(max(v3[1], v4[1]))) && return false, p0
    (y < prevfloat(min(v3[2], v4[2])) || y > nextfloat(max(v3[2], v4[2]))) && return false, p0

    point = Point2{T}(x, y)
    # don't forget to rotate the answer back
    if dorotation
        point = transpose(rotation) * point
    end

    return true, point
end

function simple_concat(vec::AbstractVector, range, endpoint::P) where P
    result = Vector{P}(undef, length(range) + 1)
    for (i, j) in enumerate(range)
        result[i] = vec[mod1(j, length(vec))]
    end
    result[end] = endpoint
    return result
end

function consecutive_pairs(arr)
    n = length(arr)
    zip(view(arr, 1:n-1), view(arr, 2:n))
end

"""
    self_intersections(points::AbstractVector{AbstractPoint})

Finds all self intersections of polygon `points`
"""
function self_intersections(points::AbstractVector{<:AbstractPoint})
    sections = similar(points, 0)
    intersections = Int[]

    wraparound(i) = mod1(i, length(points) - 1)

    for (i, (a,b)) in enumerate(consecutive_pairs(points))
        for (j, (a2, b2)) in enumerate(consecutive_pairs(points))
            is1, is2 = wraparound(i+1), wraparound(i-1)
            if i != j && is1 != j && is2 != j && !(i in intersections) && !(j in intersections)
                intersected, p = intersects(Line(a,b), Line(a2, b2))
                if intersected
                    push!(intersections, i, j)
                    push!(sections, p)
                end
            end
        end
    end
    return intersections, sections
end

"""
    split_intersections(points::AbstractVector{AbstractPoint})

Splits polygon `points` into it's self intersecting parts. Only 1 intersection
is handled right now.
"""
function split_intersections(points::AbstractVector{<:AbstractPoint})
    intersections, sections = self_intersections(points)
    if isempty(intersections)
        return [points]
    elseif length(intersections) == 2 && length(sections) == 1
        a, b = intersections
        p = sections[1]
        a,b = min(a,b), max(a,b)
        poly1 = simple_concat(points, (a+1):(b-1), p)
        poly2 = simple_concat(points, (b+1):(length(points)+a), p)
        return [poly1, poly2]
    else
        error("More than 1 intersections can't be handled currently. Found: $intersections, $sections")
    end
end


text:
using BinaryBuilder
import SHA

const version = v"1.12.0"
const URL_PREFIX = "https://github.com/felipenoris/mongo-c-driver-builder/releases/download/v$version"

download_info = Dict()

const platforms_whitelist = [ Windows(:i686), Windows(:x86_64) ]
const platforms = filter( x -> x ∉ platforms_whitelist, supported_platforms())

const artefact = "libmongoc"

@sync for platform in platforms
        filename = "$artefact.v$version.$(triplet(platform)).tar.gz"
        url = "$URL_PREFIX/$filename"
        @async begin
            run(`wget $url`)

            open(filename, "r") do f
                sha = bytes2hex(SHA.sha2_256(f))
                download_info[platform] = (url, sha)
            end
            rm(filename)
        end
end

println("==================")
println("download_info for $artefact")
println("")
println("download_info = Dict(")

for v in download_info
    println("    " * string(v[1]) * " => " * string(v[2]) * ",")
end

println(")")
println("")

max_stars_repo_name: felipenoris/mongo-c-driver-builder

max_issues_repo_name: felipenoris/mongo-c-driver-builder

max_forks_repo_name: felipenoris/mongo-c-driver-builder


text:
using StructuresKit

prop = [100 29500.00 29500.00 0.30 0.30 11346.15]

node = [1 5.0000 1.0000 1 1 1 1 -38.889
2 5.0000 0.0000 1 1 1 1 -50.000
3 2.5000 0.0000 1 1 1 1 -50.000
4 0.0000 0.0000 1 1 1 1 -50.000
5 0.0000 3.0000 1 1 1 1 -16.667
6 0.0000 6.0000 1 1 1 1 16.667
7 0.0000 9.0000 1 1 1 1 50.000
8 2.5000 9.0000 1 1 1 1 50.000
9 5.0000 9.0000 1 1 1 1 50.000
10 5.0000 8.0000 1 1 1 1 38.889]

elem = [1 1 2 0.100000 100
2 2 3 0.100000 100
3 3 4 0.100000 100
4 4 5 0.100000 100
5 5 6 0.100000 100
6 6 7 0.100000 100
7 7 8 0.100000 100
8 8 9 0.100000 100
9 9 10 0.100000 100]

springs = []
constraints = 0

BC = "S-S"

lengths = [1.00
2.00
3.00
4.00
5.00
6.00
7.00
8.00
9.00
10.00
11.00
12.00
13.00
14.00
15.00
20.00
30.00
40.00
50.00
60.00
70.00
80.00
90.00
100.00
200.00
300.00
400.00
500.00
600.00
700.00
800.00
900.00
1000.00]

m_all = ones(length(lengths))

neigs = 20

curve, shapes = CUFSM.strip(prop,node,elem,lengths,springs,constraints,neigs)


text:
using PPInterpolation
import AQFED.TermStructure:
    TSBlackModel, varianceByLogmoneyness, discountFactor, logForward, forward

#Runge-Kutta-Legendre FDM for the American option.

function makeFDMPriceInterpolation(isCall, isEuropean, model, T, strike, N, M; method = "RKL2", sDisc = "Sinh", useExponentialFitting = false, smoothing = "Kreiss", lambdaS = 0.25, Xdev = 4, Smax = 0, rklStages = 0, epsilonRKL = 0, useSqrt = true)
    upwindingThreshold = 1.0
    v0 = varianceByLogmoneyness(model, 0.0, T)
    Xspan = Xdev * sqrt(v0 * T)
    lnfK = logForward(model, log(strike), T)
    Xmin = lnfK - Xspan - 0.5 * v0 * T
    Xmax = lnfK + Xspan - 0.5 * v0 * T
    Smin = exp(Xmin)
    if Smax == 0
        Smax = exp(Xmax)
    end
    X = collect(range(Xmin, stop = Xmax, length = M))
    hm = X[2] - X[1]
    Sscale = strike * lambdaS
    if sDisc == "Exp"
        S = exp.(X)
        J = exp.(X)
        Jm = @. exp(X - hm / 2)
    elseif sDisc == "Sinh"
        u = collect(range(0.0, stop = 1.0, length = M))
        Smin = 0.0
        c1 = asinh((Smin - strike) / Sscale)
        c2 = asinh((Smax - strike) / Sscale)
        S = @. strike + Sscale * sinh((c2 - c1) * u + c1)
        hm = u[2] - u[1]
        J = @. Sscale * (c2 - c1) * cosh((c2 - c1) * u + c1)
        Jm = @. Sscale * (c2 - c1) * cosh((c2 - c1) * (u - hm / 2) + c1)
    else #if sDisc == "Linear"
        S = collect(range(0.0, stop = Smax, length = M))
        X = S
        hm = X[2] - X[1]
        J = ones(M)
        Jm = ones(M)
    end
    Smin = S[1]
    sign = 1
    if !isCall
        sign = -1
    end

    F0 = zeros(M)
    if isCall
        F0 = @. max(S - strike, 0.0)
    else
        F0 = @. max(strike - S, 0.0)
    end
    iStrike = searchsortedfirst(S, strike)

    if smoothing == "Averaging"
        if strike < (S[iStrike] + S[iStrike-1]) / 2
            iStrike -= 1
        end
        a = (S[iStrike] + S[iStrike+1]) / 2
        if !isCall
            a = (S[iStrike] + S[iStrike-1]) / 2   # int a,lnK K-eX dX = K(a-lnK)+ea-K
        end
        value = (a - strike) * (a - strike) * 0.5
        F0[iStrike] = value / (S[iStrike+1] - S[iStrike-1]) * 2
    elseif smoothing == "Kreiss"
        xmk = S[iStrike] - strike
        h = (S[iStrike+1] - S[iStrike-1]) / 2
        # F0smooth[iStrike] = 0.5*xmk + h/6 + 0.5*xmk*xmk/h*(1-xmk/(3*h))
        F0[iStrike] = 0.5 * xmk + h / 6 + 0.5 * xmk * xmk / h * (1 - xmk / (3 * h))
        if !isCall
            F0[iStrike] -= xmk # C-P = f-K
        end
        iStrike -= 1
        xmk = S[iStrike] - strike
        h = (S[iStrike+1] - S[iStrike-1]) / 2
        F0[iStrike] = 0.5 * xmk + h / 6 + 0.5 * xmk * xmk / h * (1 + xmk / (3 * h))
        if !isCall
            F0[iStrike] -= xmk # C-P = f-K
        end
    end
    F = zeros(M)
    F[1:M] = F0

    lowerBoundA = zeros(M)
    if !isEuropean
        lowerBoundA .= F
    end
    ti = T
    dt = ti / N
    if useSqrt
        ti = 0.0
        dt = sqrt(T) / N
    end
    useDirichlet = false
    lbValue = computeLowerBoundary(isCall, strike, useDirichlet, 0.0, Smin, useSqrt ? ti^2 : ti, model)
    updatePayoffExplicitTrans(F, useDirichlet, lbValue, M)
    A1ilj = zeros(M)
    A1ij = zeros(M)
    A1iuj = zeros(M)
    if method == "RKL2"
        makeSystem(model, A1ilj, A1ij, A1iuj, useExponentialFitting, upwindingThreshold, dt, S, J, Jm, hm, useDirichlet, M)
        if useSqrt
            tih = ti + dt / 2
            @. A1ij *= 2tih
            @. A1ilj *= 2tih
            @. A1iuj *= 2tih
        end
        s, a, b, w0, w1 = initRKLCoeffs(dt, A1ij, epsilonRKL = epsilonRKL, rklStages = rklStages)
        Y0 = zeros(M)
        Y1 = zeros(M)
        Y2 = zeros(M)
        for n = 1:N
            tih = useSqrt ? ti + dt / 2 : ti - dt / 2
            lbValue = computeLowerBoundary(isCall, strike, useDirichlet, 0.0, Smin, useSqrt ? T - tih^2 : tih, model)
            F .= RKLStep(s, a, b, w0, w1, A1ilj, A1ij, A1iuj, F, Y0, Y1, Y2, useDirichlet, lbValue, lowerBoundA, M)
            if useSqrt
                ti += dt
                #  println(ti^2, " ", dt)
            else
                ti -= dt
            end

            if n < N
                makeSystem(model, A1ilj, A1ij, A1iuj, useExponentialFitting, upwindingThreshold, dt, S, J, Jm, hm, useDirichlet, M)
                if useSqrt
                    tih = ti + dt / 2
                    @. A1ij *= 2tih
                    @. A1ilj *= 2tih
                    @. A1iuj *= 2tih
                    s, a, b, w0, w1 = initRKLCoeffs(dt, A1ij, epsilonRKL = epsilonRKL, rklStages = rklStages)
                end
            end
        end
    end
    spl = makeCubicPP(S, F, PPInterpolation.SECOND_DERIVATIVE, 0.0, PPInterpolation.SECOND_DERIVATIVE, 0.0, PPInterpolation.VanLeer())
    return spl
end

function initRKLCoeffs(dt, A1ij; epsilonRKL = 0.0, rklStages = 0)
    dtexplicit = dt / max(maximum(A1ij))
    dtexplicit /= 2 #lambdaS
    s = 0.0
    delta = 1 + 4 * (2 + 4 * dt / dtexplicit)
    s = ceil(Int, (-1 + sqrt(delta)) / 2)
    if s % 2 == 0
        s += 1
    end
    if epsilonRKL > 0
        s = computeRKLStages(dtexplicit, dt, epsilonRKL)
    end
    if rklStages > 0
        s = rklStages
    end
    # println("s=",s)
    a = zeros(s)
    b = zeros(s)
    w0 = 1.0
    w1 = 0.0
    if epsilonRKL == 0
        w1 = 4 / (s^2 + s - 2)
        b[1] = 1.0 / 3
        b[2] = 1.0 / 3
        a[1] = 1.0 - b[1]
        a[2] = 1.0 - b[2]
        for i = 3:s
            b[i] = (i^2 + i - 2.0) / (2 * i * (i + 1.0))
            a[i] = 1.0 - b[i]
        end
    else
        w0 = 1 + epsilonRKL / s^2
        _, tw0p, tw0p2 = legPoly(s, w0)
        w1 = tw0p / tw0p2
        b = zeros(s)
        for jj = 2:s
            _, tw0p, tw0p2 = legPoly(jj, w0)
            b[jj] = tw0p2 / tw0p^2
        end
        b[1] = b[2]
        a = zeros(s)
        for jj = 2:s
            tw0, _, _ = legPoly(jj - 1, w0)
            a[jj-1] = (1 - b[jj-1] * tw0)
        end
    end
    return s, a, b, w0, w1
end

function makeSystem(model, A1ilj::AbstractArray{T}, A1ij::AbstractArray{T}, A1iuj::AbstractArray{T}, useExponentialFitting::Bool, upwindingThreshold::Real, dt::Real, S::Vector{T}, J::Vector{T}, Jm::Vector{T}, hm::Real, useDirichlet::Bool, M::Int) where {T}
    if useDirichlet
        A1ij[1] = 0
    else
        drifti = (model.r - model.q) * S[1]
        A1iuj[1] = -dt * drifti / (Jm[2] * hm)
        A1ij[1] = -dt * (-model.r * 0.5) - A1iuj[1]
    end
    drifti = (model.r - model.q) * S[M]
    A1ilj[M] = dt * drifti / (Jm[M] * hm)
    A1ij[M] = -dt * (-model.r * 0.5) - A1ilj[M]

    @inbounds @simd for i = 2:M-1
        svi = S[i]^2 * model.vol^2 / J[i]
        drifti = (model.r - model.q) * S[i]
        if useExponentialFitting
            if abs(drifti * hm / svi) > upwindingThreshold
                svi = drifti * hm / tanh(drifti * hm / svi)
            end
        end
        svi /= hm^2
        drifti /= (2 * J[i] * hm)
        A1iuj[i] = -dt * (svi / (2Jm[i+1]) + drifti)
        A1ij[i] = -dt * (-svi / 2 * (1 / Jm[i+1] + 1 / Jm[i]) - model.r)
        A1ilj[i] = -dt * (svi / (2Jm[i]) - drifti)
    end
end

function explicitStep(a::Real, A1ilj::Vector{T}, A1ij::Vector{T}, A1iuj::Vector{T}, F::Vector{T}, Y0::Vector{T}, Y1::Vector{T}, M::Int) where {T}
    index = 1
    Y1[index] = Y0[index] - a * (A1ij[index] * F[index] + A1iuj[index] * F[index+1])
    index = M
    Y1[index] = Y0[index] - a * (A1ij[index] * F[index] + A1ilj[index] * F[index-1])
    @inbounds @simd for index = 2:M-1
        Y1[index] = Y0[index] - a * (A1ij[index] * F[index] + A1iuj[index] * F[index+1] + A1ilj[index] * F[index-1])
    end
end

function RKLStep(s::Int, a::Vector{T}, b::Vector{T}, w0::Real, w1::Real, A1ilj::Vector{T}, A1ij::Vector{T}, A1iuj::Vector{T}, F::Vector{T}, Yjm2::Vector{T}, Yjm::Vector{T}, Yj::Vector{T}, useDirichlet::Bool, lbValue::Real, lowerBound::Vector{T}, M::Int) where {T}
    mu1b = b[1] * w1
    explicitStep(mu1b, A1ilj, A1ij, A1iuj, F, F, Yjm, M)
    updatePayoffExplicitTrans(Yjm, useDirichlet, lbValue, M)
    MY0 = (Yjm - F) / mu1b
    enforceLowerBound(Yjm, lowerBound, M)
    Yjm2 .= F
    for j = 2:s
        muu = (2 * j - 1) * b[j] / (j * b[j-1])
        muj = muu * w0
        mujb = muu * w1
        gammajb = -a[j-1] * mujb
        nuj = -1.0 * b[2] / (2.0 * b[1]) #b0 = b[1]
        if j > 2
            nuj = -(j - 1) * b[j] / (j * b[j-2])
        end
        @. Yj = muj * Yjm + nuj * Yjm2 + (1 - nuj - muj) * F + gammajb * MY0 # + mujb*MYjm
        explicitStep(mujb, A1ilj, A1ij, A1iuj, Yjm, Yj, Yj, M)
        updatePayoffExplicitTrans(Yj, useDirichlet, lbValue, M)
        enforceLowerBound(Yj, lowerBound, M)
        Yjm2, Yjm = Yjm, Yjm2
        Yjm, Yj = Yj, Yjm
    end
    return Yjm
end

function enforceLowerBound(F::Vector{T}, lowerBound::Vector{T}, M::Int) where {T}
    if length(lowerBound) > 0
        @. F = max(F, lowerBound)
    end
end

function computeLowerBoundary(isCall::Bool, strike::Real, useDirichlet::Bool, B::Real, Smin::Real, ti::Real, model)
    lbValue = 0.0
    if useDirichlet && B == 0
        if !isCall
            lbValue = (strike - forward(model, Smin, ti) * discountFactor(model, ti))
        end
    end
    return lbValue
end

function updatePayoffExplicitTrans(F::Vector{T}, useDirichlet::Bool, lbValue::Real, M::Int) where {T}
    if useDirichlet
        F[1] = lbValue
    end
end

function computeRKLStages(dtexplicit, dt, ep)
    s = 1
    betaFunc = function (s::Int)
        w0 = 1 + ep / s^2
        _, tw0p, tw0p2 = legPoly(s, w0)
        beta = (w0 + 1) * tw0p2 / tw0p
        return beta - 2 * dt / (dtexplicit)
    end
    while s < 10000 && betaFunc(s) < 0
        s += 1
    end
    #s += Int(ceil(s/10))
    return s
end

function legPoly(s::Int, w0::Real)
    tjm = 1.0
    tj = w0
    if s == 1
        return tj, 1.0, 0.0
    end
    dtjm = 0.0
    dtj = 1.0
    d2tjm = 0.0
    d2tj = 0.0

    for j = 2:s
        onej = 1.0 / j
        tjp = (2 - onej) * w0 * tj - (1 - onej) * tjm
        dtjp = (2 - onej) * (tj + w0 * dtj) - (1 - onej) * dtjm
        d2tjp = (2 - onej) * (dtj * 2 + w0 * d2tj) - (1 - onej) * d2tjm
        tjm = tj
        dtjm = dtj
        d2tjm = d2tj
        tj = tjp
        dtj = dtjp
        d2tj = d2tjp
    end
    return tj, dtj, d2tj
end


text:
function Collocation_Matrix(N)

      #Radau
      if  N == 2

      t1 = 0.3333333
      t2 = 1.0

      M1 = [
            t1 1 / 2 * t1^2
            t2 1 / 2 * t2^2
            ]
      M2 = [
            1 t1
            1 t2
            ]

      M = M1 * inv(M2)

      elseif  N == 3

            t1 = 0.155051
            t2 = 0.644949
            t3 = 1.0

            M1 = [
                  t1 1 / 2 * t1^2 1 / 3 * t1^3
                  t2 1 / 2 * t2^2 1 / 3 * t2^3
                  t3 1 / 2 * t3^2 1 / 3 * t3^3
                  ]
            M2 = [
                  1 t1 t1^2
                  1 t2 t2^2
                  1 t3 t3^2
                  ]

            M = M1 * inv(M2)

      elseif N == 4

            t1 = 0.088588;
            t2 = 0.409467;
            t3 = 0.787659;
            t4 = 1;

            M1 = [
                  t1 1 / 2 * t1^2 1 / 3 * t1^3 1 / 4 * t1^4
                  t2 1 / 2 * t2^2 1 / 3 * t2^3 1 / 4 * t2^4
                  t3 1 / 2 * t3^2 1 / 3 * t3^3 1 / 4 * t3^4
                  t4 1 / 2 * t4^2 1 / 3 * t4^3 1 / 4 * t4^4
                  ]
            M2 = [
                  1 t1 t1^2 t1^3
                  1 t2 t2^2 t2^3
                  1 t3 t3^2 t3^3
                  1 t4 t4^2 t4^3
                  ]

            M = M1 * inv(M2)

      elseif N == 5

            t1 = 0.057104;
            t2 = 0.276843;
            t3 = 0.583590;
            t4 = 0.860240;
            t5 = 1;

            M1 = [
                  t1 1 / 2 * t1^2 1 / 3 * t1^3 1 / 4 * t1^4  1 / 5 * t1^5
                  t2 1 / 2 * t2^2 1 / 3 * t2^3 1 / 4 * t2^4  1 / 5 * t2^5
                  t3 1 / 2 * t3^2 1 / 3 * t3^3 1 / 4 * t3^4  1 / 5 * t3^5
                  t4 1 / 2 * t4^2 1 / 3 * t4^3 1 / 4 * t4^4  1 / 5 * t4^5
                  t4 1 / 2 * t4^2 1 / 3 * t4^3 1 / 4 * t4^4  1 / 5 * t5^5
                  ]
            M2 = [
                  1 t1 t1^2 t1^3 t1^4
                  1 t2 t2^2 t2^3 t2^4
                  1 t3 t3^2 t3^3 t3^4
                  1 t4 t4^2 t4^3 t4^4
                  1 t4 t4^2 t4^3 t5^4
                  ]

            M = M1 * inv(M2)

      end

      return M
end


text:
function ec_bar_init()
    opt = ec_init_cartesian()
    opt["tooltip"]["trigger"] = "axis"
    opt["toolbox"]["feature"]["magicType"] =
        Dict{ASCIIString,Any}(
            "show" => true,
            "type" => ["line", "bar"]
         )
    opt["calculable"] = true
    opt["yAxis"] = [
            Dict{ASCIIString,Any}(
                "splitArea" => Dict{ASCIIString,Bool}("show" => true),
                "type" => "value"
            )
        ]
    opt["xAxis"] = [
            Dict{ASCIIString,Any}("type" => "category")
        ]

    return opt
end

function ec_bar_v{S<:String, N<:Real}(legendlist::Vector{S}, data::Array{N}, labels=[1:size(data, 1)]; smooth::Bool = true)
    opt = ec_bar_init()
    opt["legend"]["data"] = legendlist

    n = size(data, 2)
    resize!(opt["series"], n)

    for i = 1:n
        opt["series"][i] = Dict(
            "name" => legendlist[i],
            "type" => "bar",
            "smooth" => smooth,
            "data" => data[:, i]
            )
    end
    opt["xAxis"][1]["data"] = labels

    return opt
end

function ec_bar_h{S<:String, N<:Real}(legendlist::Vector{S}, data::Array{N}, labels=[1:size(data, 1)]; smooth::Bool=true)
    opt = ec_bar_v(legendlist, data, labels; smooth=smooth)

    tmp = opt["xAxis"]
    opt["xAxis"] = opt["yAxis"]
    opt["yAxis"] = tmp

    return opt
end

function ec_bar_piled{S<:String, N<:Real}(legendlist::Vector{S}, data::Array{N}, labels=[1:size(data, 1)]; smooth::Bool=true)
    opt = ec_bar(legendlist, data, labels; smooth=smooth)

    for i = 1:length(opt["series"])
        opt["series"][i]["stack"] = "total"
    end

    return opt
end

function ec_bar_piled_h{S<:String, N<:Real}(legendlist::Vector{S}, data::Array{N}, labels=[1:size(data, 1)]; smooth::Bool=true)
    opt = ec_bar_h(legendlist, data, labels; smooth=smooth)

    for i = 1:length(opt["series"])
        opt["series"][i]["stack"] = "total"
    end

    return opt
end

function ec_waterfall{S<:String, N<:Real}(legendlist::Vector{S}, data::Vector{N}, labels=[1:size(data, 1)]; smooth::Bool=false)
    opt = ec_init_cartesian()
    opt["tooltip"]["axisPointer"] = Dict("type" => "shadow")
    opt["tooltip"]["trigger"] = "axis"
    opt["yAxis"] = [
        Dict("splitArea" => Dict("show" => true),
        "type" => "value")
        ]
    opt["xAxis"] = [
        Dict("type" => "category")
        ]

    opt["legend"]["data"] = legendlist

    len = length(data)

    lastsum_data = zeros(N, len)
    pos_data = Array(String, len)
    neg_data = Array(String, len)

    for i = 1:len
        if sign(data[i]) >= 0
            pos_data[i] = string(data[i])
            neg_data[i] = "-"
        else
            pos_data[i] = "-"
            neg_data[i] = string(-data[i])
        end
    end

    for i = 2:len
        if data[i-1] >=0
            lastsum_data[i] = lastsum_data[i-1] + data[i-1]
        else
            lastsum_data[i] = lastsum_data[i-1]
        end

        if data[i]< 0
            lastsum_data[i] = lastsum_data[i] + data[i]
        end
    end

    resize!(opt["series"], 3)

    opt["series"][1] = Dict(
        "name" => "",
        "type" => "bar",
        "smooth" => smooth,
        "stack" => "total",
        "data" => lastsum_data[1:end],
        "itemStyle" => Dict("normal" => Dict("borderColor"=>"rgba(0,0,0,0)","color"=>"rgba(0,0,0,0)"),
                        "emphasis" => Dict("borderColor"=>"rgba(0,0,0,0)","color"=>"rgba(0,0,0,0)"))
        )

    opt["series"][2] = Dict(
        "name" => legendlist[1],
        "type" => "bar",
        "smooth" => smooth,
        "stack" => "total",
        "data" => pos_data,
        "itemStyle" => Dict("normal" => Dict("label" => Dict("show" => true, "position" => "inside")))
        )

    opt["series"][3] = Dict(
        "name" => legendlist[2],
        "type" => "bar",
        "smooth" => smooth,
        "stack" => "total",
        "data" => neg_data,
        "itemStyle" => Dict("normal" => Dict("label" => Dict("show" => true, "position" => "inside")))
        )

    merge!(opt["xAxis"][1],Dict("data" => labels))

    return opt
end


text:
function template()
    @eval begin
        using PkgTemplates

        Template(;
        user="Ollecram-Friends",    # Take default from github.user & github.email
        # dir=                      # Take default from $JULIA_PKG_DEVDIR iff defined
        # authors=                  # Build default from `user` argument or default
        julia=v"1.5",               # Require Julia version 1.5.x

        plugins=[           # Plugins' start

            ProjectFile(;           # Creates a Project.toml
            version=v"0.1.0"                # The initial version of created packages.
            ),

            SrcDir(;                # Creates a module entrypoint
            #file=                          # Template file for src/<module>.jl
            ),

            Tests(;                 # Sets up testing for packages
            #file=                          # Template file for runtests.jl
            project=true                    # Create new project for tests (test/Project.toml)
            ),

            Readme(;                # Creates README file with badges for other included plugins
            # file=                         # Template file for the README
            # destination=                  # File destination, rel/to the repository root
            inline_badges=false,            # Badges on the same line as the package name?
            ),

            License(;               # Creates a license file
            name="ASL",                     # Apache Software License
            path=nothing,                   # Path to a custom license file.
            # destination=                  # File destination, rel/to the repository root
            ),

            Git(;                   # Creates a Git repository and a .gitignore file
            ignore=["*.jl.cov", "*.jl.*.cov", "*.jl.mem",
                    "deps/deps.jl", "deps/build.log", "deps/downloads/", "deps/usr/", "deps/src/",
                    "docs/build/", "docs/site/",
                    "Manifest.toml"],
            # name=                         # Take default from git config variable user.name
            # email=                        # Take default from git config variable user.email
            branch="main",                  # Desired name of the repos's default branch
            ssh=true,                       # Use SSH (vs HTTPS)
            jl=true,                        # Add ".jl" to the git repository URL
            manifest=false,                 # Never commit Manifest.toml
            gpgsign=false,                  # No GPG signing of Git commits
            ),

            CompatHelper(;          # Integrates with CompatHelper via GitHub Actions
            # file=                             # Template file for the workflow file
            destination="CompatHelper.yml",     # Destination rel/to .github/workflows
            cron="0 0 * * *",                   # Cron schedule interval: "0 0 * * *" means weekly, on Sunday
            ),

            TagBot(;                # Adds GitHub release support via TagBot
            # file=                             # Template file for the workflow file.
            destination="TagBot.yml",           # Destination rel/to .github/workflows
            trigger="JuliaTagBot",              # Trigger username for custom registries
            token=Secret("GITHUB_TOKEN"),       # Name of the token secret to use
            ssh=Secret("DOCUMENTER_KEY"),       # Name of the SSH private key secret to use
            ssh_password=nothing,               # Name of the SSH key password secret to use
            changelog=nothing,                  # Custom changelog template
            changelog_ignore=nothing,           # Issue/pull request labels to ignore in changelog
            gpg=nothing,                        # Name of the GPG private key secret to use
            gpg_password=nothing,               # Name of the GPG private key password secret to use
            registry=nothing,                   # Custom registry, in the format owner/repo
            branches=nothing,                   # Whether not to enable the branches option
            dispatch=nothing,                   # Whether or not to enable the dispatch option
            dispatch_delay=nothing,             # Number of minutes to delay for dispatch events
            ),

            GitHubActions(;         # Integrates your packages with GitHub Actions
            # file=                             # Template file for the workflow file
            destination="ci.yml",               # Destination rel/to .github/workflows
            linux=true,                         # Whether or not to run builds on Linux
            osx=false,                          # Whether or not to run builds on OSX (Mac OS)
            windows=false,                      # Whether or not to run builds on Windows
            x64=true,                           # Whether or not to run builds on 64-bit architecture
            x86=false,                          # Whether or not to run builds on 32-bit architecture
            coverage=true,                      # Whether or not to publish code coverage
            extra_versions=[v"1.2"],            # Extra Julia versions to test
            ),

            Codecov(;           # Sets up code coverage submission from CI to Codecov
            ),

            Documenter{GitHubActions}(; # Sets up documentation generation via Documenter.jl
            # make_jl=                          # Template file for make.jl
            # index_md=                         # Template file for index.md
            assets=String[],                    # Extra assets for the generated site
            logo=Logo(),                        # A Logo w documentation logo information
            canonical_url=nothing,              # A function to generate the site's canonical URL
            makedocs_kwargs=Dict{Symbol, Any}(),    # Extra keyword args to be inserted into makedocs
            ),

            # Develop(),          # Adds generated packages to the current environment

        ],      # Plugins' end

    )       # Template's end
    end
end

max_stars_repo_name: Ollecram-Friends/JuliaPkgTemplates

max_issues_repo_name: Ollecram-Friends/JuliaPkgTemplates

max_forks_repo_name: Ollecram-Friends/JuliaPkgTemplates


text:
using Test, Hijack

# called with lazy=true or lazy=:brutal
# @test_* are ignored at toplevel or in for/if/let/begin

@test false
@test_throws ArgumentError 1==1

for _=1:1
    let cond=true
        let
            if cond
                begin
                    @test false
                end
                while true
                    @testset "lazy inner" begin
                        @test true
                        push!(Hijack.RUN, 1)
                    end
                    break
                end
            else
                @test false
                @testset "lazy inner 2" begin
                    @test true
                    push!(Hijack.RUN, 2)
                end
            end
        end
    end
end

@testset "lazy toplevel" begin
    @test true
    push!(Hijack.RUN, 3)
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
include("tsne_imports.jl")

get_output_path(path_to_output) = let
  joinpath(splitdir(path_to_output)[1], "output.bson")
end

# results_dir = "$(ENV["SCRATCH_DIR"])/cytof/results/repfam/test-sim-6-8-6"
results_dir = "../tex/results"
csvdir = "viz/csv"
mkpath(csvdir)

function get_best_lam(samples; true_Z=nothing, I=nothing)

  if typeof(samples) == String
    path_to_results = samples
    best_lam = [Int.(readdlm("$(path_to_results)/lam$(i)_best.txt"))
                for i in 1:I]
    best_Z = [Int.(readdlm("$(path_to_results)/Z$(i)_best.txt"))
              for i in 1:I]
  else
    lams = [s[:theta__lam] for s in samples]
    Ws = [s[:theta__W] for s in samples]
    Zs = [s[:theta__Z] for s in samples]
    I = length(lams[1])

    best_idx = [PlotUtils.estimate_ZWi_index(Zs, Ws, i) for i in 1:I]
    best_lam = [lams[best_idx[i]][i] for i in 1:I]
    best_Z = [Zs[best_idx[i]] for i in 1:I]
  end

  if true_Z == nothing
    return best_lam
  else
    population = Population()
    foreach(z -> label(population, z), eachcol(true_Z))
    relabeled_lam = [[label(population, best_Z[i][:, lam_in])
                      for lam_in in best_lam[i]] for i in 1:I]

    return relabeled_lam
  end
end

for zind in 1:3
  for pmiss in (0.0, 0.2)
    for phi in (0, 1, 10, 25, 100)
      # Get simulation name
      simname = "pmiss$(pmiss)-phi$(phi)-zind$(zind)"
      println("Processing: ", simname)

      # Get path to simuldate data.
      # path_to_output = joinpath(results_dir, simname, "output.bson")

      # Get output
      # output = BSON.load(path_to_output)

      # Get relevant parameters
      # samples =  [s[1] for s in output[:samples][1]]

      # Get ture Z
      true_Z = readdlm("../Z$(zind).txt", comments=true, comment_char='#')

      # Get best lam (for each sample)
      # best_lam = vcat(get_best_lam(samples, true_Z=true_Z)...)
      path_to_results = joinpath(results_dir, simname, "img/txt")
      best_lam = vcat(get_best_lam(path_to_results, true_Z=true_Z, I=2)...)

      # Write output to file.
      open("$(csvdir)/rfam-$(simname).csv", "w") do io
        writedlm(io, best_lam)
      end
    end
  end
end


text:
const sf = ReMat(ds[:Batch])
# check alternative constructors
@test sf == ReMat(ds[:Batch],ones(length(ds[:Batch])))

const Yield = convert(Array,ds[:Yield])

# check size methods
@test size(sf) == (30, 6)
@test size(sf,1) == 30
@test size(sf,2) == 6
@test size(sf,3) == 1

dd = fill(5., 6)
@test sf'ones(30) == dd

@test ones(30)'sf == dd'
const crp = sf'sf
@test isa(crp, Diagonal{Float64})
const crp1 = copy(crp)
@test crp1 == crp
@test crp[2,6] == 0
@test crp[6,6] == 5
@test size(crp) == (6,6)
@test crp.diag == fill(5.,6)
const rhs = sf'Yield
@test rhs == [7525.0,7640.0,7820.0,7490.0,8000.0,7350.0]
@test A_ldiv_B!(crp,copy(rhs)) == [1505.,1528.,1564.,1498.,1600.,1470.]

const L = ReTerms.LT(sf)
L[:θ] = [0.5]

@test isa(scale!(L,crp),Diagonal)
@test crp.diag == fill(2.5,6)
@test copy!(crp1,crp) == crp

const sf1 = ReMat(psts[:Sample])
const sf2 = ReMat(psts[:Batch])
@test size(sf1) == (60,30)
@test size(sf2) == (60,10)

const crp11 = sf1'sf1
const pr21 = sf2'sf1
const crp22 = sf2'sf2

@test isa(crp11,Diagonal{Float64})
@test isa(crp22,Diagonal{Float64})
@test isa(pr21,SparseMatrixCSC{Float64})