text:
#  Tests of the functions in Entropy.jl
using Test
using DataFrames
using Main.CGP

P = [1,1,1,2];
Q = [2,2,3,3];
PQ = [1,1,1,1,2,2,3,3];
if Main.CGP.MyInt == UInt8
  R = [0x01,0x01,0x01,0x01,0x02,0x02,0x03,0x03];
elseif Main.CGP.MyInt == UInt16
  R = [0x001,0x001,0x001,0x001,0x002,0x002,0x003,0x003];
end
S = ["A", "B", "C", "C"]
F1 = [0.2, 0.4, 0.1, 0.3]
F2 = [0.5, 0.2, 0.1, 0.2]
@test pop_to_dist(P) == Dict(1=>0.75,2=>0.25)  # Note 0.25 and 0.75 are represented exactly
@test pop_to_dist(Q) == Dict(2=>0.5,3=>0.5)
@test pop_to_dist(PQ) == Dict(1=>0.5,2=>0.25,3=>0.25)
@test pop_to_dist(S) == Dict("A"=>0.25, "B"=>0.25, "C"=> .5 )
@test pop_to_dist(R) == Dict( R[6]=>0.25, R[8]=>0.25, R[1]=> .5 )
@test pops_to_dist(P,Q) == Dict((1,2)=>0.5, (2,3)=>0.25, (1,3)=> .25 )
@test isapprox( pops_to_tbl([P,Q]), [0.375 0.125 0.0; 0.0 0.25 0.25] )
df = DataFrame(names1=["A","B","C"], counts1=[2,5,3], names2=["A","B","D"], counts2=[4,5,1])
@test isapprox( pop_counts_to_tbl(df), [0.1 0.25 0.15 0.0; 0.2 0.25 0.0 0.05])
tbl = [0.125 0.25 0.125 0.0; 0.125 0.25 0.0 0.125]
@test table_row_to_dist(tbl,1) == Dict(1=>0.25,2=>0.5,3=>0.25,4=>0.0)
@test column_marginal(tbl) == [0.25, 0.5, 0.125, 0.125]
@test isapprox( CGP.entropy(PQ), 1.5 )
tbl = [0.25 0.5 0.125 0.125; 0.375 0.25 0.0 0.125]
@test isapprox( CGP.entropy(tbl,1), 1.75 )
@test isapprox( CGP.entropy(tbl,2), 1.4591479170272448 )
@test relative_entropy(pop_to_dist(Q),pop_to_dist(PQ)) == 1.0
@test relative_entropy(Q,PQ) == 1.0
@test isapprox( CGP.entropy(F1), 1.8464393446710154 )
@test isapprox( CGP.entropy(F2), 1.7609640474436812 )
@test isapprox( mutual_information(F1,F2), 0.40563906222956625 )

# Test the Cover and Thomas (1991) implementations of joint, conditional, relative CGP.entropy; mutual information
# Example 2.2.1 page 17 of Cover and Joy
tbl = [ 1/8 1/16 1/32 1/32; 1/16 1/8 1/32 1/32; 1/16 1/16 1/16 1/16; 1/4 0 0 0]
trans_tbl = convert(Array{Float64,2},transpose(tbl))
@test joint_entropy(tbl) ≈ 27/8  # H(X,Y)
X = column_marginal(tbl)
Y = row_marginal(tbl)
@test CGP.entropy(X) ≈ 7/4
@test CGP.entropy(Y) ≈ 2.0
@test conditional_entropy(tbl) ≈ 11/8   # H(X|Y)
@test conditional_entropy(trans_tbl) ≈ 13/8  # H(Y|X)
r = rand()
s = rand()
#p = Dict{Int64,Float64}( 0=>1.0-r, 1=>r )
#q = Dict{Int64,Float64}( 0=>1.0-q, 1=>q )
p = Main.CGP.DIST_TYPE( 0=>1.0-r, 1=>r )
q = Main.CGP.DIST_TYPE( 0=>1.0-s, 1=>s )
@test relative_entropy(p,q) ≈ (1-r)*log2((1-r)/(1-s)) + r*log2(r/s)  # Eq 2.31
@test relative_entropy(q,p) ≈ (1-s)*log2((1-s)/(1-r)) + s*log2(s/r)  # Eq 2.32
@test mutual_information(tbl) ≈ 3/8
@test mutual_information(trans_tbl) ≈ 3/8
@test mutual_information(tbl) ≈ CGP.entropy(X)-conditional_entropy(tbl)  # I(X;Y)=H(X)-H(X|Y) eq 2.43
@test mutual_information(tbl) ≈ CGP.entropy(Y)-conditional_entropy(trans_tbl)  # I(X;Y)=H(Y)-H(Y|X) eq 2.44
@test mutual_information(tbl) ≈ CGP.entropy(X)+CGP.entropy(Y)-joint_entropy(tbl)  # I(X;Y)=H(X)+H(Y)-H(X,Y0 eq 2.45
# Now test relationship of Cover mutual information to Sherwin mutual information
@test mutual_information(tbl) ≈ mutual_information(trans_tbl)  # mutual information is symmetric
@test mutual_information(tbl) ≈ sherwin_mutual_information(tbl)
@test mutual_information(trans_tbl) ≈ sherwin_mutual_information(trans_tbl)


text:
#counter=5;

load("utils.jl")

is_expr(head, ex) = (isa(ex, Expr) && ex.head == head)

# == Context ==================================================================

type SymEntry
    value::Any
    kind::Symbol
end

type Context
    symbols::HashTable{Symbol,SymEntry}  # current symbol bindings
    props::HashTable{Symbol,Any}         # context properties
    emit_function::Function              # sink for emitted instructions

    function Context(symbols::HashTable{Symbol,SymEntry}, props_seq, emit)
        # copy props (required!) and make sure it's the right type
        props = makedict(Symbol,Any,props_seq...)

        # convenience access through c.lhs, c.scalar; defalts = false, true
        lhs = setdefault(props, :lhs, false)::Bool
        scalar = setdefault(props, :scalar, true)::Bool

        new(symbols, props, emit, Symbol[], Symbol[], lhs, scalar)
    end

    inputs::Vector{Symbol}
    outputs::Vector{Symbol}

    # convenience access to props[:lhs], props[:scalar]
    lhs::Bool
    scalar::Bool
end
Context(c::Context) = Context(c.symbols, (c.props,), c.emit_function)
function Context(c::Context, props_seq...)
    Context(c.symbols, {c.props, props_seq...}, c.emit_function)
end

function child(c::Context, new_props)
# Create a child context from c:
#     Copy c
#     Replace properties given by makedict{Symbol}{Any}(new_props)
#     Set lhs=false, scalar=true if not explicitly given

#     # copy c
#     c = Context(c)
#     # modifies props: only at creation!
#     c.props[:lhs] = false
#     c.props[:scalar] = true
#     update(c.props, new_props)
#     return c
    c = Context(c, (@dict lhs=false scalar=true), new_props)
end
child(c::Context) = child(c::Context, ())

function create_argument(c::Context, name::Symbol)
    if has(context.props, :verbose)
        println("create_argument:")
        println("\tcontext = $context")
        println("\tname = $name")
        println()
    end

    if context.scalar; error("unbound scalar $name"); end

    arg_name = name
    if c.lhs # create output
        c.symbols[name] = SymEntry(arg_name, :output)
        append!(context.outputs, {arg_name})
    else     # create input
        c.symbols[name] = SymEntry(arg_name, :input)
        append!(context.inputs, {arg_name})
    end
    return arg_name
end

# -- Flatten ------------------------------------------------------------------

flatten(c::Context, ex::Any) = ex
flatten(c::Context, exprs::Vector) = { flatten(c, ex) | ex in exprs }

# like flatten(context, exprs::Vector), but only returns last value
function flatten_block(c::Context, exprs::Vector)
    value = nothing
    for ex in exprs
        value = flatten(c, ex)
    end
    value
end

function flatten(context::Context, sym::Symbol)
    if has(context.props, :verbose)
        println("flatten symbol:")
        println("\tcontext = $context")
        println("\tcontext.symbols = $(context.symbols)")
        println("\tsym = $sym")
        println()
    end

    if context.lhs
        if context.scalar # sym is an intermediate
            return sym
        else              # sym is an output
            return create_argument(context, sym)
        end
    else # rhs
        println("rhs")
        if has(context.symbols, sym)
            println("\thas\n")
            return context[sym].value             # return current symbol value
        else
            if context.scalar
                println("\tscalar:flatten\n")
                return flatten(context, expr(:call, :readinput, sym))
            else
                println("\tcreate_argument\n")
                println("\tcontext = $context")
                # create new argument named sym
                return create_argument(context, sym)
            end
        end
    end
end

function flatten(context::Context, ex::Expr)
    if has(context.props, :verbose)
        println("flatten:")
        println("\tcontext = $context")
        println("\tex = $ex")
        println()
    end

#     global counter
#     counter -= 1
#     if counter <= 0; error(); end

    if context.lhs
        if ex.head == :ref
            # indexed assignment
            dest = flatten(ex.args[1], (@dict lhs=true scalar=false))
            inds = flatten(ex.args[2:end])
            return expr(:ref, dest, inds...)
        end
        error("Unimplemented as LHS: ex = :$(ex)")
    else   # rhs Expr
        if ex.head == :block
            return flatten_block(context, ex.args)
        elseif ex.head == :line # ignore line numbers
            return nothing # CONSIDER: could this shadow the last actual value?
        elseif ex.head == :(=)  # assignment
            lhs = ex.args[1]
            rhs = flatten(context, ex.args[2])

            if has(context.props, :verbose)
                println("assignment:")
                println("\tlhs = $lhs")
                println("\trhs = $rhs")
            end

            if isa(lhs, Symbol)
                # straight assignment: just update symbol table
                context.symbols[lhs] = context.symbols[rhs]
                return rhs
            elseif (isa(lhs, Expr)) && (lhs.head == :ref)
                # indexed assignment
                lhs = flatten(child(context,(@dict lhs=true)), lhs)
                emit(context, expr(:(=), lhs, rhs))
                return rhs
            end
            error("Unimplemented: lhs = ($lhs)")
        elseif (ex.head == :call) || (ex.head == :ref)
            node = gensym()
            emit(context, expr(:(=), node, flatten_invocation(context, ex)))
            return node
        else
            # error("Unimplemented expr type: head = :$(ex.head)")
            error("Unimplemented: ex = :$(ex)")
        end
    end
end

function flatten_invocation(context::Context, ex::Expr)
    if has(context.props, :verbose)
        println("flatten_invocation:")
        println("\tcontext = $context")
        println("\tex = $ex")
        println()
    end

#     global flc = context
#     global fle = ex

    if ex.head == :call
        arg1 = ex.args[1]
        if arg1 == :readinput
            arg2 = flatten(child(context,(@dict scalar=false)), ex.args[2])
            @assert length(ex.args) == 2
            return expr(ex.head, arg1, arg2)
        else
            argsrest =  flatten(context, ex.args[2:end])
            return expr(ex.head, arg1, argsrest...)
        end
    else
        arg1 = flatten(child(context,(@dict scalar=false)), ex.args[1])
        argsrest =  flatten(context, ex.args[n_ns:end])
        return expr(ex.head, arg1, argsrest...)
    end
end

# -- Test code ----------------------------------------------------------------

code = quote
    A = B.*C + D[j,i]
    dest[2i, 2j] = A
end

symbols = HashTable{Symbol,SymEntry}()
#props = ()
props = @dict verbose=true
out_code = {}
emit = ex->append!(out_code, {ex})
context = Context(symbols, (props,), emit)

value = flatten(context, code)

ext: jl

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
# Performance tuning

v = [[1, 2, 3], [1, 2, 3], [1, 2, 3], [2, 3, 5], [1, 3, 4], [1, 2, 5], [2, 3, 4], [1, 4, 5], [3, 4, 5]]
v = [[1, 2], [1, 3], [2, 4], [1, 2, 3], [1, 2, 4], [1, 3, 4], [1, 2, 3, 4], [1, 2, 3, 5], [2, 3, 4, 6]]
# v = [rand([1, 2, 3, 4, 5], 10) for x = 1:1000];
# @code_warntype freq_itemsets(v, 0.5)

# testing freq_itemsets()
v = [[1, 2, 3], [1, 2, 3], [1, 2, 3],  [1, 2, 5], [1, 3, 4], [1, 4, 5], [2, 3, 4], [2, 3, 4], [2, 3, 5], [3, 4, 5]]
freq_itemsets(v, 0.2)

# testing gen_rules()
rules = Vector{Rule}(0)
freq = [1, 2]
consq = [Int[], [1], [2], [3], [4], [5]]
trans = [[1, 2], [1, 3], [2, 4], [1, 2, 3], [1, 2, 4], [1, 3, 4], [1, 2, 3, 4], [1, 2, 3, 5], [2, 3, 4, 6]]
ap_genrules!(freq, consq, trans, 0.20, 0.01, rules)
rules

rules = Vector{Rule}(0)
freq = [1, 2, 3, 4]
consq = [[1], [2], [3], [4]]
trans = [[1], [2], [1, 2], [1, 3], [2, 4], [1, 2, 3], [1, 2, 4], [1, 3, 4], [1, 2, 3, 4], [1, 2, 3, 5], [2, 3, 4, 6]]
ap_genrules!(freq, consq, trans, 0.2, 0.01, rules)
rules

rule1 = Vector{Rule}(0)
@time gen_onerules!([1], [[1], [2], [3]], rule1, [[1, 2], [2, 3], [1, 3]], 0.2, 0.01)

v = [rand([1, 2, 3, 4, 5, 6], 10) for x = 1:100];
rules = apriori(v, 0.01, 0.08)
show_rulestats(rules)

# Example using UCI Adult dataset

using DataFrames

d = readtable("./data/adult.csv", nastrings = ["?"])
complete_cases!(d)

dmat = convert(Array, d)
transactions = make_transactions(dmat)
@time rules = apriori(transactions[1:1000], 0.60, 0.8);

display(rules)


text:
module EventTimeSeries

export
    # types
    Event,
    EventTS,
    # extract
    timestamps,
    tags,
    values,
    # combine
    split,
    splice,
    merge_tags,
    drop_repeated,
    pretty_print

include("helpers.jl")
include("tagtype.jl")
include("timeseries.jl")

end


text:
# This example is used to generate plots in the paper, "Reflection from a multi-species material and its transmitted effective wavenumber." Proc. R. Soc. (2018): 20170864.

# Everything related to ploting has been commented so that this example can run with requiring the Plots package. That way, this example can be checked everytime the package is modified.

using EffectiveWaves

# using LaTeXStrings
# using Plots
# height=500
 # pyplot(linewidth=3, size=(2.6*height,height), border=false)

 # Plots.scalefontsizes(1.7)

## choose material

  # concrete
  medium = LimeStone
  inclusion1 = AirDry
  inclusion2 = Brick

  ωfactor = 50.0;
  ωs = 20*ωfactor.*LinRange(real(medium.c/10000),real(medium.c),400) # k from 0 to 1

  volfrac = 0.16
  r1 = 0.1/ωfactor; vol1 = 0.06
  r2 = 1.0/ωfactor; vol2 = volfrac - vol1

  sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = vol1)
  sp2 = Specie(ρ=inclusion2.ρ ,r=r2, c=inclusion2.c, volfrac = volfrac-vol1)

  # True wavenumber
  kTs  = wavenumber_low_volfrac(ωs, medium, [sp1,sp2]; tol=0.5e-5);
  # Approximate challis
  kTCs = wavenumber_challis(ωs, medium, [sp1,sp2]; hankel_order=10);

  kTs_arr = [kTs,kTCs];
  speed_arr = [ ωs./real(ks) for ks in kTs_arr];
  push!(speed_arr, 0 .*ωs .+ real(medium.c))
  atten_arr = imag.(kTs_arr)

  styles = [:solid :dashdot :dot]
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*C}$" "Lime Stone"]
  ys_arr = speed_arr;
  xs = r1.*(ωs./real(medium.c));
  m =5;
  y1 = min(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  y2 = max(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  # p1 = plot(xs, ys_arr, xlabel=L"a_S k", ylabel="sound speed (m/s)", labels=labs, line = styles)

  ys_arr = atten_arr;
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*C}$"]
  # p2 = plot(xs, ys_arr, labels=labs, xlabel=L"a_S k", ylabel="attenuation (1/m)", line=styles)

  # plot(p1,p2)
  # try mkdir("media") end
  # savefig("media/compare_concrete_large-w.png")
  # savefig("media/compare_concrete_large-w.pdf")
  # gui()

# Plots.scalefontsizes(1/1.7)


text:
using DThree

## Make scatter plot
x = linspace(0, 2*pi, 50)
y = cos(x)
dataset = "[" * join(["[$x,$y]" for (x,y) in zip(x,y)], ", ") * "]"
    width, height = 600, 400
    padding = 30

    p = D3Plot()
    d3 =D3()

    p * d3.var("svg").select("body").append("svg").attr("width", width).attr("height", height)
    ## scales
    p * d3.var("xScale").scale_linear().domain([min(x),max(x)]*1.1).range([0, width])
    p* d3.var("yScale").scale_linear().domain([min(y),max(y)]*1.1).range([height, 0])
    ## axes
    p * d3.var("xAxis").svg_axis().scale(asis("xScale")).orient("bottom")
    p * d3.receiver("svg").append("g").attr("class", "axis").attr("transform", "translate(0,$(height - padding))").call(asis("xAxis"))
    p * d3.var("yAxis").svg_axis().scale(asis("yScale")).orient("left").ticks(5)
    p * d3.receiver("svg").append("g").attr("class", "axis").attr("transform", "translate($padding, 0)").call(asis("yAxis"))

    p * d3.receiver("svg").selectAll("circle").data(asis(dataset)).enter().append("circle").
      attr("cx", asis("function(d) {return xScale(d[0])}")).
      attr("cy", asis("function(d) {return yScale(d[1])}")).
      attr("r", 5)

style = "
.axis path,
.axis line {
    fill: none;
    stroke: black;
    shape-rendering: crispEdges;
}

.axis text {
    font-family: sans-serif;
    font-size: 11px;
}
"

    browse(p, style=style)


text:
# fix eta k cchoice
# fix quadstep?
"""
  HZAW(;)

Implements the approximate Wolfe line search developed in [HZ2005].

theta controls the search in U3 of the paper. It defaults to 0.5, and as
a result uses bisection when (what is U3?)

We tweak the original algorithm slightly, by backtracking into a feasible
region if the original step length results in function values that are not
finite.
"""
struct HZAW{T} <: LineSearcher
    decrease::T
    curvature::T
    θ::T
    γ::T
end
Base.summary(::HZAW) = "Approximate Wolfe Line Search (Hager & Zhang)"
function HZAW(; decrease = 0.1, curvature = 0.9, theta = 0.5, gamma = 2 / 3)
    if !(0 < decrease ≤ curvature)
        println(
            "Decrease constant must be positive and smaller than the curvature condition.",
        )
    end
    if !(curvature < 1)
        println("Curvature constant must be smaller than one.")
    end
    HZAW(decrease, curvature, theta, gamma)
end

function find_steplength(mstyle, hzl::HZAW, φ, c, ϵk = 1e-6; maxiter = 100)
    # c = initial(k) but this is done outisde
    T = typeof(φ.φ0)
    ϵk = T(ϵk)
    δ = T(hzl.decrease)
    σ = T(hzl.curvature)
    ρ = T(5)
    φ0, dφ0 = φ.φ0, φ.dφ0
    φc, dφc = φ(c, true)

    # Backtrack into feasible region; not part of original algorithm
    ctmp, c = c, c
    iter = 0
    while !isfinite(φc) && iter <= maxiter
        iter += 1
        # don't use interpolation, this is vanilla backtracking
        ctmp, c, φc = interpolate(FixedInterp(), φ, φ0, dφ0, c, φc, T(1) / 10)
    end

    # initial convergence
    # Wolfe conditions
    if δ * dφ0 ≥ (φc - φ0) / c && dφc ≥ σ * dφ0
        return c, φc, true
    end
    # Approximate Wolfe conditions
    if (2 * δ - 1) * dφ0 ≥ dφc ≥ σ * dφ0 && φc ≤ φ0 + ϵk
        return c, φc, true
    end
    # Set up interval
    a0, b0 = bracket(hzl, c, φ, ϵk, ρ)
    j = 0
    aj, bj = a0, b0
    # Main loop
    while j < 50
        a, b = secant²(hzl, φ, aj, bj, ϵk)
        if b - a > hzl.γ * (bj - aj)
            c = (a + b) / 2
            φc, dφc = φ(c, true)
            a, b = update(hzl, a, b, c, φ, φc, dφc, ϵk)
        end

        aj, bj = a, b
        j += 1
        if _wolfe(φ0, dφ0, c, φc, dφc, δ, σ, ϵk) ||
           _approx_wolfe(φ0, dφ0, c, φc, dφc, δ, σ, ϵk)
            return c, φc, true
        end
    end
    return T(NaN), T(NaN), false
end
_wolfe(φ0, dφ0, c, φc, dφc, δ, σ, ϵk) = δ * dφ0 ≥ (φc - φ0) / c && dφc ≥ σ * dφ0
_approx_wolfe(φ0, dφ0, c, φc, dφc, δ, σ, ϵk) =
    (2 * δ - 1) * dφ0 ≥ dφc ≥ σ * dφ0 && φc ≤ φ0 + ϵk
"""
   _U3

Used to take step U3 of the updating procedure [HZ, p.123]. The other steps
are in update, but this step is separated out to be able to use it in
step B2 of bracket.
"""
function _U3(hzl::HZAW, φ, a::T, b::T, c::T, ϵk) where {T}
    # verified against paper description [p. 123, CG_DESCENT_851]
    φ0 = φ.φ0
    _a, _b = a, c
    # a)
    searching = true
    j = 1
    while searching && j < 50
        # convex combination of _a and _b; 0.5 implies bisection
        d = (1 - hzl.θ) * _a + hzl.θ * _b
        φd, dφd = φ(d, true)
        if dφd ≥ T(0) # found point of increasing objective; return with upper bound d
            _b = d
            return _a, _b
        else # now dφd < T(0)
            if φd ≤ φ0 + ϵk
                _a = d
            else # φ(d) ≥ φ0 + ϵk
                _b = d
            end
        end
        j += 1
    end
    _a, _b # throw error?
end

function update(hzl::HZ, a::T, b::T, c::T, φ, φc, dφc, ϵk) where {HZ<:HZAW,T}

    # verified against paper description [p. 123, CG_DESCENT_851]
    φ0 = φ.φ0
    #== U0 ==#
    if c ≤ a || c ≥ b # c ∉ (a, b)
        return a, b, (a = false, b = false)
    end
    #== U1 ==#
    if dφc ≥ T(0)
        return a, c, (a = false, b = true)
    else # dφc < T(0)
        #== U2 ==#
        if φc ≤ φ0 + ϵk
            return c, b, (a = true, b = false)
        end
        #== U3 ==#
        a, b = _U3(hzl, φ, a, b, c, ϵk)
        return a, b, (a = a == c, b = b == c)
    end
end
"""
  bracket

Find an interval satisfying the opposite slope condition [OSC] starting from
[0, c] [pp. 123-124, CG_DESCENT_851].
"""
function bracket(hzl::HZAW, c::T, φ, ϵk, ρ) where {T}
    # verified against paper description [pp. 123-124, CG_DESCENT_851]
    # Note, we know that dφ(0) < 0 since we're accepted that the current step is in a
    # direction of descent.
    φ0 = φ.φ0

    #== B0 ==#
    cj = c
    φcj, dφcj = φ(cj, true)
    # we only want to store a number, so we don't store all iterates
    ci, φi = T(0), φ0

    maxj = 100
    for j = 1:maxj
        #==================================================
          B1: φ is increasing at c, set b to cj as this is
              an upper bound, since φ is initially decrea-
              sing.
        ==================================================#
        if dφcj ≥ T(0)
            a, b = ci, cj
            return a, b
        else # dφcj < T(0)
            #== B2 : φ is decreasing at cj but function value is sufficiently larger than
            # φ0, use U3 to update. ==#
            if φcj > φ0 + ϵk
                a, b = _U3(hzl, φ, T(0), cj, c, ϵk)
                return a, b
            end
            #== B3 ==#
            # update ci instead of keeping all c's
            if φcj ≤ φ0 + ϵk
                ci = cj
                φci = φcj
            end
            # expand by factor ρ > 0 (shouldn't this be > 1?)
            cj = ρ * cj
            φcj, dφcj = φ(cj, true)
        end
    end
end

function secant(hzl::HZAW, a, dφa, b, dφb)
    # verified against paper description [p. 123, CG_DESCENT_851]
    #(a*dφb - b*dφa)/(dφb - dφa)
    # It has been observed that dφa can be very close to dφb,
    # so we avoid taking the difference
    a / (1 - dφa / dφb) + b / (1 - dφb / dφa)
end
function secant²(hzl::HZAW, φ, a, b, ϵk)
    # verified against paper description [p. 123, CG_DESCENT_851]
    #== S1 ==#
    φa, dφa = φ(a, true)
    φb, dφb = φ(b, true)
    c = secant(hzl, a, dφa, b, dφb)

    φc, dφc = φ(c, true)
    A, B, updates = update(hzl, a, b, c, φ, φc, dφc, ϵk)
    if updates.b # B == c
        #== S2: c is the upper bound ==#
        φB, dφB = φc, dφc
        _c = secant(hzl, b, dφb, B, dφB)
    elseif updates.a # A == c
        #== S3: c is the lower bound ==#
        φA, dφA = φc, dφc
        _c = secant(hzl, a, dφa, A, dφA)
    end
    updates
    if any(updates)
        #== S4.if: c was upper or lower bound ==#
        φ_c, dφ_c = φ(_c, true)
        _a, _b = update(hzl, A, B, _c, φ, φ_c, dφ_c, ϵk)
        return _a, _b
    else
        #== S4.otherwise: c was neither ==#
        return A, B
    end
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
using Pluto

notebook_name = length(ARGS) > 0 ? ARGS[1] : "dashboard.jl"

Pluto.run(;
    host="0.0.0.0",
    port=7777,
    launch_browser=false,
)

max_stars_repo_name: GPT-RL/generalization

max_issues_repo_name: GPT-RL/generalization

max_forks_repo_name: GPT-RL/generalization


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

robust geometric programming problem
given a convex set C in R_+^k (described by conic constraints) and B in R^{k, n},
calculate:
f(C, B) = sup_{c in C} (inf_{x in R^n, z in R_+^k} c'*z : B_i*x <= log(z_i), i = 1..k)
note the inner problem is an unconstrained GP, and C specifies coefficient uncertainty
for more details, see section 4.4 of:
"Relative entropy optimization and its applications"
(2017) by Chandrasekaran & Shah
the authors show that:
    f(C, B) = -inf_{c in C, v in R_+^k} (d(v, e*c) : B_j'*v = 0, j = 1..n)
where e = exp(1) and d(a, b) = ∑ᵢ a_i*log(a_i/b_i) is the relative entropy of a, b
=#

struct RobustGeomProgJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    n::Int
    k::Int
end

function build(inst::RobustGeomProgJuMP{T}) where {T <: Float64}
    (n, k) = (inst.n, inst.k)
    @assert n < k # want some degrees of freedom for v
    B = randn(T, k, n) # GP powers matrix
    model = JuMP.Model()
    JuMP.@variable(model, t)
    JuMP.@objective(model, Min, t)

    JuMP.@variable(model, c[1:k])
    JuMP.@variable(model, v[1:k])
    JuMP.@constraint(model, vcat(t, ℯ * c, v) in MOI.RelativeEntropyCone(1 + 2k))
    JuMP.@constraint(model, B' * v .== 0)

    # use bounded convex set C of R_+^k excluding origin
    # note that the entropy constraint already forces c >= 0
    JuMP.@constraint(model, vcat(sqrt(k) / 2, 1 .- c) in MOI.NormOneCone(1 + k))

    return model
end


text:
using FFTW,Plots;N=128;P=64N;dt=1/4N;NT=1024;W=200;w=W/P*N;
x=rand(P);v=collect(1:P.>P/2).*2 .-1;u()=(x.=mod.(x.+v/2*dt,1));n=zeros(N);
k=im.*2π*vcat(1,1:N/2,-N/2+1:-1);f(x)=Int(mod1(round(x*N),N));Σ=sum;
@gif for t in 1:NT;
  u();E=real.(ifft((E=fft((n.*=0;for j∈x;n[f(j)]+=w;end;n))./k;E[1]*=0;E)));
  u();v+=E[f.(x)]*dt;scatter(x,v,ylims=(-3,3),legend=false);
end every 8


text:
x = randobs(X[:,1:4])
@test size(x) == (4,)
@test any(X[:,i] == x for i in 1:4)

x = randobs(X[:,1:4], 2)
@test size(x) == (4, 2)
@test any(X[:,i] == x[:,1] for i in 1:4)
@test any(X[:,i] == x[:,2] for i in 1:4)


text:
module ForcesWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractGaussianContribution = @jimport org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution
const AbstractGaussianContributionContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContributionContext
const Atmosphere = @jimport org.orekit.models.earth.atmosphere.Atmosphere
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const AuxiliaryElements = @jimport org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements
const CelestialBody = @jimport org.orekit.bodies.CelestialBody
const Class = @jimport java.lang.Class
const DSSTAtmosphericDrag = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag
const DSSTForceModel = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel
const DSSTNewtonianAttraction = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction
const DSSTNewtonianAttractionContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttractionContext
const DSSTSolarRadiationPressure = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure
const DSSTTesseral = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral
const DSSTTesseralContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseralContext
const DSSTThirdBody = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody
const DSSTThirdBodyContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBodyContext
const DSSTZonal = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal
const DSSTZonalContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTZonalContext
const DragForce = @jimport org.orekit.forces.drag.DragForce
const DragSensitive = @jimport org.orekit.forces.drag.DragSensitive
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const ExtendedPVCoordinatesProvider = @jimport org.orekit.utils.ExtendedPVCoordinatesProvider
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldAbstractGaussianContributionContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.FieldAbstractGaussianContributionContext
const FieldAuxiliaryElements = @jimport org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements
const FieldDSSTNewtonianAttractionContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.FieldDSSTNewtonianAttractionContext
const FieldDSSTTesseralContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.FieldDSSTTesseralContext
const FieldDSSTThirdBodyContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.FieldDSSTThirdBodyContext
const FieldDSSTZonalContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.FieldDSSTZonalContext
const FieldEventDetector = @jimport org.orekit.propagation.events.FieldEventDetector
const FieldForceModelContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.FieldForceModelContext
const FieldOrbit = @jimport org.orekit.orbits.FieldOrbit
const FieldShortPeriodTerms = @jimport org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const ForceModelContext = @jimport org.orekit.propagation.semianalytical.dsst.forces.ForceModelContext
const Frame = @jimport org.orekit.frames.Frame
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const Object = @jimport java.lang.Object
const Orbit = @jimport org.orekit.orbits.Orbit
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const PropagationType = @jimport org.orekit.propagation.PropagationType
const RadiationSensitive = @jimport org.orekit.forces.radiation.RadiationSensitive
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const Set = @jimport java.util.Set
const ShortPeriodTerms = @jimport org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const UnnormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider

export AbstractGaussianContribution
export AbstractGaussianContributionContext
export DSSTAtmosphericDrag
export DSSTForceModel
export DSSTNewtonianAttraction
export DSSTNewtonianAttractionContext
export DSSTSolarRadiationPressure
export DSSTTesseral
export DSSTTesseralContext
export DSSTThirdBody
export DSSTThirdBodyContext
export DSSTZonal
export DSSTZonalContext
export FieldAbstractGaussianContributionContext
export FieldDSSTNewtonianAttractionContext
export FieldDSSTTesseralContext
export FieldDSSTThirdBodyContext
export FieldDSSTZonalContext
export FieldForceModelContext
export FieldShortPeriodTerms
export ForceModelContext
export ShortPeriodTerms
export equals
export get_a
export get_alpha
export get_ao_r3_pow
export get_atmosphere
export get_auxiliary_elements
export get_ax2o_a
export get_bb
export get_bbb
export get_beta
export get_bo_a
export get_bo_a_bpo
export get_body
export get_chi
export get_chi2
export get_class
export get_co2ab
export get_coefficients
export get_coefficients_key_prefix
export get_cxo2n2a2
export get_drag
export get_e2
export get_equatorial_radius
export get_events_detectors
export get_field_auxiliary_elements
export get_field_events_detectors
export get_gamma
export get_gm
export get_hk
export get_hxxx
export get_k2mh2
export get_k2mh2o2
export get_kxxx
export get_m2ao_a
export get_m_co2ab
export get_max_ar3_pow
export get_max_ecc_pow
export get_max_freq_f
export get_mean_element_rate
export get_mean_motion
export get_moa
export get_mu
export get_muo_r3
export get_muoa
export get_oo_ab
export get_oo_bpo
export get_oo_mu
export get_ooa
export get_ooab
export get_oon2a2
export get_orbit_period
export get_parameters
export get_parameters_drivers
export get_provider
export get_qns
export get_ratio
export get_rbar
export get_res_orders
export get_roa
export get_spacecraft
export get_theta
export get_ton2a
export get_x
export get_x2on2a2xp1
export get_x3on2a
export get_xon2a2
export get_xx
export get_xxx
export getb
export hash_code
export initialize
export notify
export notify_all
export register_attitude_provider
export to_string
export update_short_period_terms
export value
export wait

include("abstract_gaussian_contribution.jl")
include("abstract_gaussian_contribution_context.jl")
include("dsst_atmospheric_drag.jl")
include("dsst_force_model.jl")
include("dsst_newtonian_attraction.jl")
include("dsst_newtonian_attraction_context.jl")
include("dsst_solar_radiation_pressure.jl")
include("dsst_tesseral.jl")
include("dsst_tesseral_context.jl")
include("dsst_third_body.jl")
include("dsst_third_body_context.jl")
include("dsst_zonal.jl")
include("dsst_zonal_context.jl")
include("field_abstract_gaussian_contribution_context.jl")
include("field_dsst_newtonian_attraction_context.jl")
include("field_dsst_tesseral_context.jl")
include("field_dsst_third_body_context.jl")
include("field_dsst_zonal_context.jl")
include("field_force_model_context.jl")
include("field_short_period_terms.jl")
include("force_model_context.jl")
include("short_period_terms.jl")

end


text:
module RegularizationDP

import Optim
import Zygote

function optimize(forward_model, make_realization, regularization, x0, options)
	f = x->regularization(x) + forward_model(make_realization(x))
	opt = Optim.optimize(f, x->Zygote.gradient(f, x)[1], x0, Optim.LBFGS(), options; inplace=false)
	return make_realization(opt.minimizer), opt
end

end


text:
"""
This file demonstrates most functionality of the package.
Load Revise if you want to make changes to the source files.
"""

using Pkg; Pkg.activate(pwd())
using Revise

@time using TumorGrowth

#########################
######## SANDBOX ########
#########################

#########################
####### SIMULATION ######
#########################

@time simoutput = birth_death_pushing(3000; b=1., d=0.0, μ=0.1, ρ=Inf, dim=2, seed=1010, showprogress=false)
tumor = simoutput[:tumor]
mutations = simoutput[:mutations]

@time simoutput = nonspatial(1000000; b=1., d=0.0, μ=0.3, showprogress=true)

bumor = deepcopy(tumor)

birth_death_pushing!(bumor, length(tumor)+1; b=0.69, d=0.0, μ=0.3, dim=2,
 t=simoutput[:time], cur_id=simoutput[:index], cur_mutation=simoutput[:mutation])

birth_death_pushing(12.; b=0.69, d=0.0, μ=0.3, dim=3, seed=1010)

tumordf = tumor |> DataFrame
mutationsdf = mutations |> DataFrame

#########################
##### TUMOUR SAVING #####
#########################

using DataFrames, CSV

tumor = birth_death_pushing(5000; b=1., d=0.0, μ=0.3, ρ=6., dim=3, seed=1234)[:tumor]
CSV.write("examples/test_5000_3d.csv", DataFrame(tumor), delim='\t')
b = data_import("examples/test_5000_3d.csv")

########################
####### PLOTTING #######
########################

using Plots: palette

b = data_import("examples/test_2000_2d.csv")
@time plotting(b)
@time plotting_colored_mutations(b, colorpalette = palette(:tab20), shading=false, inline=false, show_warning = true, path="examples/test_2000_2d.png")

b = data_import("examples/test_5000_3d.csv")
@time plotting(b)
@time plotting_colored_mutations(b, shading=true, path="examples/test_5000_3d.png", size=(500,500))

# translating and layouts

using GLMakie

tumor1 = birth_death_pushing(10000; b=1.,d=0., μ=0.2, ρ=7., dim=3)[:tumor] |> DataFrame
tumor2 = birth_death_pushing(10000; b=1.,d=0., μ=0.2, ρ=Inf, dim=3)[:tumor] |> DataFrame
tumor3 = birth_death_pushing(3000; b=1.,d=0., μ=0.2, ρ=1.4, dim=2)[:tumor] |> DataFrame

scene = Scene(resolution = (1000, 1000), show_axis=false)
t1 = plotting_colored_mutations!(scene, tumor1;
        shading=true, colorpalette=cgrad(:reds)[0.4:0.1:1.] )
t2 = plotting_colored_mutations!(scene, tumor2;
        shading=true, colorpalette=cgrad(:blues)[0.4:0.1:1.] )[end]
t3 = plotting_colored_mutations!(scene, tumor3;
        shading=false, colorpalette=cgrad(:greens)[0.2:0.1:1.] )[end]
translate!(t2, Vec3f0(100, 0, 0))
translate!(t3, Vec3f0(50, 100, 0))
scene |> display

using GLMakie.MakieLayout

scene, layout = layoutscene(show_axis=false)
layout[1, 1:2]  = ax1 = LScene(scene, camera = cam2d!, resolution=(500,500))
layout[2, 1]  = ax2 = LScene(scene, camera = cam3d!, raw = false, show_axis=false, resolution=(500,500))
layout[2, 2] = ax3 = LScene(scene, camera = cam3d!, raw = false, show_axis=false, resolution=(500,500))

t1 = plotting_colored_mutations!(ax2, tumor1;
        shading=true, colorpalette=cgrad(:reds)[0.4:0.1:1.] )
t2 = plotting_colored_mutations!(ax3, tumor2;
        shading=true, colorpalette=cgrad(:blues)[0.4:0.1:1.] )
# b = deepcopy(tumor3); b.position ./= 1.8
t3 = plotting_colored_mutations!(ax1, tumor3; markersize=1.1,
        shading=false, colorpalette=cgrad(:greens)[0.2:0.1:1.] )

scene |> display

# plotting sampletumors

@time tumor = birth_death_pushing(3000; b=1., d=0.0, μ=0.5, ρ=1.6, dim=2, showprogress=false)[:tumor] |> DataFrame
samples, sampletumor = multi_region_sequencing(tumor, n=50)

plotting_colored_mutations(sampletumor; markersize=sampletumor.sample_r, show_axis=true)
plotting_sampletumor_pies(sampletumor, first=30, colorpalette = TumorGrowth.distinguishable_colors(30))
plotting_sampletumor_pies(sampletumor, mutations = [1,10,55])

########################
####### Sampling #######
########################

b = data_import("examples/test_2000_2d.csv")

fig = b |> plotting_colored_mutations
cross_section(b; x=10., width=3.) |> data -> plotting!(fig, data; color=:black)
cross_section(b; y=25., width=3.) |> data -> plotting!(fig, data; color=:black)
radial_sample(b; r=30., width=3.) |> data -> plotting!(fig, data; color=:black, path="examples/test_2000_2d.png")

fig = b |> plotting_colored_mutations
bulk(b; pos=(25.,0.), box=(20,30)) |> data -> plotting!(fig, data; color=:black)
punch(b; pos=(-20,-20), r=10) |> data -> plotting!(fig, data; color=:black)

b = data_import("examples/test_5000_3d.csv")
scene = b |> plotting_colored_mutations
cross_section(b; x=10., width=3.) |> data -> plotting!(scene, data; shading=true, color=:black, path="examples/test_5000_3d.png")
cross_section(b; y=5., width=3.) |> data -> plotting!(scene, data; color=:black)

display(scene)

fig = cross_section(b; x=10., width=6.) |> plotting_colored_mutations
radial_sample(b; r=20., width=3.) |> data->cross_section(data; x=10., width=6.) |> data -> plotting!(fig, data; color=:black)

########################
####### Analysis #######
########################

using Plots, StatsBase, Statistics

params = (b=1., d=0.2, μ=0.3, ρ=Inf, dim=3)
tumor = birth_death_pushing(4000; params...)[:tumor] |> DataFrame

tumor |> mutation_freqs
tumor |> stochastic_sequencing

function M!(fig, f; res=0.0, nBins = 100, normalize=false, plotargs...)
        hist = fit(Histogram, 1 ./ f[res.<f], nbins=nBins, closed=:left)
        Mcounts = [ sum(hist.weights[1:n]) for n=1:length(hist.weights)]
        Mcounts = normalize ? Mcounts ./ last(Mcounts) : Mcounts
        Plots.plot!(fig, midpoints(hist.edges[1]), Mcounts; plotargs...)
end
M(f; res=0.0, nBins = 100, normalize=false, plotargs...) = M!(Plots.plot(), f; res=res, nBins = nBins, normalize=normalize, plotargs...)

fig = M(mutation_freqs(tumor).frequency; res=0.0005, nBins=200,
        xlab="1/f", ylab="M(f)", lab="", marker=:o, ms=1.5)
β = params.μ * params.b/(params.b-params.d)
Plots.plot!(1:1000, n-> β*n, lab="")

f = filter(m -> m.reads > 1, stochastic_sequencing(tumor, readdepth=1000)).frequency
M!(fig, stochastic_sequencing(tumor, readdepth=2000).frequency; res=0.001, lab="", marker=:d, ms=1.5 )

########################
######## clones ########
########################

tumor = data_import("examples/test_2000_2d.csv")

fig = plotting(tumor; color=:lightblue, shading=false)
clone(tumor, nothing) |> t -> plotting!(fig, t, color=:red)
clone(tumor, 2) |> t -> plotting!(fig, t, color=:blue)

fig = plotting_colored_mutations(tumor; colorpalette=palette(:tab20), shading=false)
sort(by=size, rev=true, clones(tumor))[1] |> t -> plotting!(fig, t, color=:black)

fig = plotting(tumor; color=:lightgrey, shading=false)
sort(by=size, rev=true, clones(tumor))[1] |> clones .|> t -> plotting_colored_mutations!(fig, t; colorpalette=palette(:tab10), shading=false,autodepth=true)
fig

#########################
###### time series ######
#########################

using DataFrames, GLMakie
using Plots: palette, distinguishable_colors

N, simparams... = (N=5000, b=1., d=0.0, μ=0.3, ρ=6., dim=3, seed=1234)
tumor = birth_death_pushing(N; simparams...)[:tumor] |> DataFrame

time_series = tumor_stepper(range(0., last( tumor.t_birth), length=50); simparams...)

last(time_series) == tumor

record_growth(time_series; path="temp.gif",
        frames=1, shading=true, show_axis=false,
        colorpalette = palette(:tab20)
        )

#########################
###### density ~ b ######
#########################

using Plots

TumorGrowth.b_linear()
TumorGrowth.b_hill(6)
TumorGrowth.b_curve(1.; bup=1., ρc=1.1)

d, ρ = 0.2, 1.5
@time tumor = birth_death_pushing(5000; b=0.69, d=d, μ=0.3, ρ=ρ, dim=2)[:tumor]
tumordf = tumor |> DataFrame

# power law / surface growth
Plots.plot(log.(tumordf.t_birth), log.(1:length(tumor)), legend=:none )
Plots.plot!(log.(tumordf.t_birth), t -> 2*t, xlims=(0., 4.))

# exponential / bulk
Plots.plot(tumordf.t_birth, log.(1:5000), legend=:none )
Plots.plot!(tumordf.t_birth, t -> (log(2) - d)*t )

Plots.plot(-7.:0.01:7.,r-> TumorGrowth.w(r; σ=3.), fill=true)
Plots.plot!(-9.:0.01:9.,r-> TumorGrowth.w(r; σ=3.), xticks=-7:2:7, legend=:none)

TumorGrowth.b_linear()
TumorGrowth.b_hill(6)
Plots.plot(0.:0.01:9., ρ-> TumorGrowth.b_curve(ρ; bup=1., ρc=1.), legend=:none)
Plots.vline!([1.])

plotting_colored_mutations(tumordf, colorpalette = palette(:tab20), shading=false, inline=false)

Plots.histogram(tumordf.b, alpha=0.3, lw=0.2, lab="b")
Plots.vline!([median(tumordf.b)], lab="median", c=:blue)
Plots.vline!([d], lab="d", c=:black)

plotting(tumordf; color = tumordf.b./maximum(tumordf.b), colormap = cgrad(:reds), path="examples/surface_b.png")

#############################
### mulit region sampling ###
#############################

using Makie: meshscatter, meshscatter!, save, Point
using Statistics

@time tumor = birth_death_pushing(5000; b=0.69, d=0.1, μ=0.3, ρ=1.8, dim=2, seed=1010)[:tumor] |> DataFrame

lattice, samples = multi_region_sampling(tumor; n = 100, cells_per_sample = 20)

samples, sampletumor = multi_region_sequencing(tumor; n = 100, cells_per_sample = 20, stochastic = true)

mean(size.(samples,1))

scene = plotting_colored_mutations(tumor)
for sample in samples
    plotting!(scene, sample; color=:black)
end
display(scene)

sample_r = 2*sqrt(20/π)
cm = mean(tumor.position)
density = size(punch(tumor; pos=cm, r=10.), 1)/(π*10. ^2)
R = sqrt(size(tumor,1)/(π*density))
scene = plotting_colored_mutations(tumor)
meshscatter!(scene, [Point{2}.(lattice)...], markersize = sample_r, scale_plot = false)
meshscatter!(scene, 0:0.2:360 .|> ϕ-> Point{2}(R .*(cosd(ϕ),sind(ϕ)) .+ cm), markersize = 1. )
display(scene)

# save("multi_region_sampling.png", scene)

density = size(punch(tumor; pos = mean(tumor.position), r = 20.), 1)/(π*20^2)
scene = plotting_colored_mutations(sampletumor; markersize = sqrt(20/density/π), shading = true)
plotting_colored_mutations!(scene, tumor)

# save("multi_region_sequencing.png", scene)

##################################
## Reduce mutation artificially ##
##################################

b, d, μ = 0.69, 0.1, 1.0
@time tumor = birth_death_pushing(5000; b=b, d=d, μ=μ, ρ=Inf, dim=2)[:tumor] |> DataFrame

x = 0.3
tumor_reduced = reduced_μ(tumor, x)

using Plots, StatsBase

fig = tumor |> mutation_freqs |> df -> M(df.frequency; res = 0.001, lab="")
Plots.plot!(collect(0:1:1/0.001), finv -> μ*b/(b-d) *(finv - 1), lab="")
tumor_reduced |> mutation_freqs |> df -> M!(fig, df.frequency; res = 0.001, lab="")
Plots.plot!(collect(0:1:1/0.001), finv -> μ*x*b/(b-d) *(finv - 1), lab="")

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
module TestDataFrame
    using Base.Test
    using DataFrames

    #test_group("Operations on DataFrames that have column groupings")

    x = DataFrame(a = [1, 2, 3], b = [4, 5, 6])
    y = DataFrame(c = [1, 2, 3], d = [4, 5, 6])

    #
    # Equality
    #
    @test isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2]), b=@data([4, 5])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), c=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(b=@data([4, 5, 6]), a=@data([1, 2, 3])))
    @test !isequal(DataFrame(a=@data([1, 2, 2]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test isequal(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])))

    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2]), b=@data([4, 5])) != DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]), c=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(b=@data([4, 5, 6]), a=@data([1, 2, 3]))
    @test DataFrame(a=@data([1, 2, 2]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 3, NA]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6]))
    @test isna(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])))
    @test isna(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))

    z = deepcopy(x)

    z = hcat(x, y)

    v = DataFrame(
        a = [5,6,7],
        b = [8,9,10]
    )
    z = vcat(DataFrame[v, x])

    z = vcat(v, x)

    # Deleting columns removes any mention from groupings
    delete!(x, :a)
    @test names(x) == [:b]

    ## del calls ref, which properly deals with groupings
    z2 = z[:,[1,1,2]]
    @test names(z2) == [:a, :a_1, :b]

    #test_group("DataFrame assignment")

    #test_group("Empty DataFrame constructors")
    df = DataFrame(Int, 10, 3)
    @test size(df, 1) == 10
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Int}
    @test typeof(df[:, 3]) == DataVector{Int}

    df = DataFrame({Int, Float64, ASCIIString}, 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Float64}
    @test typeof(df[:, 3]) == DataVector{ASCIIString}

    df = DataFrame({Int, Float64, ASCIIString}, [:A, :B, :C], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Float64}
    @test typeof(df[:, 3]) == DataVector{ASCIIString}

    df = convert(DataFrame, zeros(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, ones(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, eye(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    #test_group("Other DataFrame constructors")
    df = DataFrame([{:a=>1, :b=>'c'}, {:a=>3, :b=>'d'}, {:a=>5}])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == DataVector{Int}
    @test typeof(df[:,:b]) == DataVector{Char}

    df = DataFrame([{:a=>1, :b=>'c'}, {:a=>3, :b=>'d'}, {:a=>5}], [:a, :b])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == DataVector{Int}
    @test typeof(df[:,:b]) == DataVector{Char}

    data = {:A => [1, 2], :C => [:1, :2], :B => [3, 4]}
    df = DataFrame(data)
    # Specify column_names
    df = DataFrame(data, [:C, :A, :B])

    # This assignment was missing before
    df = DataFrame(Column = [:A])
    df[1, :Column] = "Testing"

    # zero-row dataframe and subdataframe test
    df = DataFrame(x=[], y=[])
    @test nrow(df) == 0
    df = DataFrame(x=[1:3], y=[3:5])
    sdf = sub(df, df[:x] .== 4)
    @test size(sdf, 1) == 0

    @test hash(convert(DataFrame, [1 2; 3 4])) == hash(convert(DataFrame, [1 2; 3 4]))
    @test hash(convert(DataFrame, [1 2; 3 4])) != hash(convert(DataFrame, [1 3; 2 4]))
end


text:
using Revise
using LinearAlgebra, Plots, Setfield, Parameters

using PseudoArcLengthContinuation
const PALC = PseudoArcLengthContinuation

TY = Float64
AF = Array{TY}
#################################################################
if 1==1
	ENV["GKSwstype"] = "nul"    # needed for the GR backend on headless servers
	using Plots
	using SixelTerm
else
	using GR
	GR.inline("iterm")
end

using CuArrays
CuArrays.allowscalar(false)
# import LinearAlgebra: mul!, axpby!
# mul!(x::CuArray, y::CuArray, α::T) where {T <: Number} = (x .= α .* y)
# mul!(x::CuArray, y::T, α::CuArray) where {T <: Number} = (x .= α .* y)
# axpby!(a::T, X::CuArray, b::T, Y::CuArray) where {T <: Number} = (Y .= a .* X .+ b .* Y)

const TY = Float64
const AF = CuArray{TY}
#################################################################
# to simplify plotting of the solution
heatmapsol(x) = heatmap(reshape(Array(x), Nx, Ny)', color=:viridis)
norminf(x) = maximum(abs.(x))

Nx = 2^8
Ny = 2^8
lx = 8pi * 2
ly = 2*2pi/sqrt(3) * 2

X = -lx .+ 2lx/Nx * collect(0:Nx-1)
Y = -ly .+ 2ly/Ny * collect(0:Ny-1)

sol0 = [(cos(x) .+ cos(x/2) * cos(sqrt(3) * y/2) ) for x in X, y in Y]
		sol0 .= sol0 .- minimum(vec(sol0))
		sol0 ./= maximum(vec(sol0))
		sol0 = sol0 .- 0.25
		sol0 .*= 1.7
		# heatmap(sol0, color=:viridis)

using AbstractFFTs, FFTW, KrylovKit
import Base: *, \

# Making the linear operator a subtype of PALC.LinearSolver is handy as we will use it
# in the Newton iterations.
struct SHLinearOp{Treal, Tcomp, Tl1, Tplan, Tiplan} <: PALC.AbstractLinearSolver
	tmp_real::Treal         # temporary
	tmp_complex::Tcomp      # temporary
	l1::Tl1
	fftplan::Tplan
	ifftplan::Tiplan
end

struct SHEigOp{Tsh <: SHLinearOp} <: PALC.AbstractEigenSolver
	sh::Tsh
end

function SHLinearOp(Nx, lx, Ny, ly; AF = Array{TY})
	# AF is a type, it could be CuArray{TY} to run the following on GPU
	k1 = vcat(collect(0:Nx/2), collect(Nx/2+1:Nx-1) .- Nx)
	k2 = vcat(collect(0:Ny/2), collect(Ny/2+1:Ny-1) .- Ny)
	d2 = [(1-(pi/lx * kx)^2 - (pi/ly * ky)^2)^2 + 1. for kx in k1, ky in k2]
	tmpc = Complex.(AF(zeros(Nx, Ny)))
	return SHLinearOp(AF(zeros(Nx, Ny)), tmpc, AF(d2), plan_fft!(tmpc), plan_ifft!(tmpc))
end

function *(c::SHLinearOp, u)
	c.tmp_complex .= Complex.(u)
	c.fftplan * c.tmp_complex
	c.tmp_complex .= c.l1 .* c.tmp_complex
	c.ifftplan * c.tmp_complex
	c.tmp_real .= real.(c.tmp_complex)
	return copy(c.tmp_real)
end

function \(c::SHLinearOp, u)
	c.tmp_complex .= Complex.(u)
	c.fftplan * c.tmp_complex
	c.tmp_complex .=  c.tmp_complex ./ c.l1
	c.ifftplan * c.tmp_complex
	c.tmp_real .= real.(c.tmp_complex)
	return copy(c.tmp_real)
end

function (sh::SHLinearOp)(J, rhs; shift = 0., tol =  1e-9)
	u, l, ν = J
	udiag = l .+ 1 .+ 2ν .* u .- 3 .* u.^2 .- shift
	res, info = KrylovKit.linsolve( du -> -du .+ sh \ (udiag .* du), sh \ rhs, tol = tol, maxiter = 6)
	return res, true, info.numops
end

function (sheig::SHEigOp)(J, nev::Int; σ = 0.3)
	u, l, ν = J
	sh = sheig.sh
	udiag = l .+ 1 .+ 2ν .* u .- 3 .* u.^2
	N = size(u)

	A = du -> sh(J, du; shift = σ, tol = 1e-5)[1]

	# we adapt the krylov dimension as function of the requested eigenvalue number
	vals, vec, info = KrylovKit.eigsolve(A, AF(rand(eltype(u), N)), nev, :LM, tol = 1e-10, maxiter = 40, verbosity = 2, issymmetric = true, krylovdim = max(40, nev + 10))
	@show 1 ./vals .+ σ
	return 1 ./vals .+ σ, vec, true, info.numops
end

function F_shfft(u, p)
	@unpack l, ν, L = p
	return -(L * u) .+ ((l+1) .* u .+ ν .* u.^2 .- u.^3)
end

J_shfft(u, p) = (u, p.l, p.ν)

L = SHLinearOp(Nx, lx, Ny, ly, AF = AF)
Leig = SHEigOp(L) # for eigenvalues computation
# Leig((sol_hexa, -0.1, 1.3), 20; σ = 0.5)

par = (l = -0.1, ν = 1.3, L = L)

@time F_shfft(sol0, par); # 0.008022 seconds (12 allocations: 1.500 MiB)

opt_new = PALC.NewtonPar(verbose = true, tol = 1e-6, linsolver = L, eigsolver = Leig)
	sol_hexa, hist, flag = @time PALC.newton(
				F_shfft, J_shfft,
				AF(sol0), par,
				opt_new, normN = norminf)
	println("--> norm(sol) = ", norminf(sol_hexa))

heatmapsol(sol_hexa)
####################################################################################################
# trial using IterativeSolvers

# function (sh::SHLinearOp)(J, rhs)
# 	u, l, ν = J
# 	udiag = l .+ 1 .+ 2ν .* u .- 3 .* u.^2
# 	res, info = res, info = KrylovKit.linsolve( u -> -u .+ sh \ (udiag .* u), sh \ rhs, tol = 1e-9, maxiter = 6)
# 	return res, true, info.numops
# end
####################################################################################################
deflationOp = DeflationOperator(2.0, (x, y)->dot(x, y), 1.0, [sol_hexa])

opt_new = @set opt_new.maxIter = 250
outdef, _, flag, _ = @time PALC.newton(
				F_shfft, J_shfft,
				0.4 .* sol_hexa .* AF([exp(-1(x+0lx)^2/25) for x in X, y in Y]),
				par,
				opt_new, deflationOp, normN = x-> maximum(abs.(x)))
		println("--> norm(sol) = ", norm(outdef))
		heatmapsol(outdef) |> display
		flag && push!(deflationOp, outdef)

####################################################################################################
opts_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.007, ds= -0.005, pMax = 0.2, pMin = -1.0, theta = 0.5, plotEveryNsteps = 5, newtonOptions = setproperties(opt_new; tol = 1e-6, maxIter = 15), maxSteps = 88,
	computeEigenValues = true,
	detectBifurcation = 0,
	precisionStability = 1e-5,
	saveEigenvectors = false,
	nev = 11 )

	br, u1 = @time PALC.continuation(
		F_shfft, J_shfft,
		deflationOp.roots[2], par, (@lens _.l),
		opts_cont, plot = true, verbosity = 2,
		plotSolution = (x, p;kwargs...)->heatmap!(reshape(Array(x), Nx, Ny)'; color=:viridis, kwargs...),
		printSolution = (x, p) -> norm(x), normC = norminf
		)

# branches = [br]
# push!(branches, br)
# plot(branches)
#
# plot(br)


text:
@testset "Models logfunctions" begin
    rng = Random.MersenneTwister(42)
    D = 5
    p0 = MvNormal(0.5 * ones(D))
    likelihood = MvNormal(7.0 * ones(D))
    logf = x -> logpdf(likelihood, x)
    m = BayesModel(p0, logf)

    x1 = randn(rng, Float64, D)
    @test integrand(m)(x1) ≈ pdf(likelihood, x1)
    @test logprior(m)(x1) ≈ logpdf(p0, x1)
    @test logjoint(m)(x1) ≈ logpdf(p0, x1) + logpdf(likelihood, x1)
    @test BQ.p_0(m) == p0
    @test logintegrand(m) == logf
end


text:
## TODO: Reorganize this. Migrate most of the Heads handled here elsewhere.
## Which heads should be here ?

### Apply

@mkapprule Apply nodefault => true """
    Apply(f,expr)   f .% expr

replace the `Head` of `expr` with `f`.

`Apply` also works with some
Julia objects. For example `Apply(Plus, :( [1:10] ))` returns `55`. `Apply` can be used
in operator form. For example

```
m = Apply(Plus)
m(f(a,b,c))
```
"""

@curry_first Apply

@doap Apply(f,g) = mx

@doap  function Apply(head::SJSym,mxa::Mxpr)
    head == :Plus && return _apply_plus(mxa)
    if head == :Times # 4 or 5 times faster for plus on numbers, don't evaluate
#        mx = mxpr(head,copy(margs(mxa))) # we may find that we need to copy
        mx = mxpra(head,margs(mxa))
        mx = canonexpr!(mx)            # this is ok
        setcanon(mx)
    else
        mx = mxpra(head,margs(mxa))
    end
    isa(mx,Mxpr) && isempty(mx) && return 0   # do this instead. fixes bug Apply(Times, [DirectedInfinity(),0]) --> 0
    mx
end

function _apply_plus(mxa::Mxpr)
    mx = mxpra(:Plus, margs(mxa))
    mx = canonexpr!(mx)
    setcanon(mx)
    return isa(mx, Mxpr) && isempty(mx) ? 0 : mx
end

@doap Apply(h,mxa::Mxpr) = mxpra(h,margs(mxa))

@doap function Apply(h::SJSym,arr::AbstractArray{T}) where T<:Number
    h == :Plus && return convert(T,sum(arr))
    h == :Times && return convert(T,prod(arr))
    mx
end

### Hash

@mkapprule Hash nargs => 1
@doap Hash(x) = hash(x)

### Head

@sjdoc Head """
    Head(expr)

return the `Head` of `expr`.

`expr` may be a Symata expression or object of any Julia type.

!!! note
    The `Head` of a Julia expression is `Expr`, for instance,
    `Head( :( :( a = 1) ))` returns `Expr`. Note we have to quote twice, because one level of
    a quoted Julia expression is evaluated so that we can embed Julia code.
"""
@mkapprule Head  nargs =>  1
@doap Head(x) = head(x)
head(x::Mxpr) = mhead(x)
head(x) = typeof(x)
# FIXME consider making :Symbol evaluate to Symbol, as, for instance, Int64
head(s::SJSym) = :Symbol

## `Function` is the head of Symata pure functions. So we use a different head for Julia functions
# Unclear what to do here.
# typeof( (x) -> x) in
# Julia v0.4 : Function
# Julia v0.5 and later : a hash key or something

@sjdoc CompiledFunction """
    CompiledFunction

is the head of compiled functions.

Compiled functions can be written directly in the host language, Julia.
```
f = J( x -> x^2 )
```

They may also be compiled from Symata expressions

```
f = Compile([x], x^2)
```
"""
@doap Head(f::T) where {T<:Function} = :CompiledFunction

### Isa

@sjdoc Isa """
    Isa(x,type)

return `True` if `x` is of type `type`.
"""
@mkapprule Isa  nargs => 2

@doap Isa(x,T::DataType) = isa(x,T)
@doap Isa(x,T::Symbol) = isa(x,eval(T))
@doap Isa(x,T) = mx

#const Float = AbstractFloat

### ReleaseHold

#typealias Holds Union{Mxpr{:Hold}, Mxpr{:HoldForm}, Mxpr{:HoldPattern}, Mxpr{:HoldComplete}}

@mkapprule ReleaseHold nargs => 1

@sjdoc ReleaseHold """
    ReleaseHold(expr)

removes the outer layer of `Hold`, `HoldForm`, `HoldPattern`, and `HoldComplete` from `expr`.
"""

@doap function ReleaseHold(mxa::Holds)
#    length(margs(mxa)) == 0 && return mxpr(:Sequence)  ## delete theses lines after a while
    isempty(mxa) && return mxpr(:Sequence)
    length(mxa) > 1 && return  mxpra(:Sequence,margs(mxa))
#    length(margs(mxa)) > 1 && return  mxpra(:Sequence,margs(mxa))
    return mxa[1]
end

@doap ReleaseHold(ex) = ex

### Reverse

function Base.reverse(mx::Mxpr)
    mx1 = copy(mx)
    Base.reverse!(margs(mx1))
    return mx1
end

@sjdoc Reverse """
    Reverse(expr)

reverse the order of the arguments in `expr`.
"""

@mkapprule Reverse nargs => 1
@doap function Reverse(ex::Mxpr)
    isOrderless(ex) && return ex
    setfixed(mxpra(mhead(ex),reverse(margs(ex))))
end
@doap Reverse(ex::AbstractArray) = reverse(ex)

@mkapprule Reverse! nargs => 1
@doap function Reverse!(ex::Mxpr)
    isOrderless(ex) && return ex
    reverse!(margs(ex))
    ex
end
@doap Reverse!(ex::AbstractArray) = reverse!(ex)

### Map

@mkapprule Map nargs => 1:3  """
    Map(f,expr)   f % expr

return `f` applied to each element in a `expr`.

    Map(f,expr,levelspec)

map at levels specified by `levelspec`.

`levelspec` is a standard level specification.

- `n`       levels `0` through `n`.
- `[n]`     level `n` only.
- `[n1,n2]` levels `n1` through `n2`

Negative indices count backwards from the deepest level.
`Infinity` specifies the deepest level.

`expr` is copied first. `f` can be a Symata object or a Julia function.

`Map` can be used in an operator form. For example `Map(f)(expr)`.
"""

@doap function Map(f::Function,expr::Mxpr)
    args = margs(expr)
    nargs = newargs(args)
    @inbounds for i in 1:length(args)
        nargs[i] = Base.invokelatest(f,args[i]) # Probably need more evaluation
    end
    mxpra(mhead(expr),nargs)
end

# Should we return an Array or a List ? We choose List now.
@doap function Map(f::Function, a::AbstractArray)
    nargs = newargs(length(a))
    @inbounds for i in 1:length(a)
        nargs[i] = Base.invokelatest(f,a[i])
    end
    mxpra(:List,nargs)
end

# We create one Mxpr outside the loop. Old
# code (commented out) created Mxpr every time.
# This saves 30 percent of time and allocation in some tests.
@doap Map(f,expr::Mxpr) =  _Map_one(f,expr::Mxpr)

function _Map_one(f,expr::Mxpr)
    args = margs(expr)
    nargs = newargs(args)
    mx1 = mxpr(f,0) # reserve one argument
    @inbounds for i in 1:length(args)
        mx1.args[1] = args[i]  # map f of one argument
        nargs[i] = doeval(mx1)
    end
    mxpra(mhead(expr),nargs)
end

## Optimize for a very special case.
## There are many such optimizations to do.
@doap function Map(sym::Symbol, expr::PlusT)
    sym != :Length && return _Map_one(sym,expr)
    sum(x -> symlength(x), expr)
end

mutable struct MapData
    action
end

## A bug in meaning of level spec of `2` is not caught because we check if p is Null below.
## level spec 2 should mean  from 1 to 2, but the level spec code interpreted it as 0 through 2.
@doap function Map(f,expr::Mxpr, inlevspec)
    levelspec = make_level_specification(expr, inlevspec)
    ex = recursive_copy(expr)
    data = MapData(nothing)
    action = LevelAction(data,
                         function (data,expr)
                         p = data.action.parent
                           if p != Null
                             p[data.action.subind] = mxpr(f,expr)
                           end
                         end)
    data.action = action
    traverse_levels!(action, levelspec, ex)
    if has_level_zero(levelspec)
        ex = mxpr(f,ex)
    end
    ex
end

@curry_first Map

### ToExpression

@sjdoc ToExpression """
    ToExpression(str)

convert string `str` to an expression.
"""

@mkapprule ToExpression nargs => 1:3
@doap ToExpression(s::String) = symataevaluate(Meta.parse(s), EvaluateJuliaSyntaxSimple())
@doap ToExpression(x) = x
set_sysattributes("ToExpression")

### Count

@sjdoc Count """
    Count(expr,pattern)

return the number of arguments in `expr` than match `pattern`.

Only matching on one level is supported. `Count(pattern)` can be used as the head of an expression,
as an operator. For instance, `cop = Count(_^2)` defines a function that counts
the number of arguments that have the form of a square.

`Count` also works when `expr` is a Julia `Dict`.
"""

@sjexamp( Count,
         ("Count(Range(10), _Integer)", "10"),
         ("Count(_Integer)(Range(10))", "10"),
         ("Count(Range(10), 2)", "1"))

set_sysattributes(:Count)

@mkapprule Count nargs => 1:3

@doap function Count(expr,pat)
    args = margs(expr)
    jp = patterntoBlank(pat)
    m = Match()  # allocating Match here is not faster than allocating capt here.
    count( x -> match_no_capture(x,jp,m), args)
end

@doap function Count(arr::AbstractArray, pat)
    jp = patterntoBlank(pat)
    m = Match()  # allocating Match here is not faster than allocating capt here.
    count(x -> match_no_capture(x, jp, m), arr)
end

## julia count is much faster than an explicit loop
@doap Count(expr::Mxpr,pat::Union{Number,Symbol,String}) = count(x -> x == pat, margs(expr))

@curry_last Count

### Cases

@sjdoc Cases """
    Cases(expr,pattern)

return the elements in `expr` that match `pattern`.

    Cases(expr,pattern,levelspec)

return the elements in `expr` on levels specified by `levelspec` that match `pattern`.

The head of the returned object is the same as that of `expr`.

`Cases(pattern)` can be used as the `Head` of an expression, as an operator.
For example, `getints = Cases(_Integer)`.
"""

@sjexamp( Cases,
         ("Cases([1,2.0,3,\"dog\"], _Integer)", "[1,3]"))

@mkapprule Cases nargs => 1:4

@doap function Cases(expr,pat)
    jp = patterntoBlank(pat)
    m = Match()
    MListA(filter(x -> match_no_capture(x,jp,m), margs(expr)))
end

@doap function Cases(expr,pat::RuleT)
    jp = patterntoBlank(pat)
    a = margs(expr)
    r = Array{eltype(a)}(undef, 0) ## same as newargs
    capt = capturealloc()
    for ai in a
        (gotmatch,res) = replace_ptob(ai,jp,capt)
        if gotmatch
            push!(r, res)
        end
    end
    return MListA(r)
end

## maybe better that pat be anything but Mxpr
@doap Cases(expr,pat::Union{Number,Symbol,String}) = mxpra(mhead(expr),filter(x -> x == pat, margs(expr)))

@doap function Cases(expr,pat,inlevelspec)
    levelspec = make_level_specification(expr, inlevelspec)
    _doCases(levelspec,expr,pat)
end

mutable struct CasesData
    new_args
    jp
    capt
end

function _doCases(levelspec::LevelSpec, expr ,pat)
    new_args = newargs()
    jp = patterntoBlank(pat)
    capt = capturealloc()
    data = CasesData(new_args,jp,capt)
    if isa(pat,RuleT)
        action = LevelAction(data, function (data, expr)
                             (gotmatch,res) = replace(expr,data.jp)
                             gotmatch ? push!(data.new_args,res) : nothing
                             end)
    else
        action = LevelAction(data, function (data, expr)
                             (gotmatch,capt) = match_and_capt(expr,data.jp,data.capt)
                             gotmatch ? push!(data.new_args,sjcopy(expr)) : nothing
                             end)
    end
    traverse_levels!(action,levelspec,expr)
    mxpra(:List,new_args)
end

@curry_last Cases

### DeleteCases

@sjdoc DeleteCases """
    DeleteCases(expr,pattern)

deletes the elements in `expr` that match `pattern`.

The head of the returned object is the same as that of expr. Only matching on one level is supported.

`DelteCases(pattern)` can be used as the head of an expression, as an operator.
For example `noints = DeleteCases(_Integer)`.
"""

@sjexamp( DeleteCases,
         ("DeleteCases([1,2.0,3,\"dog\"], _Integer)", "[2.0,\"dog\"]"))

@mkapprule DeleteCases nargs => 1:4

@doap function DeleteCases(expr,pat)
    args = margs(expr)
    new_args = newargs()
    jp = patterntoBlank(pat)
    capt = capturealloc()
    @inbounds for i in 1:length(args)
        (gotmatch,capt) = match_and_capt(args[i],jp,capt)
        gotmatch ? nothing : push!(new_args,sjcopy(args[i])) # The difference from Cases
    end
    rmx = mxpra(mhead(expr),new_args)
    return rmx
end

@curry_last DeleteCases

### Push!

@sjdoc Push! """
    Push!(a,val)

pushes `val` onto the expression that `Symbol` `a` evaluates to.

!!!  note
   Typically, computer algebra systems follow the principle that expressions should be immutable.
   Rather than modify an expression, a copy is modified. `Push!` violates this principle. It could
   interact poorly with other parts of Symata.
"""

@sjexamp( Push!,
         ("ClearAll(a,b)",""),
         ("a = []",""),
         ("For(i=1, i < 1000, Increment(i), Push!(a,Symbol(\"b\$i\")))",""))

set_sysattributes(["Push!"],[:HoldFirst])
apprules(mx::Mxpr{:Push!}) = do_Push(mx,margs(mx)...)
do_Push(mx,args...) = mx
do_Push(mx,x::SJSym,val) = do_Push1(mx,symval(x),val)
do_Push1(mx,x,val) = mx
do_Push1(mx,x::Mxpr,val) = (push!(x.args,val); x)

### Pop!

@sjdoc Pop! """
    Pop!(expr)

pops a value from the arguments of `expr`.

`Pop!` mutates `expr`.
"""

@mkapprule Pop!

@doap function Pop!(mx1::Mxpr)
    length(mx1) < 1 && return mx1 # and warn or error
    pop!(mx1.args)
end

@doap Pop!(x) = x

@sjseealso_group(Pop!, Push!)

### ComposeList

@sjdoc ComposeList """
    ComposeList([f1,f2,...],x)

returns `[f1(x),f2(f1(x)),...]`.
"""
@mkapprule ComposeList nargs => 2

@doap function ComposeList(list::ListT,x)
    ops = reverse(margs(list))
    mout = mxpr(ops[1],x)
    nargs = newargs(1)
    nargs[1] = mout
    for i in 2:length(ops)
        mout = mxpr(ops[i],mout)
        push!(nargs,mout)
    end
    mxpra(:List,nargs)
end

evalifdelayed(r::Mxpr{:Rule}) = rhs(r)
evalifdelayed(r::Mxpr{:RuleDelayed}) = doeval(rhs(r))

### ReplacePart

@mkapprule ReplacePart  nargs => 1:2  """
    ReplacePart(expr, i => repl)

returns a copy of `expr` with the `i`th part replaced by `repl`
"""

@doap ReplacePart(expr::Mxpr,arg::T) where {T<:Union{RulesT,ListT}} = replacepart1(mx,expr,arg)

@curry_second ReplacePart

function replacepart1(mx,expr,arg)
    nexpr = deepunsetfixed(recursive_copy(expr)) ## Mma must use some kind of lazy copying
    replacepart(mx,nexpr,arg)
end

function replacepart(mx,expr,rule::RulesT)
    _rhs = evalifdelayed(rule)
    if listofpredq(lhs(rule), integerq)
       expr = replaceonepart(mx,expr,lhs(rule),evalifdelayed(rule))
    elseif listq(lhs(rule)) # lhs is a List
        foreach( x -> (expr = replaceonepart(mx,expr,x,_rhs)), lhs(rule))
    else
       expr = replaceonepart(mx,expr,lhs(rule),_rhs)
    end
    expr
end

function replacepart(mx,expr,arg::ListT)
    if listofpredq(arg,ruleq)
        foreach( r -> (expr = replaceonepart(mx,expr,lhs(r),evalifdelayed(r))), arg)
    end
    expr
end

function replaceonepart(mx,expr,_lhs::Integer,_rhs)
    if _lhs == 0
        expr = mxprnewhead(expr,_rhs)
    else
        expr[posnegi(expr,_lhs)] = _rhs
    end
    expr
end

function replaceonepart(mx,expr,_lhs::ListT,_rhs)
    if ! listofpredq(_lhs,integerq)
        @symwarn("$_lhs is not a part specification")
        return mx
    end
    setpart2!(expr, _rhs,margs(_lhs)...)
    expr
end

### Level

@mkapprule Level nargs => 2:3

@sjdoc Level """
    Level(expr,levelspec)

returns a list of all parts at `levelspec`.

    Level(expr,levelspec,f)

applies `f` to each part in the returned list.

`Level` traverses the expression breadth first. Negative indices count from the
depth of the entire expression, rather than from each leaf.
"""

mutable struct LevelData
    levellist
    action
end

@doap function Level(expr, inlevelspec)
    levelspec = make_level_specification(expr,inlevelspec)
    nargs = newargs()
    data = LevelData(nargs,nothing)
    action = LevelAction(data,
                         function(data,expr)
#                         act = data.action
                           push!(data.levellist,expr)
                         end)
    data.action = action
    traverse_levels!(action, levelspec,expr)
    ## todo level zero
    mxpra(:List,nargs)
end

@doap function Level(expr, inlevelspec, f)
    levelspec = make_level_specification(expr,inlevelspec)
    nargs = newargs()
    data = LevelData(nargs,nothing)
    action = LevelAction(data,
                         function(data,expr)
                           push!(data.levellist,mxpr(f,expr))
                         end)
    data.action = action
    traverse_levels!(action, levelspec,expr)
    ## todo level zero
    mxpra(:List,nargs)
end

### ExpandA, only a bit is implemented. Sympy Expand is more capable.

@sjdoc ExpandA """
    ExpandA(expr)

expand products in `expr`. This is only partially implemented,
mostly to test the efficiency of evaluation and evaluation control. The
Sympy version `Expand()` is more capable, but slower.
"""

apprules(mx::Mxpr{:ExpandA}) = _doexpand(mx[1])


text:
@testset "Verify vander for ProHermiteBasis: evaluation, first and second derivatives" begin
    Ne = 200
    x = randn(Ne)
    #Test all the basis
    B = ProHermiteBasis(6)
    ktab = [0; 1; 2]
    for k in ktab
        # Test evaluation
        dV = vander(B, k, x)
        dVtrue = zeros(Ne, size(B))

        for i=1:B.m
            dVtrue[:,i] .= derivative(FamilyScaledProHermite[i], k , x)
        end
        @test norm(dV - dVtrue)<1e-6
    end
end

@testset "Verify vander for PhyHermiteBasis: evaluation, first and second derivatives" begin
    Ne = 200
    x = randn(Ne)
    #Test all the basis
    B = PhyHermiteBasis(6)
    ktab = [0; 1; 2]
    for k in ktab
        # Test evaluation
        dV = vander(B, k, x)
        dVtrue = zeros(Ne, size(B))

        for i=1:B.m
            dVtrue[:,i] .= derivative(FamilyScaledPhyHermite[i], k , x)
        end
        @test norm(dV - dVtrue)<1e-6
    end
end

@testset "Verify vander for CstProHermiteBasis: evaluation, first and second derivatives" begin
    Ne = 200
    x = randn(Ne)
    #Test all the basis
    B = CstProHermiteBasis(6)
    ktab = [0; 1; 2]
    for k in ktab
        # Test evaluation
        dV = vander(B, k, x)
        dVtrue = zeros(Ne, size(B))

        for i=1:B.m
            if i==1
                Pik = derivative(FamilyProPolyHermite[1], k)
                dVtrue[:,i] .= Pik.(x)
            else
                dVtrue[:,i] .= derivative(FamilyScaledProHermite[i-1], k , x)
            end
        end
        @test norm(dV - dVtrue)<1e-6
    end
end

@testset "Verify vander for CstPhyHermiteBasis: evaluation, first and second derivatives" begin
    Ne = 200
    x = randn(Ne)
    #Test all the basis
    B = CstPhyHermiteBasis(6)
    ktab = [0; 1; 2]

    for k in ktab
        # Test evaluation
        dV = vander(B, k, x)
        dVtrue = zeros(Ne, size(B))

        for i=1:B.m
            if i==1
                Pik = derivative(FamilyProPolyHermite[1], k)
                dVtrue[:,i] .= Pik.(x)
            else
                dVtrue[:,i] .= derivative(FamilyScaledPhyHermite[i-1], k , x)
            end
        end
        @test norm(dV - dVtrue)<1e-6
    end
end

@testset "Verify vander for CstLinProHermiteBasis: evaluation, first and second derivatives" begin
    Ne = 200
    x = randn(Ne)
    #Test all the basis
    B = CstLinProHermiteBasis(6)
    ktab = [0; 1; 2]

    for k in ktab
        # Test evaluation
        dV = vander(B, k, x)
        dVtrue = zeros(Ne, size(B))

        for i=1:B.m
            if i==1
                Pik = derivative(FamilyProPolyHermite[1], k)
                dVtrue[:,i] .= Pik.(x)
            elseif i==2
                Pik = derivative(FamilyProPolyHermite[2], k)
                dVtrue[:,i] .= Pik.(x)
            else
                dVtrue[:,i] .= derivative(FamilyScaledProHermite[i-2], k , x)
            end
        end
        @test norm(dV - dVtrue)<1e-6
    end
end

@testset "Verify vander for CstLinPhyHermiteBasis: evaluation, first and second derivatives" begin
    Ne = 200
    x = randn(Ne)
    #Test all the basis
    B = CstLinPhyHermiteBasis(6)
    ktab = [0; 1; 2]

    for k in ktab
        # Test evaluation
        dV = vander(B, k, x)
        dVtrue = zeros(Ne, size(B))

        for i=1:B.m
            if i==1
                Pik = derivative(FamilyProPolyHermite[1], k)
                dVtrue[:,i] .= Pik.(x)
            elseif i==2
                Pik = derivative(FamilyProPolyHermite[2], k)
                dVtrue[:,i] .= Pik.(x)
            else
                dVtrue[:,i] .= derivative(FamilyScaledPhyHermite[i-2], k , x)
            end
        end
        @test norm(dV - dVtrue)<1e-6
    end
end
#
#
# @testset "Verify vander for basis: evaluation, first and second derivatives" begin
#     Ne = 200
#     x = randn(Ne)
#     #Test all the basis
#     B = Basis(5)#CstProHermiteBasis(5; scaled = true)
#     # Btab = [CstProHermiteBasis(5; scaled = false);
#     #         CstProHermiteBasis(5; scaled = true)];
#             # CstPhyHermiteBasis(5; scaled = false);
#             # CstPhyHermiteBasis(5; scaled = true);
#             # CstLinPhyHermiteBasis(5; scaled = false);
#             # CstLinPhyHermiteBasis(5; scaled = true);
#             # CstLinProHermiteBasis(5; scaled = false);
#             # CstLinProHermiteBasis(5; scaled = true)]
#             # Test evaluation, derivative and second derivative
#     ktab = [0; 1; 2]
#
#     for k in ktab
#         # Test evaluation
#         dV = vander(B, k, x)
#         dVtrue = zeros(Ne, size(B))
#
#         for i=1:B.m
#             if i==1
#                 Pik = derivative(FamilyProPolyHermite[1], k)
#                 dVtrue[:,i] .= Pik.(x)
#             else
#                 dVtrue[:,i] .= derivative(FamilyScaledProHermite[i-1], k , x)
#             end
#         end
#         @test norm(dV - dVtrue)<1e-6
#     end
# end

# @testset "Test defined families list" begin
#
#     B = CstPhyHermiteBasis(10; scaled =false)
#     @test size(B,1)==12
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#
#     for i=2:size(B,1)
#         @test B[i] == FamilyPhyHermite[i-1]
#     end
#
#
#     B = CstPhyHermiteBasis(10; scaled =true)
#     @test size(B,1)==12
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#
#     for i=2:size(B,1)
#         @test B[i] == FamilyScaledPhyHermite[i-1]
#     end
#
#
#     B = CstProHermiteBasis(10; scaled =false)
#     @test size(B,1)==12
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#
#     for i=2:size(B,1)
#         @test B[i] == FamilyProHermite[i-1]
#     end
#
#
#     B = CstProHermiteBasis(10; scaled =true)
#     @test size(B,1)==12
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#
#     for i=2:size(B,1)
#         @test B[i] == FamilyScaledProHermite[i-1]
#     end
#
#     B = CstLinPhyHermiteBasis(10; scaled =false)
#     @test size(B,1)==13
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#     @test B[2] ==B.f[2] == FamilyProPolyHermite[2]
#
#
#     for i=3:size(B,1)
#         @test B[i] == FamilyPhyHermite[i-2]
#     end
#
#     B = CstLinPhyHermiteBasis(10; scaled = true)
#     @test size(B,1)==13
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#     @test B[2] ==B.f[2] == FamilyProPolyHermite[2]
#
#
#     for i=3:size(B,1)
#         @test B[i] == FamilyScaledPhyHermite[i-2]
#     end
#
#     B = CstLinProHermiteBasis(10; scaled =false)
#     @test size(B,1)==13
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#     @test B[2] ==B.f[2] == FamilyProPolyHermite[2]
#
#
#     for i=3:size(B,1)
#         @test B[i] == FamilyProHermite[i-2]
#     end
#
#     B = CstLinProHermiteBasis(10; scaled = true)
#     @test size(B,1)==13
#
#     @test B[1] ==B.f[1] == FamilyProPolyHermite[1]
#     @test B[2] ==B.f[2] == FamilyProPolyHermite[2]
#
#
#     for i=3:size(B,1)
#         @test B[i] == FamilyScaledProHermite[i-2]
#     end
# end

# @testset "Test evaluation of the basis" begin
#     # @testset "Test evaluation of the basis" begin
#         Ne = 10
#         x = 0.5*randn(Ne)
#
#         B = CstPhyHermiteBasis(5; scaled =false)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         for i=0:5
#             @test abs(Bx[i+2] -  PhyHermite(i)(x[1]))< 1e-8
#         end
#
#
#         B = CstPhyHermiteBasis(5; scaled = true)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         for i=0:5
#             @test abs(Bx[i+2] -  PhyHermite(i; scaled = true)(x[1]))< 1e-8
#         end
#
#         B = CstProHermiteBasis(5; scaled =false)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         for i=0:5
#             @test abs(Bx[i+2] -  ProHermite(i)(x[1]))< 1e-8
#         end
#
#
#         B = CstProHermiteBasis(5; scaled = true)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         for i=0:5
#             @test abs(Bx[i+2] -  ProHermite(i; scaled = true)(x[1]))< 1e-8
#         end
#
#         B = CstLinPhyHermiteBasis(5; scaled =false)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         @test Bx[2] == x[1]
#
#         for i=0:5
#             @test abs(Bx[i+3] -  PhyHermite(i)(x[1]))< 1e-8
#         end
#
#         B = CstLinPhyHermiteBasis(5; scaled =true)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         @test Bx[2] == x[1]
#
#         for i=0:5
#             @test abs(Bx[i+3] -  PhyHermite(i; scaled = true)(x[1]))< 1e-8
#         end
#
#         B = CstLinProHermiteBasis(5; scaled =false)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         @test Bx[2] == x[1]
#
#         for i=0:5
#             @test abs(Bx[i+3] -  ProHermite(i)(x[1]))< 1e-8
#         end
#
#         B = CstLinProHermiteBasis(5; scaled =true)
#         m = size(B,1)
#
#         Bx = B(x[1])
#         @test Bx[1] == 1.0
#         @test Bx[2] == x[1]
#
#         for i=0:5
#             @test abs(Bx[i+3] -  ProHermite(i; scaled = true)(x[1]))< 1e-8
#         end
# end


text:
using LightFSM
using Base.Test

@enum SINGLESTATE ss_state1
@enum SINGLEEVENT stay=1

@enum TWOSTATE ts_state1 ts_state2
@enum TWOSTATEEVENTS startup=1 st1_1 st1_2 st2_1 st2_2

@enum THREESTATE rs_state1 rs_state2 rs_state3
@enum THREESTATEEVENTS rs_st1_1=1 rs_st1_2 rs_st2_1 rs_st13_2 rs_st2_2 rs_st1_3 rs_st2_3 rs_st3_3 rs_st123_3

include("Basics.jl")
include("Multistate.jl")


text:
"""
    $SIGNATURES

Informative error message for missing method.
"""
@inline function no_model_method(f, args...)
    info("You need to define `$(string(f))` with this model type.")
    throw(MethodError(f, args))
end


text:
const FuncArray{T,N,B<:BasisOrRestricted} = MulQuasiArray{T,N,<:Tuple{B,<:AbstractArray{T,N}}}
const FuncVector{T,B<:BasisOrRestricted} = FuncArray{T,1,B}
const FuncMatrix{T,B<:BasisOrRestricted} = FuncArray{T,2,B}
const FuncVecOrMat{T,B<:BasisOrRestricted} = Union{FuncVector{T,B},FuncMatrix{T,B}}

const AdjointFuncArray{T,N,B<:BasisOrRestricted} = MulQuasiArray{T,<:Any,<:Tuple{<:Adjoint{T,<:AbstractArray{T,N}},
                                                                    <:QuasiAdjoint{T,<:B}}}
const AdjointFuncVector{T,B<:BasisOrRestricted} = AdjointFuncArray{T,1,B}
const AdjointFuncMatrix{T,B<:BasisOrRestricted} = AdjointFuncArray{T,2,B}
const AdjointFuncVecOrMat{T,B<:BasisOrRestricted} = Union{AdjointFuncVector{T,B},AdjointFuncMatrix{T,B}}


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct OfferTermInfo <: SwaggerModel
    Name::Any # spec type: Union{ Nothing, String } # spec name: Name
    EffectiveDate::Any # spec type: Union{ Nothing, DateTime } # spec name: EffectiveDate

    function OfferTermInfo(;Name=nothing, EffectiveDate=nothing)
        o = new()
        validate_property(OfferTermInfo, Symbol("Name"), Name)
        setfield!(o, Symbol("Name"), Name)
        validate_property(OfferTermInfo, Symbol("EffectiveDate"), EffectiveDate)
        setfield!(o, Symbol("EffectiveDate"), EffectiveDate)
        o
    end
end # type OfferTermInfo

const _property_map_OfferTermInfo = Dict{Symbol,Symbol}(Symbol("Name")=>Symbol("Name"), Symbol("EffectiveDate")=>Symbol("EffectiveDate"))
const _property_types_OfferTermInfo = Dict{Symbol,String}(Symbol("Name")=>"String", Symbol("EffectiveDate")=>"DateTime")
Base.propertynames(::Type{ OfferTermInfo }) = collect(keys(_property_map_OfferTermInfo))
Swagger.property_type(::Type{ OfferTermInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OfferTermInfo[name]))}
Swagger.field_name(::Type{ OfferTermInfo }, property_name::Symbol) =  _property_map_OfferTermInfo[property_name]

const _allowed_OfferTermInfo_Name = ["Recurring Charge", "Monetary Commitment", "Monetary Credit"]

function check_required(o::OfferTermInfo)
    true
end

function validate_property(::Type{ OfferTermInfo }, name::Symbol, val)
    if name === Symbol("Name")
        Swagger.validate_param(name, "OfferTermInfo", :enum, val, _allowed_OfferTermInfo_Name)
    end
end


text:
struct LSSymbolServerFailure <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSSymbolServerFailure)
    print(io, ex.msg)
end


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# int kvsprintf(kstring_t *s, const char *fmt, va_list ap) KS_ATTR_PRINTF(2,0);
@defun kvsprintf(s::Ptr{kstring_t}, fmt::Cstring, ap::va_list)::Cint

# int ksprintf(kstring_t *s, const char *fmt, ...) KS_ATTR_PRINTF(2,3);
# ErrorException("unknown token: '.'")

# int kputd(double d, kstring_t *s);
@defun kputd(d::Cdouble, s::Ptr{kstring_t})::Cint

# int ksplit_core(char *s, int delimiter, int *_max, int **_offsets);
@defun ksplit_core(s::Ptr{Cchar}, delimiter::Cint, _max::Ptr{Cint}, _offsets::Ptr{Ptr{Cint}})::Cint

# char *kstrstr(const char *str, const char *pat, int **_prep);
@defun kstrstr(str::Cstring, pat::Cstring, _prep::Ptr{Ptr{Cint}})::Ptr{Cchar}

# char *kstrnstr(const char *str, const char *pat, int n, int **_prep);
@defun kstrnstr(str::Cstring, pat::Cstring, n::Cint, _prep::Ptr{Ptr{Cint}})::Ptr{Cchar}

# void *kmemmem(const void *_str, int n, const void *_pat, int m, int **_prep);
@defun kmemmem(_str::Ptr{Cvoid}, n::Cint, _pat::Ptr{Cvoid}, m::Cint, _prep::Ptr{Ptr{Cint}})::Ptr{Cvoid}

# char *kstrtok(const char *str, const char *sep, ks_tokaux_t *aux);
@defun kstrtok(str::Cstring, sep::Cstring, aux::Ptr{ks_tokaux_t})::Ptr{Cchar}

# int kgetline(kstring_t *s, kgets_func *fgets, void *fp);
@defun kgetline(s::Ptr{kstring_t}, fgets::Ptr{kgets_func}, fp::Ptr{Cvoid})::Cint

# int kgetline2(kstring_t *s, kgets_func2 *fgets, void *fp);
@defun kgetline2(s::Ptr{kstring_t}, fgets::Ptr{kgets_func2}, fp::Ptr{Cvoid})::Cint


text:
include("plume_rise.jl")

"""
    gaussian_plume_factory(scenario; downwash=false, plumerise=false)

Generates a gaussian dispersion model on the given scenario and returns a
function giving the concentration of the form
c(x, y, z[, t])

`downwash` controls whether or not stack-downwash is included, by default it
is not
`plumerise` controls whether or not plume height is adjusted, by default there
is no plume rise
"""
function gaussian_plume_factory(scenario; downwash=false, plumerise=false)

    required_params = [:mass_emission_rate, :jet_diameter, :jet_velocity,
                       :release_height, :windspeed, :pasquill_gifford]
    if all(key -> !(ismissing(getproperty(scenario,key))), required_params)
        # parameters of the jet
        Q = scenario.mass_emission_rate
        Dⱼ = scenario.jet_diameter
        uⱼ = scenario.jet_velocity
        hᵣ = scenario.release_height

        # parameters of the environment
        u = scenario.windspeed
        stability = scenario.pasquill_gifford
    else
        missing_params = [ String(i) for i in filter(key -> ismissing(getproperty(scenario,key)), required_params)]
        error_string = "These parameters cannot be missing: " * join(missing_params, ", ")
        e = MissingException(error_string)
        throw(e)
    end

    # stack-tip downwash check
    if (downwash==true) && (uⱼ < 1.5*u)
        Δh_dw = 2*Dⱼ*( (uⱼ/u) - 1.5 )
    else
        Δh_dw = 0.0
    end

    hᵣ = hᵣ + Δh_dw

    # plume rise
    Δh = plume_rise(scenario, plumerise)

    # Pasquill-Gifford dispersion
    σy = crosswind_dispersion(stability)
    σz = vertical_dispersion(stability)

    function gaussian_plume(x, y, z, t=0)
        hₑ  = hᵣ + Δh(x)
        σyₑ = √( (Δh(x)/3.5)^2 + σy(x)^2 )
        σzₑ = √( (Δh(x)/3.5)^2 + σz(x)^2 )

        gaussian_plume = ( Q/(2*π*u*σyₑ*σzₑ)
                         * exp(-0.5*(y/σyₑ)^2)
                         *( exp(-0.5*((z-hₑ)/σzₑ)^2) + exp(-0.5*((z+hₑ)/σzₑ)^2) ))

    end

    return gaussian_plume

end


text:
_clamp(x) = Point2f0(clamp(x[1], 0, 1), clamp(x[2], 0, 1))
function widget{T<:RGBA}(
        color::Signal{T}, window;
        area=(30, 30),
        kw_args...
    )
    @materialize mouse_buttons_pressed, mouseposition = window.inputs
    color_button = visualize(
        (ROUNDED_RECTANGLE, zeros(Point2f0, 1));
        scale=const_lift(x->Vec2f0(x[2]), area),
        offset=Vec2f0(0),
        color=color,
        stroke_width=1f0,
        stroke_color=RGBA{Float32}(0.9, 0.9, 0.9, 1.0),
        kw_args...
    )
    color_robj = color_button.children[]
    m2id = GLWindow.mouse2id(window)
    isoverpoint = const_lift(is_same_id, m2id, color_robj)
    key_pressed = map(mouse_buttons_pressed) do mbuttons
        if length(mbuttons) == 1
            b = first(mbuttons)
            return GLFW.MOUSE_BUTTON_LEFT == b || b == GLFW.MOUSE_BUTTON_RIGHT
        end
        false
    end
    # dragg while key_pressed. Drag only starts if isoverpoint is true
    mousedragg = GLAbstraction.dragged(mouseposition, key_pressed, isoverpoint)

    preserve(foldp((value(m2id)..., value(color)), mousedragg) do v0, dragg
        id, index, p0 = v0
        if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
            id, index = value(m2id)
            if id==color_robj.id
                p0 = value(color)
            end
        else
            if id==color_robj.id
                mbutton = first(value(mouse_buttons_pressed))
                if mbutton == GLFW.MOUSE_BUTTON_RIGHT # with the right button we control red and green
                    r, g = _clamp(Point2f0(p0.r, p0.g) + Point2f0(dragg/50f0))
                    push!(color, RGBA{eltype(T)}(r, g, p0.b, p0.alpha))
                elseif mbutton == GLFW.MOUSE_BUTTON_LEFT # with the left button we control blue and alpha
                    b, a = _clamp(Point2f0(p0.b, p0.alpha) + Point2f0(dragg/50f0))
                    push!(color, RGBA{eltype(T)}(p0.r, p0.g, b, a))
                end
            end
        end
        return id, index, p0
    end)
    return color, color_button
end


text:
export kaplanyorke_dim

"""
    kaplanyorke_dim(λs::AbstractVector)
Calculate the Kaplan-Yorke dimension, a.k.a. Lyapunov dimension[^Kaplan1970].

## Description
The Kaplan-Yorke dimension is simply the point where
`cumsum(λs)` becomes zero (interpolated):
```math
 D_{KY} = k + \\frac{\\sum_{i=1}^k \\lambda_i}{|\\lambda_{k+1}|},\\quad k = \\max_j \\left[ \\sum_{i=1}^j \\lambda_i > 0 \\right].
```

If the sum of the exponents never becomes negative the function
will return the length of the input vector.

Useful in combination with [`lyapunovspectrum`](@ref).

[^Kaplan1970]: J. Kaplan & J. Yorke, *Chaotic behavior of multidimensional difference equations*, Lecture Notes in Mathematics vol. **730**, Springer (1979)
"""
function kaplanyorke_dim(v::AbstractVector)
    issorted(v, rev = true) || throw(ArgumentError(
    "The lyapunov vector must be sorted from most positive to most negative"))

    s = cumsum(v); k = length(v)
    # Find k such that sum(λ_i for i in 1:k) is still possitive
    for i in eachindex(s)
        if s[i] < 0
            k = i-1
            break
        end
    end

    if k == 0
        return zero(v[1])
    elseif k < length(v)
        return k + s[k]/abs(v[k+1])
    else
        return typeof(v[1])(length(v))
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkGatewayConnectionStatus <: SwaggerModel

    function VirtualNetworkGatewayConnectionStatus(;)
        o = new()
        o
    end
end # type VirtualNetworkGatewayConnectionStatus

const _property_map_VirtualNetworkGatewayConnectionStatus = Dict{Symbol,Symbol}()
const _property_types_VirtualNetworkGatewayConnectionStatus = Dict{Symbol,String}()
Base.propertynames(::Type{ VirtualNetworkGatewayConnectionStatus }) = collect(keys(_property_map_VirtualNetworkGatewayConnectionStatus))
Swagger.property_type(::Type{ VirtualNetworkGatewayConnectionStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkGatewayConnectionStatus[name]))}
Swagger.field_name(::Type{ VirtualNetworkGatewayConnectionStatus }, property_name::Symbol) =  _property_map_VirtualNetworkGatewayConnectionStatus[property_name]

function check_required(o::VirtualNetworkGatewayConnectionStatus)
    true
end

function validate_property(::Type{ VirtualNetworkGatewayConnectionStatus }, name::Symbol, val)
end


text:
using MarkovTransitionMatrices

export makeSparseTransition

function makeSparseTransition(nP::Integer=31, nC::Integer=37, extendSD::Real = 1.5, minP::Real = 1e-6)
  n = nP * nC

  T = Float64

  σp = [0x1.269c263357e05p-4, 0x1.9bec08eb57abfp-5]
  σc = [0x1.2c9a3a5b8b5bap-5, 0x1.09d5d411d8a33p-7]

  Πk = [0x1.d7e13e2f38db1p-1 0x1.40f60e8638b53p-4;
        0x1.098244df0fddfp-4 0x1.decfb7641df78p-1]
  all(sum(Πk,2) .≈ 1.0) || throw(error("each row of π must sum to 1"))

  K = size(Πk,2)
  extrema_p = [0.8776572, 2.485073]
  extrema_c = [0.2437302, 1.529529]

  # state-spaces
  pspace = range(extrema_p[1]-log(extendSD), stop=extrema_p[2] + log(extendSD), length=nP)
  cspace = range(extrema_c[1]-log(extendSD), stop=extrema_c[2] + log(extendSD), length=nC)
  zspace = (pspace, cspace, 1:K,)

  # make deviations
  zrandwalk(x::Real, st::Real, σ::Real) = (x - st) / σ

  # allocate giant transition matrix
  P = zeros(T,n*K,n*K)

  # match moments for logp
  Pp1, JN, Λ, L_p1, approxErr = discreteNormalApprox(pspace, pspace, (x::Real,st::Real) -> zrandwalk(x,st,σp[1]), 15)
  Pp2, JN, Λ, L_p2, approxErr = discreteNormalApprox(pspace, pspace, (x::Real,st::Real) -> zrandwalk(x,st,σp[2]), 15)
  # plot(pspace, [L_p1,L_p2], yticks = 1:15, labels=["Hi vol", "Lo vol"], xlabel="Log p", ylabel="Moments matched")

  # match moments for logc
  Pc1, JN, Λ, L_c1, approxErr = discreteNormalApprox(cspace, cspace, (x::Real,st::Real) -> zrandwalk(x,st,σc[1]), 13)
  Pc2, JN, Λ, L_c2, approxErr = discreteNormalApprox(cspace, cspace, (x::Real,st::Real) -> zrandwalk(x,st,σc[2]), 5)
  # plot(cspace, [L_c1,L_c2], yticks = 1:15, labels=["Hi vol", "Lo vol"], xlabel="Log c", ylabel="Moments matched")

  # sparse versions
  sPp1 = MarkovTransitionMatrices.sparsify!(Pp1, minp)
  sPp2 = MarkovTransitionMatrices.sparsify!(Pp2, minp)
  sPc1 = MarkovTransitionMatrices.sparsify!(Pc1, minp)
  sPc2 = MarkovTransitionMatrices.sparsify!(Pc2, minp)

  # do the big matrix
  P12 = (kron(sPc1, sPp1), kron(sPc2, sPp2))
  for r2 in 1:K
    for r1 in 1:K
      P[(r1-1)*n+1:r1*n,  (r2-1)*n+1:r2*n] .= Πk[r1, r2] .* P12[r2]
    end
  end

  sP = sparse(P)

  return (zspace, sP)
end


text:
# source: https://github.com/salvipeter/dual-contouring
#module DualContouring

#using LinearAlgebra

const ϵ = 1e-15

"""
    isosurface_point(points, df, limits)

Computes the likely position of a point on the isosurface, based on 0-crossing on the
cell's edges given in `points`. `df` is the gradient of the function,
and `limits` is a pair of points defining the size of the cell.

When `df` is `nothing`, the isosurface point is created based only on the distance values.
"""
function isosurface_point(points, df, limits)
    m = length(points)
    (df === nothing || m < 3) && return sum(points) / m
    A = Array{Float64}(undef, m, 3)
    b = Vector{Float64}(undef, m)
    for i in 1:m
        n = df(points[i])
        A[i,:] = n
        b[i] = dot(n, points[i])
    end
    q = A \ b
    any(x < 0 for x in q - limits[1]) || any(x < 0 for x in limits[2] - q) ? sum(points) / m : q
end

"""
    compute_cell(origin, dirs, vertices, df)

Computes an approximate position of the isosurface in a cell. `origin` is the bottom-left-front
corner of the cell, with `dirs` contains 3 axis-aligned vectors going to the opposite corner.
`vertices` gives the function values at the cell corners in this fashion:
```
      7-----8
     /|    /|     Y
    3-----4 |     ^
    | 5---|-6     | Z
    |/    |/      |/
    1-----2       O----> X
```
`df` is the gradient function (when available), or `nothing`.

It is assumed that the function values are never exactly zero.
When all values have the same sign, the function returns `nothing`.
"""
function compute_cell(origin, dirs, vertices, df)
    (all(v < 0 for v in vertices) || all(v > 0 for v in vertices)) && return nothing
    edges = [0, 4, 1, 5, 2, 6, 3, 7,
             0, 2, 1, 3, 4, 6, 5, 7,
             0, 1, 2, 3, 4, 5, 6, 7]
    mean = [0, 0, 0]
    points = []
    for i in 1:12
        v1, v2 = edges[2i-1], edges[2i]
        a, b = vertices[v1+1], vertices[v2+1]
        a * b > 0 && continue
        denom = abs(b - a)
        x = denom < ϵ ? 0.5 : abs(a) / denom
        c = (i - 1) ÷ 4 + 1
        p = origin + dirs[1] * (v1 ÷ 4) + dirs[2] * (v1 % 4 ÷ 2) + dirs[3] * (v1 % 2) + dirs[c] * x
        push!(points, p)
    end
    isosurface_point(points, df, (origin, origin + sum(dirs)))
end

"""
    generate_cells(values, df, corner, delta, resolution)

Generates cells with dual points for an isosurface. `values` is a 3D matrix of dimension
`resolution[i]+1` in each direction, containing the function values. `corner` is the
bottom-left-front corner of the volume, with the `delta` vector going to the opposite corner.
`resolution` is a tuple of 3 values, describing the number of cells in each axis direction.

The return value is a `(points, cells)` pair, where `points` is a list of approximated pionts
of the isosurface, and `cells` is a 3D matrix of dimension `resolution[i]` in each direction,
containing an index to the `points` array, or 0 when empty.
"""
function generate_cells(values, df, corner, delta, resolution)
    cells = zeros(Int, resolution)
    points = []
    dirs = [[delta[1], 0, 0],
            [0, delta[2], 0],
            [0, 0, delta[3]]]
    for i in 1:resolution[1], j in 1:resolution[2], k in 1:resolution[3]
        vertices = []
        for di in 0:1, dj in 0:1, dk in 0:1
            push!(vertices, values[i+di,j+dj,k+dk])
        end
        origin = corner + [delta[1] * (i - 1), delta[2] * (j - 1), delta[3] * (k - 1)]
        surface_point = compute_cell(origin, dirs, vertices, df)
        if surface_point != nothing
            push!(points, surface_point)
            cells[i,j,k] = length(points)
        end
    end

    points, cells
end

"""
    generate_faces(values, cells, resolution)

Generate faces for the cells created by `generate_cells`, when possible. `values` and
`resolution` are defined the same way as in that function.

The return value is an array of quads, where each quad is an array containing 4 indices.
"""
function generate_faces(values, cells, resolution)
    quads = []
    for (r, value, cell) in [([1, 2, 3], (i, j, k) -> values[i,j,k], (i, j, k) -> cells[i,j,k]),
                             ([2, 3, 1], (i, j, k) -> values[j,k,i], (i, j, k) -> cells[j,k,i]),
                             ([3, 1, 2], (i, j, k) -> values[k,i,j], (i, j, k) -> cells[k,i,j])]
        for i in 1:resolution[r[1]], j in 2:resolution[r[2]], k in 2:resolution[r[3]]
            a, b, c, d = cell(i,j,k), cell(i,j-1,k), cell(i,j-1,k-1), cell(i,j,k-1)
            a * b * c * d == 0 && continue
            v1, v2 = value(i,j,k), value(i+1,j,k)
            v1 * v2 > 0 && continue
            if v1 < 0
                push!(quads, [a, b, c, d])
            else
                push!(quads, [d, c, b, a])
            end
        end
    end
    quads
end

"""
    isosurface(f, isolevel, bounding_box, resolution; df)

Generates the isosurface of `f` at the given `isolevel`. The `bounding_box` is a tuple
containing the bottom-left-front and top-right-back points. The `resolution` is given as a
tuple of 3 integers - the number of cells in each of the axis directions.

When `df` is given, it is assumed to be the gradient function.

Sphere example:

    julia> using LinearAlgebra
    julia> isosurface(p -> norm(p) - 1, 0, ([-2, -2, -2], [2, 2, 2]), (20, 20, 20))
"""
function isosurface(f, isolevel, bounding_box, resolution; df=nothing)
    delta = (bounding_box[2] - bounding_box[1]) ./ resolution
    values = zeros(resolution .+ 1)
    for i in 0:resolution[1], j in 0:resolution[2], k in 0:resolution[3]
        v = f(bounding_box[1] + delta .* [i, j, k]) - isolevel
        values[i+1,j+1,k+1] = abs(v) < ϵ ? copysign(ϵ, v) : v
    end
    points, cells = generate_cells(values, df, bounding_box[1], delta, resolution)
    faces = generate_faces(values, cells, resolution)
    points, faces
end

"""
    writeOBJ(points, faces, filename)

Writes a quad mesh in OBJ format to the given file.
"""
function writeOBJ(points, faces, filename)
    open(filename, "w") do fp
        for p in points
            println(fp, "v $(p[1]) $(p[2]) $(p[3])")
        end
        for f in faces
            println(fp, "f $(f[1]) $(f[2]) $(f[3]) $(f[4])")
        end
    end
end

#end # module