text:
export shrinkFactorGraph, getAllTrees, nnzFrontals, nnzClique, nnzTree, nnzSqrtInfoMatrix, getTreeCost_01, getTreeCost_02

"""
    $SIGNATURES

Prunes factor graph to keep up to `upto` number of variables.

Warning: uses functions that are outside of IncrementalInference (e.g.,
`isSolvable()`), so will probably need to place this elsewhere.
"""
function shrinkFactorGraph(fg; upto::Int=6)
  fgs = deepcopy(fg)

  delVars = filter(x->isSolvable(getVariable(fgs, x))==0, ls(fgs))
  todel = setdiff(lsf(fgs, solvable=0), lsf(fgs, solvable=1))
  delFcts = intersect(lsf(fgs), todel)
  allMags = filter(x->:MAGNETOMETER in getTags(getFactor(fgs, x)), lsfPriors(fgs) )
  union!(delFcts, filter(x->length(ls(fgs, x))==0, allMags) )

  union!(delVars, (ls(fgs, r"x\d") |> sortDFG)[upto:end])
  union!(delFcts, map(x->ls(fgs, x), delVars)...)

  map(x->deleteFactor!(fgs, x), delFcts)
  map(x->deleteVariable!(fgs, x), delVars)

  fgs
end

"""
    $SIGNATURES

Deterministically get all trees associated with all possible variable orderings
of `dfg` factor graph. Returns a dictionary with (tree, ordering, nnz) tuples.

Warning: factorial number of possibilities, so use carefully!
"""
function getAllTrees(fg::AbstractDFG)
    # dfg = generateCanonicalFG_Kaess(graphinit=false)
    variables = ls(fg)
    orderings = permutations(variables) |> collect

    # Dimensionality check to make sure we do not break your computer.
    max_dimension = 11 # something reasonable (11! ~ 40M).
    if length(variables) > 11
        throw(ArgumentError("You crazy! dfg is too big. Factorial explosion."))
    end

    # Produce a tree for each ordering, and store in dictionary.
    all_trees = Dict{Int, Tuple{BayesTree, Vector{Symbol}, Float64}}()
    for i in 1:length(orderings)
        tree = buildTreeReset!(fg, orderings[i])
        nnz = nnzTree(tree)
        all_trees[i] = (tree, orderings[i], nnz)
    end

    return all_trees
end

"""
    $SIGNATURES

Get number of non-zero entries for clique's frontal components. Num of non-zero
matrix entries is just the fully dense upper triangular part of square matrix.
"""
function nnzFrontals(dimension)
    if dimension == 1
        return 1
    else
        # Solved recurrence for n + (n-1) + ... + 2 + 1.
        return (dimension * (dimension + 1)) / 2
    end
end

"""
    $SIGNATURES

Get total number of non-zero entries for a clique. Num of non-zero matrix
entries is the fully dense upper triangular part (frontals) plus the
(frontal x separator)-size rectangle.
"""
function nnzClique(clique)
    frontal_dim = length(getCliqFrontalVarIds(clique))
    separator_dim = length(getCliqSeparatorVarIds(clique))
    return nnzFrontals(frontal_dim) + (frontal_dim * separator_dim)
end

"""
    $SIGNATURES

Get total number of non-zero entries for a Bayes tree. Num of non-zero matrix
entries is the sum of all non-zero entries for each individual clique.
"""
function nnzTree(tree::AbstractBayesTree)
  nnzTot = 0
  for (cliqid, cliq) in tree.cliques
    nnzTot += nnzClique(cliq)
  end
  return nnzTot
end

"""
    $SIGNATURES

Get total number of non-zero entries for a factor graph's upper triangular
square root information matrix, i.e., R matrix in A = Q[R 0]^T, using the QR's
factorization algorithm variable ordering.
"""
function nnzSqrtInfoMatrix(A::Matrix)
  q,r,p = qr(A, Val(true))
  r .= abs.(r)
  nz = 1e-5 .< r
  r[nz] .= 1
  sum(nz)
end

"""
    $SIGNATURES

Simple cost function for ranking the structure of a Bayes tree. Weighting:
    cost = (max tree depth) * (max clique dimension)^alpha
"""
function getTreeCost_01(tree::AbstractBayesTree; alpha::Float64=1.0 )
  cliqs = tree.cliques |> values |> collect
  maxdepth = map(x->getCliqDepth(tree, x)+1, cliqs) |> maximum
  maxdim = length.(map(x->getCliqVarIdsAll(x), cliqs)) |> maximum

  return maxdepth * (maxdim^alpha)
end

"""
    $SIGNATURES

Cost function for ranking the structure of a Bayes tree, putting and emphasis on
wider but shallower trees by penalizing the average number of siblings.
Weighting:
    cost = 1/(total num of child / num of parents) *
                           (max tree depth) * (max clique dimension)^alpha
"""
function getTreeCost_02(tree::AbstractBayesTree; alpha::Float64=1.0)
  # Frontal and number of children.
  ARR = Tuple{Symbol, Int}[]
  for (cliqid, vex) in tree.cliques
    afrtl = getCliqFrontalVarIds(tree.cliques[cliqid])[1]
    numch = length(getChildren(tree, afrtl))
    push!(ARR, (afrtl, numch))
  end

  numParents = filter(x->0<x[2], ARR) |> length
  totalNumChildren = (x->x[2]).(ARR) |> sum

  return getTreeCost_01(tree, alpha=alpha)/(totalNumChildren/numParents)
end

## ============================================================================
# Tools for checking the numerical performance of the solve

function mmdSolveKey(vari::DFGVariable,
                    refKey::Symbol,
                    tstKey::Symbol;
                    bw::AbstractVector{<:Real}=[0.001;] )
  #
  refVal = getBelief(vari, refKey)
  tstVal = getBelief(vari, tstKey)

  # calc mmd distance
  mmd(refVal, tstVal, getVariableType(vari), bw=bw)
end

# vari = getVariable(fg, :x1)
# kys = filter(x->!(x in [:graphinit;:default]), listSolveKeys(fg) |> collect |>sortDFG)
# X1_dist_0 = kys .|> x->mmdSolveKey(vari, :default_0, x)

# kyD = [(kys[i],kys[i+1]) for i in 1:length(kys)-1]
# X1_dist_D = kyD .|> x->mmdSolveKey(vari, x[1], x[2])

#


text:
### Implementation of observation model

function POMDPs.observation(pomdp::SingleOIPOMDP, sp::SingleOIState)
   # return SingleOIDistribution([1.0], [sp])
    if !is_observable_fixed(sp, pomdp.env) || off_the_grid(pomdp, sp.car)
        o = SingleOIObs(false, sp.ego, get_off_the_grid(pomdp))
        return SingleOIDistribution([1.0], [o])
    elseif is_crash(pomdp, sp.ego, sp.car)
        return SingleOIDistribution([1.0], [sp])
    end
    ego = sp.ego
    car = sp.car
    lane = pomdp.env.roadway[car.posF.roadind.tag]
    neighbors = Vector{VehicleState}(undef, 9)
    neighbors[1] = car_state(pomdp, lane, car.posF.s - pomdp.pos_res, car.v)
    neighbors[2] = car_state(pomdp, lane, car.posF.s + pomdp.pos_res, car.v)
    neighbors[3] = car_state(pomdp, lane, car.posF.s - pomdp.pos_res, car.v - pomdp.vel_res)
    neighbors[4] = car_state(pomdp, lane, car.posF.s + pomdp.pos_res, car.v - pomdp.vel_res)
    neighbors[5] = car_state(pomdp, lane, car.posF.s - pomdp.pos_res, car.v + pomdp.vel_res)
    neighbors[6] = car_state(pomdp, lane, car.posF.s + pomdp.pos_res, car.v + pomdp.vel_res)
    neighbors[7] = car_state(pomdp, lane, car.posF.s, car.v - pomdp.vel_res)
    neighbors[8] = car_state(pomdp, lane, car.posF.s, car.v + pomdp.vel_res)
    neighbors[9] = car_state(pomdp, lane, car.posF.s, car.v)

    states = SingleOIObs[]
    sizehint!(states, 9)
    for neighbor in neighbors
        if in_bounds_car(pomdp, neighbor) && !is_crash(pomdp, ego, neighbor)
            push!(states, SingleOIObs(false, ego, neighbor))
        end
    end
    probs = ones(length(states))
    probs = normalize!(probs, 1)
    return SingleOIDistribution(probs, states)
end
function POMDPs.observation(pomdp::SingleOIPOMDP, a::SingleOIAction, sp::SingleOIState)
    return observation(pomdp, sp)
end

"""
Check if a car in s is observable or not
"""
function is_observable_fixed(s::SingleOIState, env::SimpleInterEnv)
    m = length(env.obstacles)
    car = s.car
    ego = s.ego
    front = ego.posG + polar(VehicleDef().length/2, ego.posG.θ)
    angle = atan(car.posG.y - front.y, car.posG.x - front.x)
    ray = Projectile(VecSE2(front.x, front.y, angle), 1.0)
    for i = 1:m
        if is_colliding(ray, env.obstacles[i], car.posG)
            return false
        end
    end
    return true
end

"""
Check if a state is in bound
"""
function in_bounds_car(pomdp::SingleOIPOMDP, veh::VehicleState)
    lane = pomdp.env.roadway[veh.posF.roadind.tag]
    s = get_end(lane)
    return (veh.posF.s <= s) && (veh.v <= pomdp.env.params.speed_limit) && (veh.v >= 0.) && (veh.posF.s >= 0.)
end

ext: jl


text:
module KomePackage

function  greet(arg="Test")
	return "Hello $arg"

end

end # module


text:
const IndexTuple = Union{Tuple, NamedTuple}

"""
    IndexSpace(interval, boundingrect, nrows)

Metadata about an chunk.

- `interval`: An `Interval` object with the first and the last index tuples.
- `boundingrect`: An `Interval` object with the lowest and the highest indices as tuples.
- `nrows`: A `Nullable{Int}` of number of rows in the NDSparse, if knowable.
"""
struct IndexSpace{T<:IndexTuple}
    interval::Interval{T}
    boundingrect::Interval{T}
    nrows::Nullable{Int}
end

"""
A distributed table
"""
mutable struct DIndexedTable{T,K} <: AbstractIndexedTable
    # primary key columns
    pkey::Vector{Int}
    # extent of values in the pkeys
    domains::Vector{IndexSpace}
    chunks::Vector
end

noweakref(w::WeakRefString) = string(w)
noweakref(x) = x
function Dagger.domain(t::IndexedTable)
    ks = pkeys(t)
    T = eltype(ks)

    if isempty(t)
        return t.pkey => EmptySpace{T}()
    end

    wrap = T<:NamedTuple ? IndexedTables.namedtuple(fieldnames(T)...)∘tuple : tuple

    interval = Interval(map(noweakref, first(ks)), map(noweakref, last(ks)))
    cs = astuple(columns(ks))
    extr = map(extrema, cs[2:end]) # we use first and last value of first column
    boundingrect = Interval(wrap(noweakref(first(cs[1])), noweakref.(map(first, extr))...),
                            wrap(noweakref(last(cs[1])), noweakref.(map(last, extr))...))
    return t.pkey => IndexSpace(interval, boundingrect, Nullable{Int}(length(t)))
end

# if one of the input vectors is a Dagger operation / array
# chose the distributed implementation.
IndexedTables._impl(::Val, ::Dagger.ArrayOp, x...) = Val(:distributed)

function table(::Val{:distributed}, tup::Tup; chunks=nothing, kwargs...)
    if chunks === nothing
        # this means the vectors are distributed.
        # pick the first distributed vector and distribute
        # all others similarly
        idx = findfirst(x->isa(x, ArrayOp), tup)
        if idx == 0
            error("Don't know how to distribute. specify `chunks`")
        end
        darr = compute(tup[idx])
        chunks = domainchunks(darr)
    end

    darrays = map(x->distribute(x, chunks), tup)

    if isempty(darrays)
        error("Table must be constructed with at least one column")
    end

    nchunks = length(darrays[1].chunks)
    cs = Array{Any}(undef, nchunks)
    names = isa(tup, NamedTuple) ? [keys(tup)...] : nothing
    f = delayed((cs...) -> table(cs...; names=names, kwargs...))
    for i = 1:nchunks
        cs[i] = f(map(x->x.chunks[i], darrays)...)
    end
    fromchunks(cs)
end

# Copying constructor
function table(t::Union{IndexedTable, DIndexedTable};
               columns=columns(t),
               pkey=t.pkey,
               presorted=false,
               copy=true, kwargs...)

    table(columns;
          pkey=pkey,
          presorted=presorted,
          copy=copy, kwargs...)
end

Base.eltype(dt::DIndexedTable{T}) where {T} = T
function colnames(t::DIndexedTable{T}) where T
    fieldnames(T)
end

function trylength(t)::Nullable{Int}
    len = 0
    for l in map(x->x.nrows, t.domains)
        if !isnull(l)
            len = len + get(l)
        else
            return nothing
        end
    end
    return len
end

function Base.length(t::DIndexedTable)
    l = trylength(t)
    if isnull(l)
        error("The length of the DNDSparse is not yet known since some of its parts are not yet computed. Call `compute` to compute them, and then call `length` on the result of `compute`.")
    else
        get(l)
    end
end

"""
    fromchunks(cs)

Construct a distributed object from chunks. Calls `fromchunks(T, cs)`
where T is the type of the data in the first chunk. Computes any thunks.
"""
function fromchunks(cs::AbstractArray, args...; output=nothing, fnoffset=0, kwargs...)
    if output !== nothing
        if !isdir(output)
            mkdir(output)
        end
        cs = Any[(
            fn = lpad(idx+fnoffset, 5, "0");
            delayed(Dagger.savechunk, get_result=true)(
                c, output, fn);
        ) for (idx, c) in enumerate(cs)]

        vec_thunk = treereduce(delayed(vcat, meta=true), cs)
        if length(cs) == 1
            cs = [compute(get_context(), vec_thunk)]
        else
            cs = compute(get_context(), vec_thunk)
        end
        return fromchunks(cs)
    elseif any(x->isa(x, Thunk), cs)
        vec_thunk = delayed((refs...) -> [refs...]; meta=true)(cs...)
        cs = compute(get_context(), vec_thunk)
    end
    T = chunktype(first(cs))
    fromchunks(T, cs, args...; kwargs...)
end

function fromchunks(::Type{<:AbstractVector},
                    cs::AbstractArray, args...; kwargs...)
    lengths = length.(domain.(cs))
    dmnchunks = DomainBlocks((1,), (cumsum(lengths),))
    T = promote_eltype_chunktypes(cs)
    DArray(T, ArrayDomain(1:sum(lengths)),
           dmnchunks, cs, dvcat)
end

function fromchunks(::Type{<:IndexedTable}, chunks::AbstractArray;
                    domains::AbstractArray = last.(domain.(chunks)),
                    pkey=first(domain(first(chunks))),
                    K = promote_eltypes(domains),
                    T = promote_eltype_chunktypes(chunks))

    nzidxs = findall(!isempty, domains)
    domains = domains[nzidxs]

    DIndexedTable{T, K}(pkey, domains, chunks[nzidxs])
end

function promote_eltypes(ts::AbstractArray)
    nonempty_domains = Iterators.filter(!has_empty_values, ts)
    mapreduce(eltype, _promote_type, nonempty_domains, init = Union{})
end

function promote_eltype_chunktypes(ts::AbstractArray)
    nonempty_chunks = Iterators.filter(!has_empty_values∘domain, ts)
    mapreduce(eltype∘chunktype, _promote_type, nonempty_chunks, init = Union{})
end

"""
    distribute(t::Table, chunks)

Distribute a table in `chunks` pieces. Equivalent to `table(t, chunks=chunks)`.
"""
distribute(t::IndexedTable, chunks) = table(t, chunks=chunks, copy=false)

compute(t::DIndexedTable; kwargs...) = compute(get_context(), t; kwargs...)

function compute(ctx, t::DIndexedTable; output=nothing)
    if any(Dagger.istask, t.chunks)
        fromchunks(IndexedTable, cs, output=output)
    else
        map(Dagger.unrelease, t.chunks) # don't let this be freed
        foreach(Dagger.persist!, t.chunks)
        t
    end
end

distribute(t::DIndexedTable, cs) = table(t, chunks=cs)

collect(t::DIndexedTable) = collect(get_context(), t)

function collect(ctx::Context, dt::DIndexedTable{T}) where T
    cs = dt.chunks
    if length(cs) > 0
        collect(ctx, treereduce(delayed(_merge), cs))
    else
        table(similar(IndexedTables.arrayof(T), 0), pkey=dt.pkey)
    end
end

# merging two tables

function _merge(f, a::IndexedTable, b::IndexedTable)
    if isempty(a.pkey) && isempty(b.pkey)
        return table(vcat(rows(a), rows(b)))
    end

    @assert a.pkey == b.pkey
    ia = pkeys(a)
    ib = pkeys(b)

    if isempty(a)
        b
    elseif isempty(b)
        a
    elseif last(ia) < first(ib)
        # can vcat
        table(map(vcat, columns(a), columns(b)), pkey=a.pkey, presorted=true, copy=false)
    elseif last(ib) < first(ia)
        _merge(f, b, a)
    else
        f(a,b)
    end
end

_merge(f, x::IndexedTable) = x
function _merge(f, x::IndexedTable, y::IndexedTable, ys::IndexedTable...)
    treereduce((a,b)->_merge(f, a, b), [x,y,ys...])
end

_merge(x::IndexedTable, y::IndexedTable...) = _merge((a,b) -> merge(a, b), x, y...)

function Base.show(io::IO, big::DIndexedTable)
    h, w = displaysize(io)
    showrows = h - 5 # This will trigger an ellipsis when there's
                     # more to see than the screen fits
    t = first(Iterators.partition(big, showrows))
    len = trylength(big)
    vals = isnull(len) ? "" : " with $(get(len)) rows"
    header = "Distributed Table$vals in $(length(big.chunks)) chunks:"
    cstyle = Dict(i=>:bold for i in t.pkey)
    showtable(io, t; header=header, ellipsis=:end, cstyle=cstyle)
end


text:
type BlackCallableFixedRateBondEngine{V <: CallableBondVolatilityStructure} <: PricingEngine
  volatility::V
end

function BlackCallableFixedRateBondEngine(fwdYieldVol::Quote)
  vol = CallableBondConstantVol(0, NullCalendar(), fwdYieldVol, Actual365())

  return BlackCallableFixedRateBondEngine(vol)
end


text:
module Thermomodel

export zhang2002model!,dMdtzhang2002model,dynamicsmodel,dMdtdynamicsmodel,wallmodel,liquidmodel

using ..Systems,..Tools

"""
    This function is required by "DifferentialEquation.jl" Package.
        du :: an empty state vector to be derived at the end of this function.
        u  :: the state vector input from the last time step.
        p  :: some important parameters which do not belong to u
        t  :: the current time

    This approach solves the same set of dimensionless governing equations as
    "https://doi.org/10.1016/S0017-9310(01)00348-9"(Zhang et al. 2002)
    Instead of solving them piece by piece. Our current approach solve them as a whole system.
"""

function zhang2002model!(du::Array{Float64,1},u::Array{Float64,1},p::PHPSystem,t::Float64)

    (Xp,dXdt0,M)=vectoXM(u)

    numofliquidslug =  Integer( (length(u) - 1)/5 )
    sys0 = p

    γ = sys0.liquid.γ
    ω0 = sys0.liquid.ω0
    ℘ = sys0.liquid.℘
    Lvaporplug = XptoLvaporplug(Xp,sys0.tube.L,sys0.tube.closedornot)

    height = getheight(Xp,sys0.tube.L2D,sys0.tube.angle)
    Xpvapor = getXpvapor(Xp,sys0.tube.L,sys0.tube.closedornot)

    # get P from M and γ
    # P = (M./Lvaporplug).^(γ)
    # P = zeros(size(M))
    # for i in length(M)
    #        P[i] = M[i] > 0 ? (M[i]./Lvaporplug[i]).^(γ) : -(-M[i]./Lvaporplug[i]).^(γ)
    #    end
    P = real.((M./Lvaporplug .+ 0im).^(γ))

    # get θ from P and γ
    # θ = zeros(size(P))
    # for i in length(P)
    #        θ[i] = P[i] > 0 ? P[i].^((γ-1)/γ) : -(-P[i]).^((γ-1)/γ)
    #    end
    θ = real.((P .+ 0im).^((γ-1)/γ))
    # θ = P.^((γ-1)/γ)

    for i = 1:numofliquidslug
        du[2*i-1] = u[2*numofliquidslug+2*i-1]
        du[2*i] = du[2*i-1]

        du[2*numofliquidslug + 2*i-1] = -32*u[2*numofliquidslug + 2*i-1] - (ω0[i]^2)*(0.5*(height[i][end]-height[i][1])) + ℘[i]*(P[i]-P[i+1])
        du[2*numofliquidslug + 2*i] = du[2*numofliquidslug + 2*i-1]

    end

# not sure which one is better
        du[4*numofliquidslug+1:5*numofliquidslug+1] .= dMdtzhang2002model(Xpvapor,θ,sys0)
        # du[4*numofliquidslug+1:5*numofliquidslug+1] .= dMdtconstantH(Xpvapor,θ,sys0)

    return du

end

"""
    This function is a sub-function of zhang2002model. This function is the phase change part,
        Xpvapor ::  storing the location of two ends of each vapor
        θ       ::  the temperature in each vapor.
        sys0    ::  the struct that stores every needed initial conditions and boundary conditions

    In this case. He and Hc are preset constants.
"""

function dMdtzhang2002model(Xpvapor::Array{Tuple{Float64,Float64},1},θ::Array{Float64,1},sys0::PHPSystem)

    dMdt=zeros(length(Xpvapor))

    Xe = sys0.evaporator.Xe
    He = sys0.evaporator.He
    θe = sys0.evaporator.θe

    Xc = sys0.condenser.Xc
    Hc = sys0.condenser.Hc
    θc = sys0.condenser.θc

    Levapoverlap=XpvaportoLoverlap(Xpvapor,Xe)
    Lcondoverlap=XpvaportoLoverlap(Xpvapor,Xc)

    # May not be right for multi liquid flow
    for i = 1:length(Xpvapor)
        if Lcondoverlap[i] < 1e-8
            dMdt[i] = He * Levapoverlap[i] * (θe-θ[i])
        else
            dMdt[i] = -Hc * Lcondoverlap[i] * (θ[i]-θc)
        end
    end

    return dMdt

end

function dMdtconstantH(Xpvapor::Array{Tuple{Float64,Float64},1},θ::Array{Float64,1},sys0::PHPSystem)

    dMdt=zeros(length(Xpvapor))

    Xe = sys0.evaporator.Xe
    He = sys0.evaporator.He
    θe = sys0.evaporator.θe

    Xc = sys0.condenser.Xc
    Hc = sys0.condenser.Hc
    θc = sys0.condenser.θc

    Levapoverlap=XpvaportoLoverlap(Xpvapor,Xe)
    Lcondoverlap=XpvaportoLoverlap(Xpvapor,Xc)

    for i = 1:length(Xpvapor)
        dMdt[i] = He * Levapoverlap[i] * (θe-θ[i]) -Hc * Lcondoverlap[i] * (θ[i]-θc)
    end
    return dMdt

end

function dynamicsmodel(u::Array{Float64,1},p::PHPSystem)

if p.tube.closedornot == false

    du = zero(deepcopy(u))

    (Xp,dXdt0,M,δ)=vectoXMδ(u)

    numofliquidslug =  Integer( (length(u) - 2)/6 )
    sys = deepcopy(p)

    γ = sys.liquid.γ
    ω0 = sys.liquid.ω0
    ℘ = sys.liquid.℘
    Lvaporplug = XptoLvaporplug(Xp,sys.tube.L,sys.tube.closedornot)
    height = getheight(Xp,sys.tube.L2D,sys.tube.angle)
    Xpvapor = getXpvapor(Xp,sys.tube.L,sys.tube.closedornot)

    P = real.((M./Lvaporplug .+ 0im).^(γ))
    θ = real.((P .+ 0im).^((γ-1)/γ))

    for i = 1:numofliquidslug
        du[2*i-1] = u[2*numofliquidslug+2*i-1]
        du[2*i] = du[2*i-1]
        du[2*numofliquidslug + 2*i-1] = -32*u[2*numofliquidslug + 2*i-1] - (ω0[1]^2)*(0.5*(height[i][end]-height[i][1])) + ℘[1]*(P[i]-P[i+1])
        # use ℘[1] and ω0[1] for now, in the future change them to one value rather than a array.
        du[2*numofliquidslug + 2*i] = du[2*numofliquidslug + 2*i-1]

    end

# not sure which one is better
    du[4*numofliquidslug+1:5*numofliquidslug+1] .= dMdtdynamicsmodel(Xpvapor,θ,sys)
    du[5*numofliquidslug+2:end] .= [0.0]

    return du
end

if p.tube.closedornot == true

        du = zero(deepcopy(u))

        (Xp,dXdt0,M,δ)=vectoXMδ(u)

        numofliquidslug =  Integer( (length(u))/6 )
        sys = deepcopy(p)

        γ = sys.liquid.γ
        ω0 = sys.liquid.ω0
        ℘ = sys.liquid.℘
        Lvaporplug = XptoLvaporplug(Xp,sys.tube.L,sys.tube.closedornot)
        height = getheight(Xp,sys.tube.L2D,sys.tube.angle)
        Xpvapor = getXpvapor(Xp,sys.tube.L,sys.tube.closedornot)

        P = real.((M./Lvaporplug .+ 0im).^(γ))
        θ = real.((P .+ 0im).^((γ-1)/γ))

        for i = 1:numofliquidslug
            du[2*i-1] = u[2*numofliquidslug+2*i-1]
            du[2*i] = du[2*i-1]

            if i != numofliquidslug
                du[2*numofliquidslug + 2*i-1] = -32*u[2*numofliquidslug + 2*i-1] - (ω0[1]^2)*(0.5*(height[i][end]-height[i][1])) + ℘[1]*(P[i]-P[i+1])
            else
                du[2*numofliquidslug + 2*i-1] = -32*u[2*numofliquidslug + 2*i-1] - (ω0[1]^2)*(0.5*(height[i][end]-height[i][1])) + ℘[1]*(P[i]-P[1])
            end
            # use ℘[1] and ω0[1] for now, in the future change them to one value rather than a array.
            du[2*numofliquidslug + 2*i] = du[2*numofliquidslug + 2*i-1]

        end

    # not sure which one is better
        du[4*numofliquidslug+1:5*numofliquidslug] .= dMdtdynamicsmodel(Xpvapor,θ,sys)
        du[5*numofliquidslug+1:end] .= [0.0]

        return du
end

end

function dMdtdynamicsmodel(Xpvapor::Array{Tuple{Float64,Float64},1},θ::Array{Float64,1},sys::PHPSystem)

    dMdt=zeros(length(Xpvapor))

    He = sys.evaporator.He

    Hc = sys.condenser.Hc

        dx = sys.wall.Xarray[2]-sys.wall.Xarray[1]

    for i = 1:length(Xpvapor)
        indexes = findall( x -> (mod(x[1],length(Xpvapor)) == mod(i,length(Xpvapor)) && (x[end] == -1)), sys.mapping.walltoliquid)

            for j in length(indexes)

                # print("index",indexes,"\n")
                # print("θarray",sys.wall.θarray,"\n")

                dMdt[i] += He*dx*(sys.wall.θarray[indexes[j]] - θ[i])

            end

    end

    return dMdt

end

function wallmodel(θarray::Array{Float64,1},p::PHPSystem)
    sys = deepcopy(p)

    du = zero(deepcopy(θarray))

    γ = sys.vapor.γ
    Hₗ = sys.liquid.Hₗ
    He = sys.evaporator.He
    dx = sys.wall.Xarray[2]-sys.wall.Xarray[1]

    Xarray = sys.wall.Xarray
    Wearray = getwallWearray(Xarray,sys)

    Hwc = sys.condenser.Hwc
    θc  = sys.condenser.θc
    Xc  = sys.condenser.Xc
    hevisidec=ifamong.(Xarray,[Xc])

    H = zero(deepcopy(θarray))
    θarray_temp_flow = zero(deepcopy(θarray))

        for i = 1:length(θarray)

            index = sys.mapping.walltoliquid[i]

            if index[2] == -1

                if sys.tube.closedornot == false
                    P = sys.vapor.P[index[1]]
                end

                if sys.tube.closedornot == true
                    if index[1] > length(sys.vapor.P)
                        P = sys.vapor.P[index[1]-length(sys.vapor.P)]
                    else
                        P = sys.vapor.P[index[1]]
                    end
                end

                θarray_temp_flow[i] = real.((P .+ 0im).^((γ-1)/γ))

                H = He
            else
                θliquidarrays = sys.liquid.θarrays
                θarray_temp_flow[i] = θliquidarrays[index[1]][index[2]]

                H = Hₗ
            end
        end

        du = sys.wall.α .* laplacian(θarray,sys.tube.closedornot) ./ dx ./ dx + H .* (θarray_temp_flow - θarray) .* dx + Wearray .* dx + hevisidec .* Hwc .* (θc .- θarray) .* dx
        return du
end

function liquidmodel(θarrays,p::PHPSystem)
    sys = deepcopy(p)

    du = zero.(deepcopy(θarrays))

    γ = sys.vapor.γ
    Hₗ = sys.liquid.Hₗ

    θarray_temp_wall = zero.(deepcopy(θarrays))
    for i = 1:length(θarrays)

        dx = sys.wall.Xarray[2]-sys.wall.Xarray[1]

        indexes = sys.mapping.liquidtowall[i]

        for j = 1:length(indexes)
            θarray_temp_wall[i][j] = sys.wall.θarray[indexes[j]]
        end

        du[i] = sys.wall.α .* laplacian(θarrays[i]) ./ dx ./ dx + Hₗ .* (θarray_temp_wall[i] - θarrays[i]) .* dx
    end

    return du
end

end


text:
using Underscores
using FunctionalCollections
using Utils

const Str = SubString
const Storage = Set{Str}

struct Problem
    ingredients::Storage
    at_least::Storage
end

function parse_file(f)
    splitter = " (contains "

    re = Problem[]
    for line in eachline(f)
        @assert occursin(splitter, line)
        ing, allergens = split(line, splitter)
        push!(re, Problem(Set(split(ing)), Set(split(allergens[1:end-1], ", "))))
    end
    re
end

function load()
    open(joinpath(@__DIR__, "21.txt"), "r") do f
        parse_file(f)
    end
end

function ex1(problem::Array{Problem})
    all_ing = Storage()
    all_all = Storage()
    for p in problem
        union!(all_ing, p.ingredients)
        union!(all_all, p.at_least)
    end

    unsafe = Storage()

    for allergen in all_all
        int = Storage()
        for p in problem
            allergen in p.at_least || continue
            combiner = isempty(int) ? union! : intersect!
            combiner(int, p.ingredients)
        end
        union!(unsafe, int)
    end

    safe = setdiff(all_ing, unsafe)

    sum(safe) do ing
        count(p -> ing in p.ingredients, problem)
    end
end

function ex2(problem::Array{Problem}, mapped::PersistentArrayMap{Str,Str}, allergens)::Union{Nothing,PersistentArrayMap{SubString,SubString}}
    length(mapped) == allergens && return mapped

    for p in problem
        for allergen in p.at_least
            allergen in keys(mapped) && continue
            # Now we have a candidate for mapping and find all possible ingredients.
            candidates = @_ filter(allergen::Str in _.at_least, problem) |> mapreduce(Set(_.ingredients), intersect, __)
            # But not the ones already mapped.
            setdiff!(candidates, values(mapped))

            for candidate in candidates
                sub = ex2(problem, assoc(mapped, allergen, candidate), allergens)
                !isnothing(sub) && return sub
            end
        end
    end
end

function ex2(problem::Array{Problem}, allergens=8)
    increasing_allergens = sort(problem, by=p -> length(p.at_least))
    mappings = ex2(increasing_allergens, PersistentArrayMap{Str,Str}(), allergens)
    @_ mappings |> keys |> collect |> sort |> map(mappings[_], __) |> join(__, ',')
end

problem = load()

println("Appearance count: ", ex1(problem))
println("Canonical dangerous ingredient list: ", ex2(problem))

using Test

@testset "Allergen Assessment" begin
    @testset "example 1" begin
        input1 = """
        mxmxvkd kfcds sqjhc nhms (contains dairy, fish)
        trh fvjkl sbzzf mxmxvkd (contains dairy)
        sqjhc fvjkl (contains soy)
        sqjhc mxmxvkd sbzzf (contains fish)
        """

        example1 = parse_file(IOBuffer(input1))

        @test ex1(example1) == 5
        @test ex2(example1, 3) == "mxmxvkd,sqjhc,fvjkl"
    end

    @testset "results" begin
        @test ex1(problem) == 2556
        @test ex2(problem) == "vcckp,hjz,nhvprqb,jhtfzk,mgkhhc,qbgbmc,bzcrknb,zmh"
    end
end

max_stars_repo_name: CmdQ/AoC


text:
#=
Template Matching Spike Sorting
=#

function slider_release_template(widget::Ptr,param_tuple,user_data::Tuple{Gui_Handles})

    han, = user_data

    han.sort_widgets.slider_active = false

    nothing
end

function template_slider(widget::Ptr,user_data::Tuple{Gui_Handles})

    han, = user_data

    if !han.sort_widgets.slider_active
        han.sort_widgets.slider_active = true
    end

    myval=getproperty(han.sc.adj_sort, :value, Float64) # primary display

    clus=han.buf.selected_clus

    if clus>0

        han.sc.temp.tol[clus] = myval

        if ((han.buf.count>0)&(han.sc.pause))
            SpikeSorting.template_cluster(han.sc,clus,han.sc.temp.templates[:,clus],han.sc.temp.sig_min[:,clus],han.sc.temp.sig_max[:,clus],han.sc.temp.tol[clus])
            SpikeSorting.replot_all_spikes(han.sc)

            #incremental plot
        end
    end
    han.c_changed=true
    nothing
end

function draw_templates_c3(sc::SpikeSorting.Single_Channel)

    ctx=Gtk.getgc(sc.c3)

    myheight=height(ctx)
    mywidth=width(ctx)

    total_clus = max(sc.total_clus+1,5)

    for clus = 1:sc.total_clus

        s=sc.s * .25
        o=sc.o

        Cairo.translate(ctx,0.0,50.0)
        Gtk.scale(ctx,mywidth/(sc.wave_points*total_clus),s)

        startx=(clus-1)*(sc.wave_points)+1
        move_to(ctx,1.0+startx,sc.temp.templates[1,clus]-o)

        for i=2:size(sc.temp.sig_max,1)
            y=sc.temp.templates[i,clus]-o
            line_to(ctx,i+startx,y)
        end

        select_color(ctx,clus+1)
        set_line_width(ctx,1.0)
        stroke(ctx)

        identity_matrix(ctx)
    end

    nothing
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
# documentation files (the “Software”), to deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of
# the Software. THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT
# LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT
# SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
# CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.

using Raytracing, Raytracing.Interpreter
using Test, LinearAlgebra, StaticArrays
import ColorTypes:RGB

##########################################################################################92

@testset "test_RGB" begin
	include("test_RGB.jl")
end

@testset "test_HDRimage" begin
	include("test_HDRimage.jl")
end

@testset "test_ReadingWriting" begin
	include("test_ReadingWriting.jl")
end

@testset "test_ToneMapping" begin
	include("test_ToneMapping.jl")
end

##########################################################################################92

@testset "test_Geometry" begin
	include("test_Geometry.jl")
end

@testset "test_Transformation" begin
	include("test_Transformation.jl")
end

@testset "test_Rays-Cameras-ImageTracer" begin
	include("test_Rays-Cameras-ImageTracer.jl")

end

##########################################################################################92

@testset "test_AABB" begin
	include("test_AABB.jl")
end

@testset "test_Sphere" begin
	include("test_Sphere.jl")
end

@testset "test_Plane" begin
	include("test_Plane.jl")
end

@testset "test_Cube" begin
	include("test_Cube.jl")
end

@testset "test_Triangle" begin
	include("test_Triangle.jl")
end

@testset "test_Shapes-World" begin
	include("test_Shapes-World.jl")
end

##########################################################################################92

@testset "test_Pigment" begin
	include("test_Pigment.jl")
end

@testset "test_Renderers" begin
	include("test_Renderers.jl")
end

@testset "test_PCG" begin
	include("test_PCG.jl")
end

@testset "test_OrthoNormalBasis" begin
	include("test_OrthoNormalBasis.jl")
end

@testset "test_Antialiasing" begin
	include("test_antialiasing.jl")
end

##########################################################################################92

@testset "test_RangeTesters" begin
	include("test_RangeTesters.jl")
end

@testset "test_Interpreter" begin
	include("test_Interpreter.jl")
end

@testset "test_Demo-DemoAnimation.jl" begin
	include("test_Demo-DemoAnimation.jl")
end

@testset "test_Render-RenderAnimation.jl" begin
	include("test_Render-RenderAnimation.jl")
end

@testset "test_Raytracer.jl" begin
	include("test_Raytracer.jl")
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

using LinearAlgebra
# using ProgressMeter
# using Plots

# +
import Base.Threads: @spawn, @sync, nthreads, threadid
include("../src/mpiinterface.jl")
# include("../src/nompiinterface.jl")
include("../src/advection.jl")
include("../src/rotation.jl")
include("../src/spline.jl")
include("../src/bspline.jl")
include("../src/bsplinelu.jl")
include("../src/bsplinefft.jl")
include("../src/lagrange.jl")
include("../src/interpolation.jl")

using DoubleFloats

function printout(
    advd::Advection1dData{T,Nsp,Nv,Nsum,timeopt},
    str,
) where {T,Nsp,Nv,Nsum,timeopt}
    if timeopt != MPIOpt || advd.adv.mpid.ind == 1
        println(str)
    end
end

"""

   exact(tf, mesh1, mesh2)

   Julia function to compute exact solution "

```math
\frac{d f}{dt} +  (y \frac{df}{delta1} - x \frac{df}{delta2}) = 0
```

"""
function exact!(f, mesh1::UniformMesh{T}, mesh2::UniformMesh{T}, tf::T) where {T}

    for (i, x) in enumerate(mesh1.points), (j, y) in enumerate(mesh2.points)
        xn = cos(tf) * x - sin(tf) * y
        yn = sin(tf) * x + cos(tf) * y
        f[i, j] = exp(-13 * ((xn)^2 + (yn + T(6 // 5))^2))
    end
    f
end

function trace_diffrotation(
    advd::Advection1dData{T,Nsp,Nv,Nsum,timeopt},
    t::T,
    n,
) where {T,Nsp,Nv,Nsum,timeopt}

    if t == 0
        printout(advd, "#nt\ttime\tnorm_2\tnorm_inf")
    end
    ptr = pointer(advd.bufdata)
    f = unsafe_wrap(Array, ptr, sizeall(advd.adv))
    exact!(f, advd.adv.t_mesh_sp[1], advd.adv.t_mesh_v[1], t)
    n2 = Float64(norm(advd.data - f))
    ninf = Float64(norm(advd.data - f, Inf))
    printout(advd, "$n\t$(Float32(t))\t$n2\t$ninf")

end

function rotation(advd::Advection1dData, nbdt)

    # global cl_obs
    # clockreset(cl_obs)

    dt = advd.adv.dt_base
    trace_diffrotation(advd, zero(dt), 0)
    for i = 1:nbdt
        while advection!(advd)
        end
        trace_diffrotation(advd, i * dt, i)
    end
    println("#  end")
    # printall(cl_obs)
end
function rotation1_1(T::DataType, nbdt, timeopt; sz = (64, 64), interp = Lagrange(T, 21))
    dt = T(2big(pi) / nbdt)

    spmin, spmax, nsp = T(-5), T(5), sz[1]
    vmin, vmax, nv = -T(5), T(5), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    adv = Advection1d(
        (mesh_sp,),
        (mesh_v,),
        (interp,),
        (interp,),
        dt,
        timeopt = timeopt,
        tab_fct = [tan, sin, tan],
    )

    data = zeros(T, sizeall(adv))

    exact!(data, mesh_sp, mesh_v, zero(T))

    pvar = getrotationvar(adv)

    advd = Advection1dData(adv, data, pvar)

    printout(advd, "# dt=$(Float64(dt)) nbdt=$nbdt size_x=$nsp size_v=$nv")
    printout(advd, "# sp : from $(Float64(mesh_sp.start)) to $(Float64(mesh_sp.stop))")
    printout(advd, "# v : from $(Float64(mesh_v.start)) to $(Float64(mesh_v.stop))")
    printout(advd, "# interpolation : $(get_type(interp)) order=$(get_order(interp))")
    printout(advd, "# type=$T precision = $(precision(T))")
    printout(advd, "# timeopt=$timeopt")
    if timeopt == SimpleThreadsOpt || timeopt == SplitThreadsOpt
        printout(advd, "# nb threads : $(Threads.nthreads())")
    elseif timeopt == MPIOpt
        printout(advd, "# nb process : $(adv.mpid.nb)")
    else
        printout(advd, "# monothread version")
    end
    printout(advd, "typeof(data)=$(typeof(data)) size(data)=$(size(data))")

    # advdata = Advection1dData(adv, data, pvar)

    rotation(advd, nbdt)
end

T = Double64
nbdt = 1000
dt = T(2big(pi)) / nbdt
# landau1_1(T, 50, NoTimeOpt, sz=(64,128))
#rotation1_1(T, nbdt, MPIOpt, sz=(256,256), interp=Lagrange(T, 101))
rotation1_1(T, nbdt, MPIOpt, sz = (256, 256), interp = Lagrange(T, 5))


text:
#=
Written by Rasmus Kyng.
Based closely on Matlab code by Anup Rao, Sushant Sachdeva & Rasmus Kyng.
Matlab code available at https://github.com/sachdevasushant/Isotonic

Interior point method for computing the Isotonic regression (under Euclidean norm) on arbitrary directed acyclic graphs. The code accompanies the paper 'Fast, Provable Algorithms for Isotonic Regression in all l_p-norms' at NIPS 2015

This repository presently contains code for computing isotonic regression in directed acyclic graphs (DAG).

isotonicIPM takes as input the adjacency matrix ‘A’ of a DAG, followed by a vector ‘v’ giving the initial values on the vertices. v should be a vector of length equal to the number of vertices.

On input ‘A’ and ‘v’, the code computes a vector ‘x’ of same length as v, such that x minimizes the quantity: sum_i (x(i) - v(i))^2 subject to x(i) <= x(j) if (i,j) is an entry in ‘a’.

=#

"""Compute isotonic regression of v with constraints given by A,
and multiplicative error param eps,
using a PD SDD linear system solver, given by argument *solver*."""
function isotonicIPMrelEps(A::SparseMatrixCSC{Tv,Ti},
                           v::Array{Tv,1},
                           eps::Real=0.1,
                           solver=(H -> augTreeSddm(H,tol=1e-1,maxits=1000))) where {Tv,Ti}
    return isotonicIPM(A,v,eps,solver,true)
end

"""Compute isotonic regression of v with constraints given by A,
and additive error param eps,
using a PD SDD linear system solver, given by argument *solver*.
Setting relGapTermination=true means the IPM continues until a multiplicative error of eps is reached."""
function isotonicIPM(A::SparseMatrixCSC{Tv,Ti},
                     v::Array{Tv,1},
                     eps::Tv=0.1,
                     solver=(H -> augTreeSddm(H,tol=1e-1,maxits=1000)),
                     relGapTermination::Bool=false) where {Tv,Ti}
    n = size(A)[1]
    topoOrder = toposort(A)
    permMat = (sparse(I,n,n))[1:n,topoOrder]
    Bt = -dirEdgeVertexMat(A) #sign convention implies all edges should get positive values
    m = size(Bt)[1]

    Theta = m
    eps0 = eps
    beta2 = 10.0

    x = collect(1.0:n) # our initial guess
    x = permMat*x #this seems to be the right direction for permMat

    if notFeasible(x, Bt)
        error("Input graph is not a DAG!")
    end

    eTaStartFn(y) = ( (y-v)'*Bt'*(1 ./ (Bt*y))/norm(y-v)^2 )[1]
    #This seems to be the formula for the minimizer?

    #TODO
    #eTaStartFn(y) = ( (y-v)'*Bt'*(1./(Bt*y))/norm(y-v) )[1]
    #This is the formula Anup used?
    #Seems to give better results?

    eTa1 = eTaStartFn(x)
    #@printf "eTa1 = %f\n" eTa1
    while eTa1<0
        x = 2*x
        eTa1 = eTaStartFn(x)
        #Note: y'*Bt'*(1./(Bt*y)) is always pos when y is isotonic,
        #We can show this loop will complete
        #println("doubling")
    end

    #@printf "eTa1 = %f\n" eTa1

    if eTa1 < 0
        error("bad initialization!")
    end
    eTa2 = Theta/eps0
    mu0 = eTa1

    relGap = Inf
    i=0;
    muTooSmall = mu0 < eTa2

    # numSteps = 0

    while muTooSmall
        F(y) = mu0*norm(y-v)^2 - sum(log.(Bt*y))
        centered = false
        while !centered
            # numSteps += 1
            # if numSteps % 30 == 1
            #     println("So far ", numSteps, " steps were made")
            # end

            (H,xNewton) = l2NewtonStep( Bt,x,v,mu0, solver )

            s = Bt*x
            gradF = 2*mu0*(x-v) - Bt'*(1 ./ s)
            x = backtrackLineSearch(F,gradF,Bt,xNewton,x)
            centMeasure = (xNewton'*H*xNewton)[1]
            centered = centMeasure < 10.0^-2
            i = i+1
        end

        dval = dualVal(Bt,x,v,mu0)
        relGap = (sum((x-v).^2) - dval)/sum((x-v).^2)
        if relGapTermination
            mu0 = mu0*beta2
            if relGap < eps
                break
            end
            #muTooSmall stays false: we only terminate if relGap is good
        else
            #this is the old termination condition:
            #stop if within 1 % relGap or if eta is too large?
            if relGap < 10.0^-2
                break
            end
            muTooSmall = mu0 < eTa2
            if muTooSmall
                mu0 = mu0*beta2
            end
        end
    end

    if notFeasible(x, Bt)
        error("output is not isotonic!")
    end

    return (x,relGap,i)
end

function l2NewtonStep( Bt::SparseMatrixCSC{Tv,Ti},
                       x::Array{Tv,1},
                       v::Array{Tv,1},
                       mu0::Real,
                       solver ) where {Tv,Ti}
    m = size(Bt)[1]
    n = size(Bt)[2]
    d1 = 2 * mu0 * sparse(I,n,n)
    s = Bt * x
    d2 = sparse(1:m, 1:m, s .^ -2, m, m)
    grad = 2 * mu0 * (x - v) - Bt'*(1 ./ s)
    H1 = sparse(Bt'*d2*Bt)

    H = H1 + d1

    invMinEntryH = maximum(abs.(s)) ^ 2

    #println("Maximum eigenvalue of the matrix we are computing: ", eigs(invMinEntryH * H;nev=1,which=:LM,tol=1e-1)[1][1])
    #println("Smallest eigenvalues of the matrix we are computing: ", eigs(invMinEntryH * H;nev=2,which=:SM,tol=1e-1)[1][1])
    F = solver(invMinEntryH * H)
    xNewton = -F(grad)*invMinEntryH

    return (H,xNewton)
end

function dualVal( Bt::SparseMatrixCSC{Tv,Ti},
                  x::Array{Tv,1},
                  v::Array{Tv,1},
                  mu0::Real) where {Tv,Ti}
    s = Bt*x # all pos entries
    f = 1/mu0 * 1 ./ s
    dval = (-1/4 * sum((Bt'*f).^2) - f'*Bt*v)[1]
    return dval
end

function backtrackLineSearch(F,
                      gradF::Array{Tv,1},
                      Bt::SparseMatrixCSC{Tv,Ti},
                      xNewton::Array{Tv,1},
                      x::Array{Tv,1}) where {Tv,Ti}
    #track line search as in Boyd's book
    A0 = 0.01;
    B0 = 0.5;

    t = 1;
    iter = 0;

    xt = x + t * xNewton;

    while (notFeasible(xt, Bt)) || (F(xt) > F(x) + A0 * t * (gradF'*xNewton)[1])
        t = B0 * t;
        iter = iter + 1;
        if iter >= 100
            error("Number of iterations in line search exceeded the limit")
        end
        xt = x + t * xNewton;
    end

    return xt
end

function notFeasible(x::Array{Tv,1},
              Bt::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    return minimum(Bt * x) <= 0 #not feasible if some entry is non-positive
end


text:
module HeresController
using Genie, Genie.Renderer, Genie.Requests, Genie.Router
# using JSON
using SearchLight
using GeoJSONEx
using JSON2

function index()
    rs = SearchLight.query("select id, name, ST_ASGeoJSON(lnglat) as lnglat from heres")
    # respond(JSON.json(to_featurecollection(rs, :lnglat)), "application/json")
    respond(JSON2.write(to_featurecollection(rs, :lnglat)), "application/json")
end

function show()
    id = @params(:id)
    rs = SearchLight.query("select id, name, ST_ASGeoJSON(lnglat) as lnglat from heres where id=$id")
    if size(rs, 1) == 1
        fs = to_features(rs, :lnglat)
        respond(JSON2.write(fs[1]), "application/json")
    else
        respond("""{"error": "No value with $id"}""", "application/json")
    end
end

function create()
    json = jsonpayload()
    name = json["name"]
    lnglat = json["lnglat"] |> JSON2.write
    SearchLight.query("insert into heres (name, lnglat) values (
        '$name',
        ST_SetSRID(ST_GeomFromGeoJSON('$lnglat'), 4326)
    ) RETURNING id")[1][1]
end

function delete()
    id = @params(:id)
    SearchLight.query("delete from heres where id=$id")
end

end

max_stars_repo_name: huaxk/rlservice

max_issues_repo_name: huaxk/rlservice

max_forks_repo_name: huaxk/rlservice


text:
function do_pass()
    # pass
end

function inline_pass()
    # pass
end

function indexing()::Int64
    sum = 0
    a::Array{Int64} = []
    for i in (0:10-1)
        push!(a, i)
        sum += a[i+1]
    end
    return sum
end

function infer_bool(code::Int64)::Bool
    return code in [1, 2, 4]
end

function show()
    a1 = 10
    b9 = 2
    b2 = 2
    @assert(b2 == b9)
    a2::Float64 = 2.1
    println(join([a2], " "))
    for i in (0:10-1)
        println(join([i], " "))
    end
    for i in (0:2:10-1)
        println(join([i], " "))
    end
    a3 = -(a1)
    a4 = (a3 + a1)
    println(join([a4], " "))
    sum1 = indexing()
    println(join([sum1], " "))
    a5 = [1, 2, 3]
    println(join([length(a5)], " "))
    a9::Array{String} = ["a", "b", "c", "d"]
    println(join([length(a9)], " "))
    a6 = Set([1, 2, 3, 4])
    println(join([length(a6)], " "))
    a7 = Dict("a" => 1, "b" => 2)
    println(join([length(a7)], " "))
    a8 = true
    if a8
        println(join(["true"], " "))
    else

        if a4 > 0
            println(join(["never get here"], " "))
        end
    end
    if a1 == 11
        println(join(["false"], " "))
    else

        println(join(["true"], " "))
    end
    if 1 != nothing
        println(join(["World is sane"], " "))
    end
    do_pass()
    inline_pass()
    s = "1    2"
    println(join([s], " "))
end

function main()
    show()
end

main()


text:
using OneClassActiveLearning
using SVDD
using Test
using Random
using JuMP, Ipopt
using MLDataUtils, MLLabelUtils, MLKernels
using ValueHistories
using PyCall
using Statistics
using LinearAlgebra
using Dates
using JSON, Unmarshal
using InteractiveUtils
using Distributions
using GaussianMixtures
using Serialization

TEST_SOLVER = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
TEST_DATA_FILE = joinpath(@__DIR__, "test.csv")
TEST_OUTPUT_FILE = joinpath(@__DIR__, "output.tmp")
TEST_DATA_NUM_DIMENSIONS, TEST_DATA_NUM_OBSERVATIONS = size(load_data(TEST_DATA_FILE)[1])

Random.seed!(0)

@testset "OneClassActiveLearning" begin
    include("QueryStrategies/pool_qs_test.jl")
    include("QueryStrategies/batch_qs_test.jl")
    include("QueryStrategies/query_synthesis_utils_test.jl")
    include("QueryStrategies/query_synthesis_optimization_test.jl")
    include("QueryStrategies/query_synthesis_test.jl")
    include("QueryStrategies/subspace_qs_test.jl")
    include("oracle/oracle_test.jl")
    include("data_util_test.jl")
    include("evaluate_test.jl")
    include("setup_test.jl")
    include("serialize_test.jl")
    include("result_test.jl")
end


text:
# This file contains a function to
# reset an instance to default values

"""
    reset!(Result::ResultStruct)

Reset an instance of [`ResultStruct`](@ref) to default values.

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.gen_example(file)
s.Gen(file)
s.gen!()
# Generate result data
s.Result.reset!()

# output

```
"""
function reset!(Result::ResultStruct)

    # Reset all elements
    Result.Δt = Result.q = Result.threshold = 0.0
    Result.t = Result.x = Result.X_FFT_ABS = []
    Result.ν = Result.D = Result.c = Result.cw = Result.Dw = []

    return nothing

end


text:
"""
    loadtrain(path::String; kwargs...)

load training set on `path` with `normalize` option in kwargs; output Dataset struture
"""
function loadtrain(path::String; kwargs...)
    data = CSV.read(path, DataFrame; normalizenames=true)
    return Dataset(data; kwargs...), Labels(data)
end

"""
    loadtest(path::String; kwargs...)

load testing set on `path` with `normalize` option in kwargs; output Dataset struture
"""
function loadtest(path::String; kwargs...)
    data = CSV.read(path, DataFrame; normalizenames=true)
    return Dataset(data; kwargs...)
end

"""
    savepredictions(path::String, labels::Labels)

save `labels` predicted by a classifier into file on `path` in correct format
"""
function savepredictions(path::String, labels::Labels)
    ids, survived = out(labels)
    df = DataFrame(PassengerId = ids, Survived = survived)
    CSV.write(path, df)
end


text:
"""
SAC.jl provides routines for dealing with SAC-formatted time series files,
including reading, writing, filtering, mean removal, rotating, and so on.
Sister library SACPlot.jl can be used for plotting.
"""
module SAC

import Dates
using LinearAlgebra
import Statistics: mean, covm, varm
import DSP, Dierckx
import Glob
import Base: ==, copy, getindex, setindex!, time, write
import AbstractFFTs: fft
import FFTW: rfft

export
    SACtr,
    add!,
    add,
    bandpass!,
    bp!,
    copy,
    cut!,
    differentiate!,
    divide!,
    divide,
    envelope!,
    fft,
    flip_component!,
    highpass!,
    hp!,
    integrate!,
    interpolate!,
    lp!,
    lowpass!,
    modify!,
    modify,
    multiply!,
    multiply,
    normalise!,
    normalise,
    read_wild,
    rmean!,
    rotate_through!,
    rotate_through,
    rotate_to_gcp!,
    rotate_to_gcp,
    rtrend!,
    sample,
    stalta,
    subtract!,
    subtract,
    taper!,
    time,
    tshift!,
    window_average,
    write

include("constants.jl")
include("types.jl")
include("great_circle.jl")
include("io.jl")
include("operations.jl")
include("filtering.jl")
include("util.jl")
include("statistics.jl")

# Submodules
include("Stack.jl")
using .Stack
export stack

# Build all copying routines
"""Dict with keys given by name of each function to have a copying version.
   Where an abbreviated version exists, that is given as the value; otherwise
   the value is `nothing`"""
const copying_funcs = Dict(
    :bandpass! => :bp!,
    :cut! => nothing,
    :differentiate! => nothing,
    :envelope! => nothing,
    :flip_component! => :flip!,
    :highpass! => :hp!,
    :integrate! => :int!,
    :interpolate! => nothing,
    :lowpass! => :lp!,
    :rmean! => nothing,
    :rtrend! => nothing,
    :taper! => nothing,
    :tshift! => nothing,
    )
for (name, abbrev) in copying_funcs
    new_name = Symbol(string(name)[1:end-1])
    @eval begin
        function ($new_name)(s::Union{SACtr,AbstractArray{SACtr}}, args...; kwargs...)
            s_new = deepcopy(s)
            $(name)(s_new, args...; kwargs...)
            s_new
        end
        @doc """
        $(@doc $name)
            $($new_name)(s::Union{SACtr,AbstractArray{SACtr}}, args...; kwargs...) -> s_new

        Copying version of `$($name)` which returns modified version(s) of the trace(s)
        in `s`, leaving the originals unaltered.
        """ $new_name
        export $new_name
    end
    if abbrev != nothing
        new_abbrev = Symbol(string(abbrev)[1:end-1])
        @eval begin
            const $new_abbrev = $new_name
            export $new_abbrev
        end
    end
end

# Build methods which can be used with chaining (`|>`)
"""Dict with keys given by name of each function to have a chaining version.
   Where an abbreviated version exists, that is given as the value; otherwise
   the value is nothing."""
const chaining_funcs = Dict(
    copying_funcs...,
    :normalise! => nothing
    )

for (name, abbrev) in copying_funcs
    new_name = Symbol(string(name)[1:end-1])
    @eval begin
        ($name)(args...; kwargs...) = s -> ($name)(s, args...; kwargs...)
        ($new_name)(args...; kwargs...) = s -> ($new_name)(s, args...; kwargs...)
    end
end

include("precompile.jl")

end # module SAC


text:
# parser.jl, Deniz Yuret, March 30, 2015
#
# Transition based parser based on:
#
# [GN13]  Goldberg, Yoav; Nivre, Joakim. Training Deterministic Parsers with Non-Deterministic Oracles. TACL 2013.
# [H13]   http://honnibal.wordpress.com/2013/12/18/a-simple-fast-algorithm-for-natural-language-dependency-parsing
# [KGS11] Kuhlmann, Marco, Carlos Gómez-Rodríguez, and Giorgio Satta. Dynamic programming algorithms for transition-based dependency parsers. ACL 2011
# [N03]   Nivre, Joakim. An efficient algorithm for projective dependency parsing. IWPT 2003.

# We define Parser{T,V} as a parametric type where T is :ArcHybrid,
# :ArcEager etc. and V is a variant symbol such as :GN13, :R1.  This
# allows specialization and inheritence of methods.  Thanks to
# julia-users members Toivo Henningsson and Simon Danisch for suggesting
# this design.

type Parser{T,V}
    nword::Int            # number of words in sentence
    ndeps::Int            # number of dependency labels (ROOT=1)
    nmove::Int            # number of possible moves (set by init!)
    wptr::Position        # index of first word in buffer
    sptr::Int             # index of last word (top) of stack
    stack::Pvec           # nword vector for stack of indices
    head::Pvec            # nword vector of heads
    deprel::Dvec          # nword vector of dependency labels

    function Parser(nword::Int,ndeps::Int)
        init!(new(nword,ndeps,0,1,0,Pzeros(nword),Pzeros(nword),Droots(nword)))
    end

    function Parser(s::Sentence)
        nword = length(s.word)
        ndeps = length(s.vocab.deprels)
        init!(new(nword,ndeps,0,1,0,Pzeros(nword),Pzeros(nword),Droots(nword)))
    end

end # Parser

# A parser provides four functions:
#
# nmoves(p): number of possible moves for parser p
# nmoves(p,s): number of moves it takes for parser p to parse s (a sentence or a corpus)
# anyvalidmoves(p): whether there are any valid moves in the current state
# move!(p,m): execute move m with parser p
# movecosts(p,h,d,c): the number of gold arcs that become unreachable for each move
#
# Default definitions are provided below.  A new parser Parser{T} will
# fall back to these if it doesn't override them.

# By default our initial state is [w1][w2,...,wn] and our final state is
# [wroot][].  This is achieved by always performing a SHIFT during
# initialization and ensures 2n-2 moves for each sentence.

init!(p::Parser)=(p.nmove=(p.ndeps<<1);shift(p);p)
anyvalidmoves(p::Parser)=(shiftok(p)||reduceok(p)||leftok(p)||rightok(p))
nmoves(p::Parser)=p.nmove
nmoves(p::Parser,s::Sentence)=((wcnt(s)-1)<<1)
nmoves(p::Parser,c::Corpus)=(n=0; for s in c; n += nmoves(p,s); end; n)

function move!(p::Parser, m::Move)
    if !(1 <= m <= p.nmove); error("Move $m is not supported"); end
    if (m == shiftmove(p));      if shiftok(p);  shift(p);            else; error("Bad move"); end
    elseif in(m, rightmoves(p)); if rightok(p);  right(p,label(p,m)); else; error("Bad move"); end
    elseif in(m, leftmoves(p));  if leftok(p);   left(p,label(p,m));  else; error("Bad move"); end
    elseif (m == reducemove(p)); if reduceok(p); reduce(p);           else; error("Bad move"); end
    else error("Move $m is not supported")
    end
end

function movecosts(p::Parser, head::Pvec, deprel::Dvec,
                   cost::Pvec=Array(Cost,p.nmove))
    (length(head) == p.nword)   ||error("Bad head")
    (length(deprel) == p.nword) ||error("Bad deprel")
    (length(cost) == p.nmove)   ||error("Bad cost")
    fill!(cost, typemax(Cost))
    n0 = p.wptr
    s0 = (p.sptr > 0 ? p.stack[p.sptr] : 0)
    s1 = (p.sptr > 1 ? p.stack[p.sptr-1] : 0)
    n0l=0; for i=1:p.sptr; si=p.stack[i]; head[si]==n0 && p.head[si]==0 && (n0l+=1); end
    s0r=0; for i=p.wptr:p.nword; (head[i]==s0) && (s0r += 1); end

    if shiftok(p)
        cost[shiftmove(p)]  =  shiftcost(p, head, n0l, s0r)
    end
    if reduceok(p)
        cost[reducemove(p)] = reducecost(p, head, n0l, s0r)
    end
    if leftok(p)
        (lh,ld) = leftmovepair(p)                               # left adds the arc (lh,ld)
        lcost = leftcost(p, head, n0l, s0r)
        if (head[ld] == lh)                                     # if this is the correct head
            cost[leftmoves(p)] = lcost + 1                      # +1 for the wrong labels
            cost[leftmove(p,deprel[ld])] -= 1                   # except for the correct label
        else
            cost[leftmoves(p)] = lcost				# otherwise we incur leftcost
        end
    end
    if rightok(p)
        (rh,rd) = rightmovepair(p)                              # right adds the arc (rh,rd)
        rcost = rightcost(p, head, n0l, s0r)
        if (head[rd] == rh)                                     # if this is the correct head
            cost[rightmoves(p)] = rcost+1  			# +1 for the wrong labels
            cost[rightmove(p,deprel[rd])] -= 1                  # except for the correct label
        else                                                    #
            cost[rightmoves(p)] = rcost 			# otherwise we incur rightcost
        end
    end
    return cost
end # movecosts

function truecost(p::Parser, s::Sentence)
    cost = 0
    @inbounds for i=1:p.nword
        if (p.head[i] != s.head[i]) || (p.deprel[i] != s.deprel[i])
            cost += 1
        end
    end
    return cost
end

################################ ArcEager,GN13 ##################
# We provide default fallback definitions for Parser based on the
# ArcEager system from [N03,GN13].

typealias ArcEager Parser{:ArcEager}
typealias ArcEager13 Parser{:ArcEager,:GN13}

# In the arc-eager system (N03), a configuration c= (σ,β,A) consists of
# a stack σ, a buffer β, and a set A of dependency arcs.

# There are four types of moves:
# SHIFT[(σ, b|β, A)] = (σ|b, β, A)
# REDUCE[(σ|s, β, A)] = (σ, β, A)
# RIGHT_lb[(σ|s, b|β, A)] = (σ|s|b, β, A∪{(s,lb,b)})
# LEFT_lb[(σ|s, b|β, A)] = (σ, b|β, A∪{(b,lb,s)})

leftmovepair(p::Parser)=(p.wptr,p.stack[p.sptr])
rightmovepair(p::Parser)=(p.stack[p.sptr],p.wptr)

shift(p::Parser)=(p.stack[p.sptr+=1]=p.wptr; p.wptr+=1)
reduce(p::Parser)=(p.sptr-=1)
left(p::Parser, l::DepRel)=(arc!(p, p.wptr, p.stack[p.sptr], l); reduce(p))
right(p::Parser, l::DepRel)=(arc!(p, p.stack[p.sptr], p.wptr, l); shift(p))

# Moves are represented by integers 1..p.nmove
# Default order is SHIFT,L2,R2,..,L[ndeps],R[ndeps],REDUCE

shiftmove(p::Parser)=1
reducemove(p::Parser)=p.nmove
rightmoves(p::Parser)=(3:2:(p.nmove-1))
leftmoves(p::Parser)=(2:2:(p.nmove-2))

# Dependency labels (deprel) are represented by integers 1..p.ndeps.
# The special ROOT deprel is represented by 1 and has no associated move.

rightmove(p::Parser,l::DepRel)=(l<<1-1)
leftmove(p::Parser,l::DepRel)=(l<<1-2)
label(p::Parser,m::Move)=convert(DepRel,m>>1+1)

# GN13 has the following preconditions for moves:
#
# "There is a precondition on the RIGHT and SHIFT transitions to be
# legal only when b != ROOT (p.wptr <= p.nword), and for LEFT, RIGHT and
# REDUCE to be legal only when the stack is non-empty (p.sptr >
# 0). Moreover, LEFT is only legal when s does not have a parent in A,
# and REDUCE when s does have a parent in A."
#
# We implement the ROOT-LEFT as a REDUCE, our LEFT gets an extra
# condition (p.wptr <= p.nword) and REDUCE gets an extra option.

shiftok(p::Parser)=(p.wptr <= p.nword)
reduceok(p::Parser)=((p.sptr > 0) && (s0head(p) || s0head0(p)))
rightok(p::Parser)=((p.wptr <= p.nword) && (p.sptr > 0))
leftok(p::Parser)=((p.wptr <= p.nword) && (p.sptr > 0) && !s0head(p))

s0head(p::Parser)=(p.head[p.stack[p.sptr]] != 0)
s0head0(p::Parser)=((p.wptr > p.nword) && (p.sptr > 1))

# movecosts() counts gold arcs that become impossible after each move.
# A token starts life without any arcs in the buffer.  It moves to the
# head of the buffer (n0) with shift or right moves (each right ends
# with a shift).  Left deps are acquired first while at n0 using left
# moves (each left move ends with a reduce) possibly interspersed with
# other reduces (to get rid of s0 that already have heads).  The token
# moves to s0 with a right or shift.  Then rdeps are acquired while at
# s0 using right moves.  Head is acquired as n0 before rdeps (which
# moves the token to s0, so buffer words never have heads) or as s0
# after rdeps.

function shiftcost(p::Parser, head::AbstractArray, n0l::Int, s0r::Int)
    # eager: n0 gets no more ldeps or lhead
    (n0l + (0 < findprev(p.stack, head[p.wptr], p.sptr)))
end

function reducecost(p::Parser, head::AbstractArray, n0l::Int, s0r::Int)
    s0r # eager: s0 gets no more rdeps
end

function leftcost(p::Parser, head::AbstractArray, n0l::Int, s0r::Int)
    # eager: s0 gets no more rdeps, rhead>n0, 0head
    s0 = p.stack[p.sptr]; s0h = head[s0]
    (s0r + (s0h > p.wptr) + (s0h == 0))
end

function rightcost(p::Parser, head::AbstractArray, n0l::Int, s0r::Int)
    # eager: n0 gets no more ldeps, rhead, 0head, or lhead<s0
    n0 = p.wptr; n0h = head[n0]
    (n0l + (n0h > n0) + (n0h == 0) + (0 < findprev(p.stack, n0h, p.sptr-1)))
end

################################ ArcEagerR1 ##########################
# ArcEagerR1 is a modification of ArcEager to ensure a single rootword

typealias ArcEagerR1 Parser{:ArcEager,:R1}

# In GN13 the initial configuration has an empty stack, and a buffer
# with special symbol ROOT to the right of all the words at w[n+1], i.e.
# [][w1,w2,...,wn,ROOT]. The final state in GN13 is only the ROOT token
# in the buffer and nothing on the stack, i.e. [][ROOT].  The special
# ROOT token can have more than one child, thus the sentence does not
# always have a unique rootword.

# GN13 ends up with multiple rootwords by shifting them to the stack and
# popping them with left moves when the only token left in the buffer is
# ROOT.  We can prevent multiple rootwords by:
# (1) avoid SHIFT on last buffer word unless the stack is empty.
# (2) avoid RIGHT on last buffer word if multiple headless words in stack.

shiftok(p::ArcEagerR1)=((p.wptr < p.nword) || ((p.wptr == p.nword) && (p.sptr == 0)))
rightok(p::ArcEagerR1)=((p.sptr > 0) && ((p.wptr < p.nword) || ((p.wptr == p.nword) && (headless(p)==1))))
leftok(p::ArcEagerR1)=((p.sptr > 0) && (p.wptr <= p.nword) && (p.head[p.stack[p.sptr]] == 0))
reduceok(p::ArcEagerR1)=((p.sptr > 0) && (p.head[p.stack[p.sptr]] != 0))
headless(p::ArcEagerR1)=(h=0;for i=1:p.sptr; si=p.stack[i]; (p.head[si]==0) && (h+=1); end; h)

# Proof: Buffer words never have heads, some stack words do.  As long as
# we have words in the buffer we can recover a single root.  Once the
# buffer is empty, nothing can go back into it.  With an empty buffer,
# REDUCE is the only possible move, (we don't allow left moves with an
# empty buffer) now new non-root arcs can be added.  So before we get to
# the empty buffer we need to make sure the stack will contain a single
# headless rootword.  Stack words are either headless or are next to
# their parents, i.e. a non-empty stack always contains at least one
# headless word.  SHIFTing the last word into a non-empty stack would
# add another, so can't do it.  RIGHT with the last word would not
# create another headless stack word, but we have to make sure there is
# only one before doing RIGHT on the last word.

# WolframAlpha suggested the following succinct form for anyvalidmoves(p).
# REDUCE if ((p.sptr > 0) && (p.head[p.stack[p.sptr]] != 0))
# otherwise (p.wptr <= p.nword) is true
# SHIFT if (p.sptr == 0)
# otherwise (p.sptr > 0) is true
# so (p.head[p.stack[p.sptr]] != 0) is false
# and LEFT is possible

anyvalidmoves(p::ArcEagerR1)=((p.wptr <= p.nword) || ((p.sptr > 0) && (p.head[p.stack[p.sptr]] != 0)))

################################ ArcHybrid ####################
# ArcHybrid13 is the ArcHybrid system described in [KGS11,GN13]

typealias ArcHybrid Parser{:ArcHybrid}
typealias ArcHybrid13 Parser{:ArcHybrid,:GN13}

# In the arc-hybrid system (KGS11), a configuration c= (σ,β,A) consists
# of a stack σ, a buffer β, and a set A of dependency arcs.

# There are three types of moves:
# SHIFT[(σ, b|β, A)] = (σ|b, β, A)
# RIGHT_lb[(σ|s1|s0, β, A)] = (σ|s1, β, A ∪ {(s1, lb, s0)})
# LEFT_lb[(σ|s, b|β, A)] = (σ, b|β, A ∪ {(b, lb, s)})

# The hybrid system has a different RIGHT action linking (s1,s0).

right(p::ArcHybrid,l::DepRel)=(arc!(p, p.stack[p.sptr-1], p.stack[p.sptr], l); reduce(p))
rightmovepair(p::ArcHybrid)=(p.stack[p.sptr-1], p.stack[p.sptr])

# GN13 has the following preconditions for ArcHybrid moves:
#
# "There is a precondition on RIGHT to be legal only when the stack has
# at least two elements, and on LEFT to be legal only when the stack is
# non-empty and s != ROOT."

shiftok(p::ArcHybrid)=(p.wptr <= p.nword)
rightok(p::ArcHybrid)=(p.sptr > 1)
leftok(p::ArcHybrid)=((p.wptr <= p.nword) && (p.sptr > 0))

# We introduce an additional REDUCE move when the stack has a single
# word to represent the ROOT linkage.  We terminate with the last word
# in stack, so REDUCE has a precondition of a non-empty buffer.

reduceok(p::ArcHybrid)=((p.sptr == 1) && (p.wptr <= p.nword))
anyvalidmoves(p::ArcHybrid)=((p.wptr <= p.nword) || (p.sptr > 1))

# movecosts() counts gold arcs that become impossible after possible
# moves.  Tokens start their lifecycle in the buffer without links.
# They move to the top of the buffer (n0) with SHIFT moves.  There they
# acquire left dependents using LEFT moves.  After that a single SHIFT
# moves them to the top of the stack (s0).  There they acquire right
# dependents with SHIFT-RIGHT pairs.  Finally from s0 they acquire a
# head with a LEFT or RIGHT move.  Any token from the buffer may become
# the right head but only s1 from the stack may become a left head for
# s0.  The parser terminates with a single word at s0 whose head is ROOT
# (represented as head=deprel=0).
#
# 1. SHIFT moves n0 to s0: n0 cannot acquire left dependents after a
# shift.  Also it can no longer get a head from the stack to the left
# of s0 or get a root head if there is s0: (0+s\s0,n0) + (n0,s)
#
# 2. R0MOVE pops s0 linking it to ROOT: s0 cannot acquire a head or
# dependent from the buffer after right: (s0,b) + (b,s0)
#
# 3. RIGHT adds (s1,s0): s0 cannot acquire a head or dependent from
# the buffer after right: (s0,b) + (b,s0)
#
# 4. LEFT adds (n0,s0): s0 cannot acquire s1 or 0 (if there is no s1)
# or ni (i>0) as head.  It also cannot acquire any more right
# children: (s0,b) + (b\n0,s0) + (s1 or 0,s0)

function shiftcost(p::ArcHybrid, head::AbstractArray, n0l::Int, s0r::Int)
    # n0 gets no more ldeps or lhead<s0 or root head if there is s0
    n0 = p.wptr; n0h = head[n0]
    (n0l + (findprev(p.stack, n0h, p.sptr-1) > 0) + ((n0h==0) && (p.sptr>0)))
end

function reducecost(p::ArcHybrid, head::AbstractArray, n0l::Int, s0r::Int)
    # s0 gets no more rdeps or rhead
    s0 = p.stack[p.sptr]; s0h = head[s0]
    (s0r + (s0h >= p.wptr))
end

function leftcost(p::ArcHybrid, head::AbstractArray, n0l::Int, s0r::Int)
    # s0 gets no more rdeps, rhead>n0, s1 head or 0head if alone
    s0 = p.stack[p.sptr]; s0h = head[s0]
    (s0r + ((s0h > p.wptr) || ((p.sptr == 1) && (s0h == 0)) || ((p.sptr > 1) && (s0h == p.stack[p.sptr-1]))))
end

function rightcost(p::ArcHybrid, head::AbstractArray, n0l::Int, s0r::Int)
    # s0 gets no more rdeps or rhead
    s0 = p.stack[p.sptr]; s0h = head[s0]
    (s0r + (s0h >= p.wptr))
end

################################ ArcHybridR1 #############################
# ArcHybridR1 is a modification of ArcHybrid to ensure a single rootword

typealias ArcHybridR1 Parser{:ArcHybrid,:R1}

# It turns out we can ensure a single rootword by simply not allowing
# the reduce move (i.e. the ROOT-RIGHT move)

reduceok(p::ArcHybridR1)=false
reducemove(p::ArcHybridR1)=nothing

# This gives us one less legal move.

init!(p::ArcHybridR1)=(p.nmove=(p.ndeps<<1-1);shift(p);p)

# Moves are represented by integers 1..p.nmove, 0 is not valid.
# They correspond to SHIFT,L2,R2,...,L[ndeps],R[ndeps]

leftmoves(p::ArcHybridR1)=(2:2:(p.nmove-1))
rightmoves(p::ArcHybridR1)=(3:2:p.nmove)

################################################################
# Some general utility functions:

function arc!(p::Parser, h::Position, d::Position, l::DepRel)
    p.head[d] = h
    p.deprel[d] = l
end

function Base.:(==)(a::Parser, b::Parser)
    if typeof(a)!=typeof(b); return false; end
    for f in fieldnames(a)
        if getfield(a,f) != getfield(b,f); return false; end
    end
    return true
end

function Base.copy!(dst::Parser, src::Parser)
    (dst.nword == src.nword && dst.ndeps == src.ndeps && dst.nmove == src.nmove) || error("Incompatible parsers")
    dst.wptr = src.wptr
    dst.sptr = src.sptr
    copy!(dst.stack, src.stack)
    copy!(dst.head, src.head)
    copy!(dst.deprel, src.deprel)
    dst
end # copy!

Base.copy{T<:Parser}(src::T)=copy!(T(src.nword, src.ndeps), src)

function reset!(p::Parser)
    p.wptr = 1
    p.sptr = 0
    p.stack[:] = 0
    p.head[:] = 0
    p.deprel[:] = 0
    init!(p)
end

reset!{T<:Parser}(pa::Vector{T})=(for p in pa; reset!(p); end)

Base.show(io::IO, p::Parser)=print(io, map(Int,p.stack[1:p.sptr]), Int[p.wptr], map(Int,p.head))

function Base.hash(p::Parser, h::UInt)
    h += UInt(0xac816432c42f6df9)
    for n in fieldnames(p)
        h = hash(p.(n), h)
    end
    return h
end


text:
using Kinetic, Plots, Flux
using KitBase.JLD2
using KitBase.ProgressMeter: @showprogress

cd(@__DIR__)

set = Setup(case = "sod", space = "1d2f1v", maxTime = 0.15)
ps = PSpace1D(0.0, 1.0, 200, 1)
vs = VSpace1D(-5.0, 5.0, 100)
gas = Gas(Kn = 1e-4, K = 2.0, γ = 5/3)
ib = IB1F(ib_sod(set, ps, vs, gas)...)
ks = SolverSet(set, ps, vs, gas, ib)
ctr, face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ctr[1].w)
avg = zero(res)
regime = ones(Int, axes(ks.ps.x))
regime0 = deepcopy(regime)

@time for iter = 1:nt
    #reconstruct!(ks, ctr)
    evolve!(ks, ctr, face, dt)
    update!(ks, ctr, face, dt, res)

    t += dt
end

plot(ks, ctr)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using MD5
using Base.Test

const truth = Dict{String, Vector{UInt8}}()
truth[""] = UInt8[0xd4, 0x1d, 0x8c, 0xd9, 0x8f, 0x00, 0xb2, 0x04, 0xe9, 0x80, 0x09, 0x98, 0xec, 0xf8, 0x42, 0x7e]
truth["3asdf"] = UInt8[0x26, 0xcf, 0xe6, 0x62, 0xea, 0x61, 0xa3, 0x9d, 0x25, 0xf1, 0x61, 0x3f, 0xda, 0x34, 0xa0, 0x5a]
truth["a"] = UInt8[0x0c, 0xc1, 0x75, 0xb9, 0xc0, 0xf1, 0xb6, 0xa8, 0x31, 0xc3, 0x99, 0xe2, 0x69, 0x77, 0x26, 0x61]

@testset "compare with truth" begin
    for (s, h) ∈ truth
        @test md5(s) == h
    end
end

include("nettle.jl")


text:
# Run advection tests

push!(LOAD_PATH, abspath(joinpath(pwd(), "..")))

using MPI
using PDESolver
#using Base.Test
using Base.Test
using ODLCommonTools
using Utils
using PdePumiInterface  # common mesh interface - pumi
using SummationByParts  # SBP operators
using AdvectionEquationMod
using NonlinearSolvers   # non-linear solversa
using OptimizationInterface
using ArrayViews
import ODLCommonTools.sview
import ArrayViews.view
using Input
using LinearSolvers
using PETSc2
using Jacobian

function clean_dict(collection)
  for i in keys(collection)
    delete!(collection, i)
  end
end

#------------------------------------------------------------------------------
# define tests and tags
using TestSystem
# define tags that will be used
include("../tags.jl")

# test list
global const AdvectionTests = TestList()

include("test_frontend.jl")
include("test_linearsolver.jl")
include("test_lowlevel.jl")
include("test_eqn_deepcopy.jl")
include("test_3d.jl")
include("test_utils.jl")
include("test_gamma.jl")
include("test_mms.jl")
include("test_jac.jl")
include("test_GLS2.jl")
# test_newton.jl?
include("test_dg.jl")
include("test_staggered.jl")
include("test_adjoint.jl")
include("test_parallel_serialpart.jl")
include( "./energy/runtests.jl")
#cd("./Nonlinearsolvers/")
include(joinpath(pwd(), "Nonlinearsolvers", "runtests_serial.jl"))
#cd("../")

#------------------------------------------------------------------------------
# run tests
@testset "----- Running Advection tests -----" begin
  runTestSystem(AdvectionTests, solvePDE, ARGS)
end

println("finished running tests")

#------------------------------------------------------------------------------
# cleanup


text:
using YaoBase

"""
    content(x)

Returns the content of `x`.
"""
content(x::AbstractContainer) = x.content

"""
    chcontent(x, blk)

Create a similar block of `x` and change its content to blk.
"""
chcontent(x::AbstractContainer, blk) = chsubblocks(x, blk)

subblocks(x::AbstractContainer) = (content(x),)
# NOTE: there's only one block inside, so we expand the iterator
# this would error if there's more than one block in it. But will
# work if there exactly one block.
chsubblocks(x::AbstractContainer, it) = chsubblocks(x, it...)

# throw better error msg when no chsubblocks is overloaded
# for this container block
# since every AbstractContainer should overload this method
chsubblocks(x::AbstractContainer, it::AbstractBlock) =
    throw(NotImplementedError(:chsubblocks, (x, it)))

# TODO:
#   - use simple traits instead
#   - each property should have a trait
# NOTE: this is a holy trait, no overhead, don't use methods on this
"""
    PropertyTrait
    PropertyTrait(::AbstractContainer) -> PropertyTrait

Typically, it can be `PreserveAll()` for those containers that do not change `ishermitian`, `isunitary` and `isreflexive` properties, otherwise fallback to `PreserveNothing`.
"""
abstract type PropertyTrait end
struct PreserveAll <: PropertyTrait end
struct PreserveNothing <: PropertyTrait end

PropertyTrait(c::AbstractContainer) = PreserveNothing()

for METHOD in (:ishermitian, :isreflexive, :isunitary)
    @eval begin
        # forward to trait
        YaoBase.$METHOD(x::AbstractContainer) = $METHOD(PropertyTrait(x), x)
        # forward parent block property
        YaoBase.$METHOD(::PreserveAll, c::AbstractContainer) = $METHOD(content(c))
        # forward to default property by calculating the matrix
        YaoBase.$METHOD(::PreserveNothing, c::AbstractContainer) = $METHOD(mat(c))
    end
end

function Base.:(==)(lhs::AbstractContainer{BT,N}, rhs::AbstractContainer{BT,N}) where {BT,N}
    return content(lhs) == content(rhs)
end

include("chain.jl")
include("kron.jl")
include("control.jl")
include("put_block.jl")
include("repeated.jl")
include("subroutine.jl")
include("reduce.jl")
include("unitary_channel.jl")

chsubblocks(x::ChainBlock, it::AbstractBlock) = chsubblocks(x, (it,))
chsubblocks(x::KronBlock, it::AbstractBlock) = chsubblocks(x, (it,))
chsubblocks(x::Add, it::AbstractBlock) = chsubblocks(x, (it,))

# tag blocks
include("tag/tag.jl")
include("tag/cache.jl")
include("tag/dagger.jl")
include("tag/scale.jl")


text:
#%%
using CSV
using DataFrames
#%%
# df = DataFrame()
#%%

statsncb = CSV.File("ncb-1000.csv") |> DataFrame
statslst = CSV.File("lst-1000.csv") |> DataFrame
statslbf = CSV.File("lbfgs-1000.csv") |> DataFrame

#%%
#seleciona colunas
statsncb_uteis = statsncb[:, [:name, :nvar, :status, :elapsed_time, :iter]]
statslst_uteis = statslst[:, [:name, :nvar, :status, :elapsed_time, :iter]]
statslbf_uteis = statslbf[:, [:name, :nvar, :status, :elapsed_time, :iter]]

#%%
# filtra por status
res1 = statsncb_uteis[statsncb_uteis[:status] .== "first_order", :]
res2 = statslst_uteis[statslst_uteis[:status] .== "first_order", :]
res3 = statslbf_uteis[statslbf_uteis[:status] .== "first_order", :]
#%%

#%%
#ordena
naores1 = statsncb_uteis[statsncb_uteis[:status] .!= "first_order", :]
sort!(naores1, [:nvar, :name], rev = [true, false])
println(naores1)
#%%
naores2 = statslst_uteis[statslst_uteis[:status] .!= "first_order", :]
sort!(naores2, [:nvar, :name], rev = [true, false])
println(naores2)
#%%
naores3 = statslbf_uteis[statslbf_uteis[:status] .!= "first_order", :]
sort!(naores3, [:nvar, :name], rev = [true, false])
println(naores3)
#%%
# CSV.write("naoncb-5.csv", naores1)
# CSV.write("naolst-5.csv", naores2)
# CSV.write("naolbfgs-5.csv", naores3)
#%%

#%%
#contagens
conta_1 = by(statsncb_uteis, :status, c -> DataFrame(count = nrow(c)))
print(conta_1)
#%%
conta_2 = by(statslst_uteis, :status, c -> DataFrame(count = nrow(c)))
print(conta_2)
#%%
conta_3 = by(statslbf_uteis, :status, c -> DataFrame(count = nrow(c)))
print(conta_3)
#%%
using Latexify
print(latexify(naores2, env = :table, latex = false))
#%%
A1 = naores1[:,:name]
A2 = naores2[:,:name]
A3 = naores3[:,:name]

m=1
n=1
for i in A3
    verifica = :false
    for j in A2
        if i == j
            verifica = :true
            for k in A1
                if i == k
                    println("$i não resolvido por nenhum dos 3")
                end
            end
        end
    end
    if verifica == :false
        for j in A1
            if i == j
                println("Não por A1, nem por A3: $m:- $i")
                m +=1
                verifica = :true
            end
        end
    end
    if verifica == :false
        #não resolveu A1
        println("##  APENAS: Não por A3:$n- $i")
        n+=1
    end
end
#%%
m=1
n=1
for i in A1, j in A2
    verifica = :false
    for k in A3
        if i == k || j == k
            println("$i não resolvido por nenhum dos 3")
        end
    end
end


text:
function _scatter(V; kwargs...)
	size(V,2)==2 && return scatter(;x=V[:,1], y=V[:,2], kwargs...)
	scatter3d(;x=V[:,1], y=V[:,2], z=V[:,3], kwargs...)
end

function plotsimplices(V, sg, colorBy, colorDict;
	                   drawPoints=true, drawLines=true, drawTriangles=true,
	                   title="",
	                   opacity=0.3, markerSize=5, lineWidth=2,
	                   shapeBy=nothing, shapeDict=nothing,
	                   xLabel="x", yLabel="y", zLabel="z",
	                   legendTitle="",
	                   cameraAttr=attr())
	@assert size(V,2) in 2:3

	traces = GenericTrace[]

	# colorscale setup
	colorScale = nothing
	tOffs, tScale = 0.0, 1.0
	if colorDict==nothing
		# this is the "RdBu" colorscale used by PlotlyJS
		colorScale = ColorScale([0, 0.35, 0.5, 0.6, 0.7, 1], [RGB(5/255,10/255,172/255), RGB(106/255,137/255,247/255), RGB(190/255,190/255,190/255), RGB(220/255,170/255,132/255), RGB(230/255,145/255,90/255), RGB(178/255,10/255,28/255)])

		if !isempty(skipmissing(colorBy))
			m,M = Float64.(extrema(skipmissing(colorBy)))
			tOffs,tScale = M>m ? (m, 1.0./(M-m)) : (m-0.5, 1.0) # if there is only one unique value, map it to the middle of the colorscale
		end
	end

	if drawTriangles && size(V,2)>2 # Triangles not supported in 2d plot for now
		TRIANGLE_LIMIT = 2_000_000

		triangleInds = Int[]
		for c=1:size(sg.G,2)
			ind = findall(sg.G[:,c])
			isempty(ind) && continue

			length(ind)<3 && continue # no triangles

			# slow and ugly solution
			for tri in IterTools.subsets(ind,3)
				append!(triangleInds, sort(tri))
			end
		end
		triangleInds = reshape(triangleInds,3,:)
		triangleInds = unique(triangleInds,dims=2) # remove duplicates
		triangleInds .-= 1 # PlotlyJS wants zero-based indices

		if size(triangleInds,2)>TRIANGLE_LIMIT
			@warn "More than $TRIANGLE_LIMIT triangles to plot, disabling triangle plotting for performance reasons."
		else
			if colorDict==nothing
				vertexColor = [ ismissing(t) ? RGB(0.,0.,0.) : lookup(colorScale, (t-tOffs)*tScale) for t in colorBy ]
			else
				vertexColor = getindex.((colorDict,), colorBy)
			end
			push!(traces, mesh3d(; x=V[:,1],y=V[:,2],z=V[:,3],i=triangleInds[1,:],j=triangleInds[2,:],k=triangleInds[3,:], vertexcolor=vertexColor, opacity=opacity, lighting=attr(ambient=1.0,diffuse=0.0,specular=0.0,fresnel=0.0), showlegend=false))

		end
	end

	if drawLines
		LINE_LIMIT = 300_000

		GK = sg.G*sg.G'
		nbrLines = div(nnz(GK) - sum(i->GK[i,i]>0, 1:size(GK,1)), 2) # number of elements above the diagonal

		if nbrLines > LINE_LIMIT
			@warn "More than $LINE_LIMIT lines to plot, disabling line plotting for performance reasons."
		else
			# every third point is "nothing" to make lines disconnected
			VL = repeat(Union{Float64,Nothing}[nothing], nbrLines*3, size(V,2))
			colorsRGB = repeat([RGB(0.,0.,0.)], nbrLines*3)

			lineInd = 1
			rows = rowvals(GK)
			for c = 1:size(GK,2)
				for k in nzrange(GK, c)
					r = rows[k]
					r>=c && break # just use upper triangular part
					VL[lineInd,   :] .= V[r,:]
					VL[lineInd+1, :] .= V[c,:]

					if colorDict==nothing
						t1,t2 = colorBy[r],colorBy[c]
						col1 = ismissing(t1) ? RGB(0.,0.,0.) : lookup(colorScale, (t1-tOffs)*tScale)
						col2 = ismissing(t2) ? RGB(0.,0.,0.) : lookup(colorScale, (t2-tOffs)*tScale)
					else
						col1 = colorDict[colorBy[r]]
						col2 = colorDict[colorBy[c]]
					end

					colorsRGB[lineInd]   = col1
					colorsRGB[lineInd+1] = col2

					lineInd+=3
				end
			end
			@assert lineInd == nbrLines*3+1

			if size(V,2)==2
				# PLOTLY DOESN'T HANDLE PER VERTEX COLORS FOR 2D LINES. THIS IS A STUPID WORKAROUND TO MAKE LINES BLACK IN 2D.
				push!(traces, _scatter(VL; mode="lines", line=attr(color=colorant"black", width=lineWidth), showlegend=false))
			else
				push!(traces, _scatter(VL; mode="lines", line=attr(color=colorsRGB, width=lineWidth), showlegend=false))
			end
		end
	end

	# plot each group in different colors
	if drawPoints
		# TODO: merge the two cases below?

		if colorDict != nothing
			for cb in unique(colorBy)
				ind = findall(colorBy.==cb )
				col = colorDict[cb]

				extras = []
				length(colorDict)==1 && push!(extras, (;showlegend=false))
				shapeBy!=nothing && shapeDict!=nothing && push!(extras, (marker_symbol=[shapeDict[k] for k in shapeBy[ind]],))
				isempty(extras) || (extras = pairs(extras...))

				points = _scatter(V[ind,:]; mode="markers", marker_color=col, marker_size=markerSize, marker_line_width=0, name=string(cb), extras...)
				push!(traces, points)
			end
		else
			extras = []
			shapeBy!=nothing && shapeDict!=nothing && push!(extras, (marker_symbol=[shapeDict[k] for k in shapeBy],))
			isempty(extras) || (extras = pairs(extras...))

			V2 = V
			c = colorBy

			# handle missing values by plotting them in another trace (with a different color)
			if any(ismissing, c)
				mask = ismissing.(c)
				pointsNA = _scatter(V2[mask,:]; mode="markers", marker=attr(color=colorant"black", size=markerSize, line_width=0), name="", showlegend=false, extras...)
				push!(traces, pointsNA)

				V2 = V2[.!mask,:]
				c = disallowmissing(c[.!mask])
			end

			points = _scatter(V2; mode="markers", marker=attr(color=c, colorscale=to_list(colorScale), showscale=true, size=markerSize, line_width=0, colorbar=attr(title=legendTitle)), name="", showlegend=false, extras...)
			push!(traces, points)
		end
	end

	if size(V,2)==2
		layout = Layout(title=title,
		                xaxis=attr(title=xLabel, constrain="domain"),
		                yaxis=attr(title=yLabel, constrain="domain", scaleanchor='x'),
		                legend=attr(title_text=legendTitle, itemsizing="constant"),
		                hovermode="closest")
	else
		layout = Layout(title=title,
		                scene=attr(xaxis=attr(title=xLabel), yaxis=attr(title=yLabel), zaxis=attr(title=zLabel), camera=cameraAttr, aspectmode="data"),
		                legend=attr(title_text=legendTitle, itemsizing="constant"))
	end
	traces, layout # return plot args rather than plot because of threading issues.
end

_distinguishable_colors(n) = distinguishable_colors(n+1,colorant"white")[2:end]

function colordict(x::AbstractArray)
	k = unique(x)
	Dict(s=>c for (s,c) in zip(k,_distinguishable_colors(length(k))))
end

struct ColorScale
	t::Vector{Float64} # increasing values between 0 and 1. First must be 0 and last must be 1.
	colors::Vector{RGB{Float64}} # must be same length as t
end
to_list(colorScale::ColorScale) = [ [t,c] for (t,c) in zip(colorScale.t, colorScale.colors)]

function lookup(colorScale::ColorScale, t::Float64)
	i2 = searchsortedfirst(colorScale.t, t)
	i2>length(colorScale.t) && return colorScale.colors[end]
	i2==1 && return colorScale.colors[1]
	i1 = i2-1
	α = (t-colorScale.t[i1]) / (colorScale.t[i2]-colorScale.t[i1])
	weighted_color_mean(α, colorScale.colors[i2], colorScale.colors[i1]) # NB: order because α=1 means first color in weighted_color_mean
end

const SHAPES = ["circle","square","diamond","cross","x","triangle-up","triangle-down","triangle-left","triangle-right","triangle-ne","triangle-se","triangle-sw","triangle-nw","pentagon","hexagon","hexagon2","octagon","star","hexagram","star-triangle-up","star-triangle-down","star-square","star-diamond","diamond-tall","diamond-wide","hourglass","bowtie","circle-cross","circle-x","square-cross","square-x","diamond-cross","diamond-x","cross-thin","x-thin","asterisk","hash","y-up","y-down","y-left","y-right","line-ew","line-ns","line-ne","line-nw","circle-open","square-open","diamond-open","cross-open","x-open","triangle-up-open","triangle-down-open","triangle-left-open","triangle-right-open","triangle-ne-open","triangle-se-open","triangle-sw-open","triangle-nw-open","pentagon-open","hexagon-open","hexagon2-open","octagon-open","star-open","hexagram-open","star-triangle-up-open","star-triangle-down-open","star-square-open","star-diamond-open","diamond-tall-open","diamond-wide-open","hourglass-open","bowtie-open","circle-cross-open","circle-x-open","square-cross-open","square-x-open","diamond-cross-open","diamond-x-open","cross-thin-open","x-thin-open","asterisk-open","hash-open","y-up-open","y-down-open","y-left-open","y-right-open","line-ew-open","line-ns-open","line-ne-open","line-nw-open","","circle-dot","square-dot","diamond-dot","cross-dot","x-dot","triangle-up-dot","triangle-down-dot","triangle-left-dot","triangle-right-dot","triangle-ne-dot","triangle-se-dot","triangle-sw-dot","triangle-nw-dot","pentagon-dot","hexagon-dot","hexagon2-dot","octagon-dot","star-dot","hexagram-dot","star-triangle-up-dot","star-triangle-down-dot","star-square-dot","star-diamond-dot","diamond-tall-dot","diamond-wide-dot","hash-dot","circle-open-dot","square-open-dot","diamond-open-dot","cross-open-dot","x-open-dot","triangle-up-open-dot","triangle-down-open-dot","triangle-left-open-dot","triangle-right-open-dot","triangle-ne-open-dot","triangle-se-open-dot","triangle-sw-open-dot","triangle-nw-open-dot","pentagon-open-dot","hexagon-open-dot","hexagon2-open-dot","octagon-open-dot","star-open-dot","hexagram-open-dot","star-triangle-up-open-dot","star-triangle-down-open-dot","star-square-open-dot","star-diamond-open-dot","diamond-tall-open-dot","diamond-wide-open-dot","hash-open-dot"]
shapedict(x::AbstractArray) = Dict(v=>SHAPES[mod(i-1,length(SHAPES))+1] for (i,v) in enumerate(unique(x)))


text:
"""
    fr2ped(dir, list, ped, allele)
---

# Description
Given a `list` of final report files, this function merge them into one file
`ped`.

The function will analyse the first file in the list to see if `allele` format
e.g., `Top`, `AB`", `Forward`, `Design`, or `Plus`, is available.  The column
number will be recorded. If not found, it raises an error.  By default, `allele`
is of `Top` format.

# Example
```
dir = "path-to-final-reports"
list = readdir(dir)
fr2ped(dir, list, "path-to/plink.ped", allele="Top")
```
"""
function fr2ped(dir::AbstractString,
                list,
                ped::AbstractString,
                allele::AbstractString = "Top")
    item("Analyse fields")
    nlc = x = y = 0
    open(joinpath(dir, list[1]), "r") do io
        # skip header
        for line in eachline(io)
            occursin("[Data]", line) && break
            occursin("Num SNPs", line) && ( nlc = parse(Int, split(line, '\t')[2]) )
        end

        # find allele columns
        header = split(readline(io), '\t')
        fields = split(readline(io), '\t') # first record
        n = 0
        for f in fields
            n += 1
            n <= length(header) || break
            (length(f) > 0) && println(n, '\t', header[n], '\t', f)
            if occursin("Allele1", header[n]) && split(header[n], ' ')[3] == allele
                x, y = n, n+1
            end
        end
        x > 0 || error("Allele $allele designation not found in these final reports")
        length(fields[x]) > 0 || error("The columns found for '$allele' are empty")
        n = 1
        for _ in eachline(io)
            n += 1
        end
        n == nlc || error("Number of SNP available, $n is different from $nlc stated in the header")
    end
    message("Genotypes in $allele format are in column $x and $y.\n" *
            "A total of $nlc loci found in this report")

    item("Merge the reports in '$dir' to '$ped'")
    function isvalid(c::Char)
        valid = Set(['A', 'C', 'G', 'T', 'B'])
        return (c in valid)
    end
    cnt = 0
    open(ped, "w") do io
        for report in list
            cnt += 1
            id = ""::String
            gt = fill('0', nlc*2) # default plink missing value is '0'
            ilc = 0
            print("\r Dealing with $report, number $cnt                    ")
            open(joinpath(dir, report), "r") do rp
                for line in eachline(rp) # skip header
                    occursin("SNP Name", line) && break
                end
                for line in eachline(rp)
                    ilc += 1
                    id, a, b = [split(line, '\t')[i] for i in [2, x, y]]
                    isvalid(a[1]) && (gt[ilc*2-1] = a[1])
                    isvalid(b[1]) && (gt[ilc * 2] = b[1])
                end
            end
            write(io, "dummy $id 0 0 0 -9 ", join(gt, ' '), '\n')
        end
    end
    done()
end


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Finance

# **Originally Contributed by**: Arpit Bhatia

# Optimization models play an increasingly important role in financial
# decisions. Many computational finance problems can be solved efficiently using
# modern optimization techniques.
#
# In this tutorial we will discuss two such examples taken from the book
# [Optimization Methods in Finance](https://doi.org/10.1017/9781107297340).

# This tutorial uses the following packages
using JuMP
import GLPK

# ## Short Term Financing

# Corporations routinely face the problem of financing short term cash
# commitments such as the following:

# **Month**        |Jan    |Feb    |Mar    |Apr    |May    |Jun
# :-----:          |:-----:|:-----:|:-----:|:-----:|:-----:|:-----:
# **Net Cash Flow**|-150   |-100   |200    |-200   |50     |300

# Net cash flow requirements are given in thousands of dollars. The company has
# the following sources of funds:

# - A line of credit of up to \$100K at an interest rate of 1% per month,
# - In any one of the first three months, it can issue 90-day commercial paper
#   bearing a total interest of 2% for the 3-month period,
# - Excess funds can be invested at an interest rate of 0.3% per month.

# Our task is to find out the most economical way to use these 3 sources such
# that we end up with the most amount of money at the end of June.
#
# We model this problem in the following manner:

# We will use the following decision variables:
# - the amount $u_{i}$ drawn from the line of credit in month $i$
# - the amount $v_{i}$ of commercial paper issued in month $i$
# - the excess funds $w_{i}$ in month $i$

# Here we have three types of constraints:
# 1. for every month, cash inflow = cash outflow for each month
# 2. upper bounds on $u_{i}$
# 3. nonnegativity of the decision variables $u_{i}$, $v_{i}$ and $w_{i}$.

# Our objective will be to simply maximimse the company's wealth in June, which
# say we represent with the variable $m$.

financing = Model(GLPK.Optimizer)

@variables(financing, begin
    0 <= u[1:5] <= 100
    0 <= v[1:3]
    0 <= w[1:5]
    m
end)

@objective(financing, Max, m)

@constraints(financing, begin
    u[1] + v[1] - w[1] == 150 # January
    u[2] + v[2] - w[2] - 1.01u[1] + 1.003w[1] == 100 # February
    u[3] + v[3] - w[3] - 1.01u[2] + 1.003w[2] == -200 # March
    u[4] - w[4] - 1.02v[1] - 1.01u[3] + 1.003w[3] == 200 # April
    u[5] - w[5] - 1.02v[2] - 1.01u[4] + 1.003w[4] == -50 # May
    -m - 1.02v[3] - 1.01u[5] + 1.003w[5] == -300 # June
end)

optimize!(financing)

objective_value(financing)

# ## Combinatorial Auctions

# In many auctions, the value that a bidder has for a set of items may not be
# the sum of the values that he has for individual items.
#
# Examples are equity trading, electricity markets, pollution right auctions and
# auctions for airport landing slots.
#
# To take this into account, combinatorial auctions allow the bidders to submit
# bids on combinations of items.

# Let $M=\{1,2, \ldots, m\}$ be the set of items that the auctioneer has to
# sell. A bid is a pair $B_{j}=\left(S_{j}, p_{j}\right)$ where
# $S_{j} \subseteq M$ is a nonempty set of items and $p_{j}$ is the price offer
# for this set.
#
# Suppose that the auctioneer has received $n$ bids
# $B_{1}, B_{2}, \ldots, B_{n}.$ The goal of this problem is to help an
# auctioneer determine the winners in order to maximize his revenue.

# We model this problem by taking a decision variable $y_{j}$ for every bid. We
# add a constraint that each item $i$ is sold at most once. This gives us the
# following model:

# ```math
# \begin{aligned}
# \max && \sum_{i=1}^{n} p_{j} y_{j} \\
# \text { s.t. }  && \sum_{j : i \in S_{j}} y_{j} \leq 1 && \forall i=\{1,2 \ldots m\} \\
# && y_{j} \in\{0,1\} && \forall j \in\{1,2 \ldots n\}
# \end{aligned}
# ```

bid_values = [6 3 12 12 8 16]
bid_items = [[1], [2], [3 4], [1 3], [2 4], [1 3 4]]

auction = Model(GLPK.Optimizer)
@variable(auction, y[1:6], Bin)
@objective(auction, Max, sum(y' .* bid_values))
for i in 1:6
    @constraint(auction, sum(y[j] for j in 1:6 if i in bid_items[j]) <= 1)
end

optimize!(auction)

objective_value(auction)

#-

value.(y)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# # Introduction

# ElTopo.jl is a interface to a C++ library [eltopo](https://github.com/tysonbrochu/eltopo). ElTopo the C++ library is a package for tracking dynamic surfaces represented as triangle meshes in 3D. It robustly handles topology changes such as merging and pinching off, while adaptively maintaining a tangle-free, high-quality triangulation. Suitable to simulate time dependant boundary problems, such as droplet behaviour uppon external field, interface dynamics of liquids in a container, etc.

# Currently, the wrapper works only on Linux (Help wanted! See [ElTopoBuilder](https://github.com/akels/ElTopoBuilder) and [Travis](https://travis-ci.org/akels/ElTopo.jl/jobs/574602453 ) ) and does not expose internal triangular mesh data structure methods. Nevertheless, it provides the mesh stabilization method with expossed parameters to control it. And together with `SurfaceTopology.jl` should be useful for soft surface dynamics calculations.

# ##  The API

# The wrapper focuses on the high-level functionality provided in [ElTopo](https://github.com/tysonbrochu/eltopo). There are three important functions in it `SurfTrack::improve_mesh`, `SurfTrack::topology_changes` and `SurfTrack::integrate`. The first two methods are combined under the name `stabilize` whereas the latter one, which modifies tangential velocities to avoid vertex collisions, is not wrapped [^1].

# ```@autodocs
# Modules = [ElTopo]
# ```

# Let's first load ElTopo library and load the triangulation to stabilize. In this case, a sphere made out of two subdivisions of an icosahedron.

using ElTopo
include("sphere.jl")
msh = unitsphere(2)

# Berore executing stabilization method we need to initialize parameters. See docstrings for fields and also [article](http://www.cs.ubc.ca/labs/imager/tr/2009/eltopo/sisc2009.pdf).

scale = 0.2
par = SurfTrack(
    min_edge_length = 0.7*scale,
    max_edge_length = 1.5*scale,
    max_volume_change = 0.1*scale^3,
    merge_proximity_epsilon = 0.5*scale,
)

# Then putting parameters and mesh in we can execute stabilization:

newmsh = stabilize(msh,par)

# [^1]: I had issues with using Cxx to represent `SurfTrack` instance in julia. Thus the best I was able to do was to use a single `icxx` block in julia function to do all the stuff. In spite of that, I think that covers most of the use cases.