text:
# Notational note
# JMW's dx_history <=> NW's S
# JMW's dg_history <=> NW's Y

# Here alpha is a cache that parallels betas
# It is not the step-size
# q is also a cache
function twoloop!(s,
                  gr,
                  rho,
                  dx_history,
                  dg_history,
                  m::Integer,
                  pseudo_iteration::Integer,
                  alpha,
                  q,
                  scaleinvH0::Bool,
                  precon)
    # Count number of parameters
    n = length(s)

    # Determine lower and upper bounds for loops
    lower = pseudo_iteration - m
    upper = pseudo_iteration - 1

    # Copy gr into q for backward pass
    copyto!(q, gr)

    # Backward pass
    for index in upper:-1:lower
        if index < 1
            continue
        end
        i   = mod1(index, m)
        dgi = dg_history[i]
        dxi = dx_history[i]
        @inbounds alpha[i] = rho[i] * real(dot(dxi, q))
        @inbounds q .-= alpha[i] .* dgi
    end

    # Copy q into s for forward pass
    # apply preconditioner if precon != nothing
    # (Note: preconditioner update was done outside of this function)
    if scaleinvH0 == true && pseudo_iteration > 1
        # Use the initial scaling guess from
        # Nocedal & Wright (2nd ed), Equation (7.20)

        #=
        pseudo_iteration > 1 prevents this scaling from happening
        at the first iteration, but also at the first step after
        a reset due to invH being non-positive definite (pseudo_iteration = 1).
        TODO: Maybe we can still use the scaling as long as iteration > 1?
        =#
        i = mod1(upper, m)
        dxi = dx_history[i]
        dgi = dg_history[i]
        scaling = real(dot(dxi, dgi)) / sum(abs2, dgi)
        @. s = scaling*q
    else
        ldiv!(s, precon, q)
    end
    # Forward pass
    for index in lower:1:upper
        if index < 1
            continue
        end
        i = mod1(index, m)
        dgi = dg_history[i]
        dxi = dx_history[i]
        @inbounds beta = rho[i] * real(dot(dgi, s))
        @inbounds s .+= dxi .* (alpha[i] - beta)
    end

    # Negate search direction
    rmul!(s, eltype(s)(-1))

    return
end

struct LBFGS{T, IL, L, Tprep<:Union{Function, Nothing}} <: FirstOrderOptimizer
    m::Int
    alphaguess!::IL
    linesearch!::L
    P::T
    precondprep!::Tprep
    manifold::Manifold
    scaleinvH0::Bool
end
"""
# LBFGS
## Constructor
```julia
LBFGS(; m::Integer = 10,
alphaguess = LineSearches.InitialStatic(),
linesearch = LineSearches.HagerZhang(),
P=nothing,
precondprep = (P, x) -> nothing,
manifold = Flat(),
scaleinvH0::Bool = true && (typeof(P) <: Nothing))
```
`LBFGS` has two special keywords; the memory length `m`,
and the `scaleinvH0` flag.
The memory length determines how many previous Hessian
approximations to store.
When `scaleinvH0 == true`,
then the initial guess in the two-loop recursion to approximate the
inverse Hessian is the scaled identity, as can be found in Nocedal and Wright (2nd edition) (sec. 7.2).

In addition, LBFGS supports preconditioning via the `P` and `precondprep`
keywords.

## Description
The `LBFGS` method implements the limited-memory BFGS algorithm as described in
Nocedal and Wright (sec. 7.2, 2006) and original paper by Liu & Nocedal (1989).
It is a quasi-Newton method that updates an approximation to the Hessian using
past approximations as well as the gradient.

## References
 - Wright, S. J. and J. Nocedal (2006), Numerical optimization, 2nd edition. Springer
 - Liu, D. C. and Nocedal, J. (1989). "On the Limited Memory Method for Large Scale Optimization". Mathematical Programming B. 45 (3): 503–528
"""
function LBFGS(; m::Integer = 10,
                 alphaguess = LineSearches.InitialStatic(), # TODO: benchmark defaults
                 linesearch = LineSearches.HagerZhang(),  # TODO: benchmark defaults
                 P=nothing,
                 precondprep = (P, x) -> nothing,
                 manifold::Manifold=Flat(),
                 scaleinvH0::Bool = true && (typeof(P) <: Nothing) )
    LBFGS(Int(m), alphaguess, linesearch, P, precondprep, manifold, scaleinvH0)
end

Base.summary(::LBFGS) = "L-BFGS"

mutable struct LBFGSState{Tx, Tdx, Tdg, T, G} <: AbstractOptimizerState
    x::Tx
    x_previous::Tx
    g_previous::G
    rho::Vector{T}
    dx_history::Tdx
    dg_history::Tdg
    dx::Tx
    dg::Tx
    u::Tx
    f_x_previous::T
    twoloop_q
    twoloop_alpha
    pseudo_iteration::Int
    s::Tx
    @add_linesearch_fields()
end

function initial_state(method::LBFGS, options, d, initial_x)
    T = real(eltype(initial_x))
    n = length(initial_x)
    initial_x = copy(initial_x)
    retract!(method.manifold, initial_x)

    value_gradient!!(d, initial_x)

    project_tangent!(method.manifold, gradient(d), initial_x)
    LBFGSState(initial_x, # Maintain current state in state.x
              copy(initial_x), # Maintain previous state in state.x_previous
              similar(gradient(d)), # Store previous gradient in state.g_previous
              Vector{T}(undef, method.m), # state.rho
              [similar(initial_x) for i = 1:method.m], # Store changes in position in state.dx_history
              [similar(gradient(d)) for i = 1:method.m], # Store changes in position in state.dg_history
              similar(initial_x), # Buffer for new entry in state.dx_history
              similar(initial_x), # Buffer for new entry in state.dg_history
              similar(initial_x), # Buffer stored in state.u
              real(T)(NaN), # Store previous f in state.f_x_previous
              similar(initial_x), #Buffer for use by twoloop
              Vector{T}(undef, method.m), #Buffer for use by twoloop
              0,
              similar(initial_x), # Store current search direction in state.s
              @initial_linesearch()...)
end

function update_state!(d, state::LBFGSState, method::LBFGS)
    n = length(state.x)
    # Increment the number of steps we've had to perform
    state.pseudo_iteration += 1

    project_tangent!(method.manifold, gradient(d), state.x)

    # update the preconditioner
    method.precondprep!(method.P, state.x)

    # Determine the L-BFGS search direction # FIXME just pass state and method?
    twoloop!(state.s, gradient(d), state.rho, state.dx_history, state.dg_history,
             method.m, state.pseudo_iteration,
             state.twoloop_alpha, state.twoloop_q, method.scaleinvH0, method.P)
    project_tangent!(method.manifold, state.s, state.x)

    # Save g value to prepare for update_g! call
    copyto!(state.g_previous, gradient(d))

    # Determine the distance of movement along the search line
    lssuccess = perform_linesearch!(state, method, ManifoldObjective(method.manifold, d))

    # Update current position
    state.dx .= state.alpha .* state.s
    state.x .= state.x .+ state.dx
    retract!(method.manifold, state.x)

    lssuccess == false # break on linesearch error
end

function update_h!(d, state, method::LBFGS)
    n = length(state.x)
    # Measure the change in the gradient
    state.dg .= gradient(d) .- state.g_previous

    # Update the L-BFGS history of positions and gradients
    rho_iteration = one(eltype(state.dx)) / real(dot(state.dx, state.dg))
    if isinf(rho_iteration)
        # TODO: Introduce a formal error? There was a warning here previously
        return true
    end
    idx = mod1(state.pseudo_iteration, method.m)
    @inbounds state.dx_history[idx] .= state.dx
    @inbounds state.dg_history[idx] .= state.dg
    @inbounds state.rho[idx] = rho_iteration
end

function trace!(tr, d, state, iteration, method::LBFGS, options, curr_time=time())
  common_trace!(tr, d, state, iteration, method, options, curr_time)
end


text:
module JSON

using Compat

export json # returns a compact (or indented) JSON representation as a String

include("Parser.jl")

import .Parser.parse

const INDENT=true
const NOINDENT=false
const unescaped = Bool[isprint(c) && !iscntrl(c) && !(c in ['\\','"']) for c in '\x00':'\x7F']

type State{I}
    indentstep::Int
    indentlen::Int
    prefix::AbstractString
    otype::Array{Bool, 1}
    State(indentstep::Int) = new(indentstep,
                                 0,
                                 "",
                                 Bool[])
end
State(indentstep::Int=0) = State{indentstep>0}(indentstep)

function set_state(state::State{INDENT}, operate::Int)
    state.indentlen += state.indentstep * operate
    state.prefix = " "^state.indentlen
end

set_state(state::State{NOINDENT}, operate::Int) = nothing

suffix(::State{INDENT})   = "\n"
suffix(::State{NOINDENT}) = ""

prefix(s::State{INDENT})  = s.prefix
prefix(::State{NOINDENT}) = ""

colon(s::State{INDENT})  = ": "
colon(::State{NOINDENT}) = ":"

# short hand for printing suffix then prefix
printsp(io::IO, state::State{INDENT}) = Base.print(io, suffix(state), prefix(state))
printsp(io::IO, state::State{NOINDENT}) = nothing

function start_object(io::IO, state::State{INDENT}, is_dict::Bool)
    push!(state.otype, is_dict)
    Base.print(io, is_dict ? "{": "[", suffix(state))
    set_state(state, 1)
end

function start_object(io::IO, state::State{NOINDENT}, is_dict::Bool)
    Base.print(io, is_dict ? "{": "[")
end

function end_object(io::IO, state::State{INDENT}, is_dict::Bool)
    set_state(state, -1)
    pop!(state.otype)
    printsp(io, state)
    Base.print(io, is_dict ? "}": "]")
end

function end_object(io::IO, state::State{NOINDENT}, is_dict::Bool)
    Base.print(io, is_dict ? "}": "]")
end

function print_escaped(io::IO, s::AbstractString)
    for c in s
        c <= '\x7f' ? (unescaped[@compat(Int(c))+1]  ? Base.print(io, c) :
                       c == '\\'       ? Base.print(io, "\\\\") :
                       c == '"'        ? Base.print(io, "\\\"") :
                       8 <= c <= 10    ? Base.print(io, '\\', "btn"[@compat(Int(c))-7]) :
                       c == '\f'       ? Base.print(io, "\\f") :
                       c == '\r'       ? Base.print(io, "\\r") :
                                         Base.print(io, "\\u", hex(c, 4))) :
                      Base.print(io, c) #JSON is UTF8 encoded
    end
end

function _print(io::IO, state::State, s::AbstractString)
    Base.print(io, '"')
    JSON.print_escaped(io, s)
    Base.print(io, '"')
end

function _print(io::IO, state::State, s::Symbol)
    _print(io, state, string(s))
end

function _print(io::IO, state::State, s::Union(Integer, FloatingPoint))
    if isnan(s) || isinf(s)
        Base.print(io, "null")
    else
        Base.print(io, s)
    end
end

function _print(io::IO, state::State, n::Nothing)
    Base.print(io, "null")
end

function _print(io::IO, state::State, a::Associative)
    if length(a) == 0
        Base.print(io, "{}")
        return
    end
    start_object(io, state, true)
    first = true
    for (key, value) in a
        first ? (first = false) : Base.print(io, ",", suffix(state))
        Base.print(io, prefix(state))
        JSON._print(io, state, string(key))
        Base.print(io, colon(state))
        JSON._print(io, state, value)
    end
    end_object(io, state, true)
end

function _print(io::IO, state::State, a::Union(AbstractVector,Tuple))
    if length(a) == 0
        Base.print(io, "[]")
        return
    end
    start_object(io, state, false)
    Base.print(io, prefix(state))
    i = start(a)
    !done(a,i) && ((x, i) = next(a, i); JSON._print(io, state, x); )

    while !done(a,i)
        (x, i) = next(a, i)
        Base.print(io, ",")
        printsp(io, state)
        JSON._print(io, state, x)
    end
    end_object(io, state, false)
end

function _print(io::IO, state::State, a)
    start_object(io, state, true)
    range = typeof(a).names
    if length(range) > 0
        Base.print(io, prefix(state), "\"", range[1], "\"", colon(state))
        JSON._print(io, state, a.(range[1]))

        for name in range[2:end]
            Base.print(io, ",")
            printsp(io, state)
            Base.print(io, "\"", name, "\"", colon(state))
            JSON._print(io, state, a.(name))
        end
    end
    end_object(io, state, true)
end

function _print(io::IO, state::State, f::Function)
    Base.print(io, "\"function at ", f.fptr, "\"")
end

function _print(io::IO, state::State, d::DataType)
    Base.print(io, "\"", d, "\"")
end

function _print(::IO, ::State, m::Module)
    throw(ArgumentError("cannot serialize Module $m as JSON"))
end

# Note: Arrays are printed in COLUMN MAJOR format.
# i.e. json([1 2 3; 4 5 6]) == "[[1,4],[2,5],[3,6]]"
function _print{T, N}(io::IO, state::State, a::AbstractArray{T, N})
    lengthN = size(a, N)
    if lengthN >= 0
        start_object(io, state, false)
        newdims = ntuple(N - 1, i -> 1:size(a, i))
        Base.print(io, prefix(state))
        JSON._print(io, state, slice(a, newdims..., 1))

        for j in 2:lengthN
            Base.print(io, ",")
            printsp(io, state)
            newdims = ntuple(N - 1, i -> 1:size(a, i))
            JSON._print(io, state, slice(a, newdims..., j))
        end
        end_object(io, state, false)
    else
        Base.print(io, "[]")
    end
end

function print(io::IO, a, indent=0)
    JSON._print(io, State(indent), a)
    if indent > 0
        Base.print(io, "\n")
    end
end

function print(a, indent=0)
    JSON._print(STDOUT, State(indent), a)
    if indent > 0
        println()
    end
end

json(a, indent=0) = sprint(JSON.print, a, indent)

function determine_bracket_type(io::IO)
    open_bracket = close_bracket = nothing
    while open_bracket == nothing
        eof(io) && throw(EOFError())
        c = read(io, Char)
        if c == '{'
            open_bracket = '{'
            close_bracket = '}'
        elseif c == '['
            open_bracket = '['
            close_bracket = ']'
        elseif c == '\0'
            throw(EOFError())
        end
    end
    open_bracket, close_bracket
end

###
# Consume a string (even if it is invalid), with ack to Douglas Crockford.
# On entry we must already have consumed the opening quotation double-quotation mark
# Add the characters of the string to obj
function consumeString(io::IO, obj::IOBuffer)
    c = '"'

    # When parsing for string values, we must look for " and \ characters.
    while true
        eof(io) && throw(EOFError())
        c = read(io, Char)
        if c == '"'
            write(obj, c)
            return
        end
        if c == '\\'
            write(obj, c)
            eof(io) && throw(EOFError())
            c = read(io, Char)
        end
        write(obj, c)
    end
    throw(EOFError())
end

function parse(io::IO; ordered::Bool=false)
    open_bracket = close_bracket = nothing
    try
        open_bracket, close_bracket = determine_bracket_type(io)
    catch exception
        isa(exception, EOFError) && return
    end
    num_brackets_needed = 1

    obj = IOBuffer()
    write(obj, open_bracket)

    while num_brackets_needed > 0
        eof(io) && throw(EOFError())
        c = read(io, Char)
        write(obj, c)

        if c == open_bracket
            num_brackets_needed += 1
        elseif c == close_bracket
            num_brackets_needed -= 1
        elseif c == '"'
            consumeString(io, obj)
        end
    end
    JSON.parse(takebuf_string(obj), ordered=ordered)
end

function parsefile(filename::AbstractString; ordered::Bool=false, use_mmap=true)
    sz = filesize(filename)
    open(filename) do io
        s = use_mmap ? UTF8String(mmap_array(Uint8, (sz,), io)) : readall(io)
        JSON.parse(s, ordered=ordered)
    end
end

end # module

ext: jl


text:
using Base: @get!

@inline tuple_va(N, xs) = xs
@inline tuple_va(N, x, xs...) = (x, tuple_va(N, xs...)...)
@inline tuple_va(::Val{N}, ::Nothing) where N = ntuple(_ -> nothing, Val(N))

iscall(x, m::Module, n::Symbol) = isexpr(x, :call) && x.args[1] == GlobalRef(m, n)

function isassert(ir, i)
  ex = ir.stmts[i+3]
  iscall(ex, Zygote, :typeassert)
end

# TODO: Move this to Base
function append_node!(ir, @nospecialize(typ), @nospecialize(node), line)
  @assert isempty(ir.new_nodes)
  push!(ir.stmts, node)
  push!(ir.types, typ)
  push!(ir.lines, line)
  push!(ir.flags, 0)
  last_bb = ir.cfg.blocks[end]
  ir.cfg.blocks[end] = BasicBlock(StmtRange(first(last_bb.stmts):length(ir.stmts)),
      last_bb.preds,
      last_bb.succs)
  return SSAValue(length(ir.stmts))
end

function merge_returns(ir)
  any(x -> x == unreachable, ir.stmts) && error("Unsupported control flow")
  rs = findall(x -> x isa ReturnNode && isdefined(x, :val), ir.stmts)
  length(rs) <= 1 && return ir
  bs = blockidx.(Ref(ir), rs)
  xs = Any[]
  bb = length(ir.cfg.blocks)+1
  @assert length(unique(bs)) == length(bs)
  push!(ir.cfg.blocks, BasicBlock(StmtRange(length(ir.stmts)+1, length(ir.stmts)), bs, []))
  push!(ir.cfg.index, length(ir.stmts) + 1)
  r = append_node!(ir, Any, nothing, ir.lines[end])
  append_node!(ir, Any, ReturnNode(r), ir.lines[end])
  for r in rs
    x = ir.stmts[r].val
    x = insert_node!(ir, r, Any, x)
    push!(xs, x)
    ir.stmts[r] = GotoNode(bb)
  end
  for b in bs
    push!(ir.cfg.blocks[b].succs, bb)
  end
  ir.stmts[r.id] = PhiNode(bs, xs)
  ir = compact!(ir)
  return ir
end

struct Alpha
  id::Int
end

Base.show(io::IO, x::Alpha) = print(io, "@", x.id)

alpha(x) = x
alpha(x::SSAValue) = Alpha(x.id)

gradindex(x, i) = x[i]
gradindex(::Nothing, i) = nothing
xgetindex(x, i...) = Expr(:call, GlobalRef(Base, :getindex), x, i...)
xgradindex(x, i) = xcall(Zygote, :gradindex, x, i)

function record_branches!(ir::IRCode)
  ir = IncrementalCompact(ir)
  offset = 0
  for (i, x) in ir
    bi = findfirst(x -> x == i+1-offset, ir.ir.cfg.index)
    bi == nothing && continue
    preds = ir.ir.cfg.blocks[bi+1].preds
    length(preds) > 1 || continue
    insert_node_here!(ir, PhiNode(sort(preds), Int8.(1:length(preds))), Int8, ir.result_lines[i])
    offset += 1
  end
  return finish_dc(ir)
end

istrackable(x) =
  x isa GlobalRef && x.mod ∉ (Base, Core) &&
  !(isconst(x.mod, x.name) && typeof(getfield(x.mod, x.name)) <: Union{Function,Type})

function record_globals!(ir::IRCode)
  for i = 1:length(ir.stmts)
    ex = ir[SSAValue(i)]
    # TODO general globalrefs
    if isexpr(ex, :call)
      for j = 1:length(ex.args)
        istrackable(ex.args[j]) || continue
        ex.args[j] = insert_node!(ir, i, Any, xcall(Zygote, :unwrap, ex.args[j]))
      end
    end
  end
  return compact!(ir)
end

ignored_f(f) = f in (GlobalRef(Base, :not_int),
                     GlobalRef(Core.Intrinsics, :not_int),
                     GlobalRef(Core, :(===)),
                     GlobalRef(Core, :apply_type),
                     GlobalRef(Core, :typeof),
                     GlobalRef(Core, :throw),
                     GlobalRef(Base, :kwerr),
                     GlobalRef(Core, :kwfunc))
ignored_f(ir, f) = ignored_f(f)
ignored_f(ir, f::SSAValue) = ignored_f(ir[f])

ignored(ir, ex) = isexpr(ex, :call) && ignored_f(ir, ex.args[1])
ignored(ir, ex::SSAValue) = ignored(ir, ir[ex])

function valid_usages(ir)
  r = Dict()
  for (x, us) in usages(ir)
    x isa Union{SSAValue,Argument} || continue
    us′ = filter(i -> !ignored(ir, i), us)
    isempty(us′) || (r[x] = us′)
  end
  return r
end

reachable(ir) = keys(valid_usages(ir))

# Hack to work around fragile constant prop through overloaded functions
is_literal_getproperty(ex) = iscall(ex, Base, :getproperty) && ex.args[3] isa QuoteNode

# TODO: remove this once we don't mess with type inference
function _forward_type(Ts)
  usetyped || return Any
  all(T -> isconcretetype(T) || T <: DataType, Ts) || return Any
  T = Core.Compiler.return_type(_forward, Tuple{Context,Ts...})
  return T == Union{} ? Any : T
end

isvalidtype(jT, yT) = jT <: Tuple && length(jT.parameters) == 2 && jT.parameters[1] <: yT

function record!(ir::IRCode)
  pushfirst!(ir.argtypes, typeof(_forward), Context)
  xs = reachable(ir)
  for i = 1:length(ir.stmts)
    ex = argmap(x -> Argument(x.n+2), ir[SSAValue(i)])
    isexpr(ex, :new) && (ex = ir[SSAValue(i)] = xcall(Zygote, :__new__, ex.args...))
    is_literal_getproperty(ex) &&
      (ex = ir[SSAValue(i)] = xcall(Zygote, :literal_getproperty, ex.args[2], Val(ex.args[3].value)))
    if isexpr(ex, :call) && !ignored(ir, ex)
      yT = widenconst(types(ir)[i])
      T = _forward_type(exprtype.(Ref(ir), ex.args))
      if yT == Any || isvalidtype(T, yT)
        yJ = insert_node!(ir, i, T, xcall(Zygote, :_forward, Argument(2), ex.args...))
        ir[SSAValue(i)] = xgetindex(yJ, 1)
        insert_node!(ir, i, T == Any ? Any : T.parameters[2], xgetindex(yJ, 2), true)
      else
        yJ = insert_node!(ir, i, Any, xcall(Zygote, :_forward, Argument(2), ex.args...))
        y =  insert_node!(ir, i, Any, xgetindex(yJ, 1))
        J =  insert_node!(ir, i, Any, xgetindex(yJ, 2))
        ir[SSAValue(i)] = xcall(Zygote, :typeassert, y, yT)
      end
    else
      ir[SSAValue(i)] = ex
    end
  end
  ir, m = _compact!(ir)
  return ir, Set(x isa Argument ? Argument(x.n+2) : x for x in rename(xs, m))
end

# Backwards Pass

function reverse_cfg(cfg, perm)
  newidx(i) = invperm(perm)[i]
  CFG([BasicBlock(StmtRange(1,0),newidx.(b.succs),newidx.(b.preds)) for b in cfg.blocks[perm]])
end

function reverse_order(cfg)
  n = length(cfg.blocks)
  perm = n:-1:1
  guess = reverse_cfg(cfg, perm)
  dt = construct_domtree(guess)
  perm[sortperm(1:n, by = x -> dt.nodes[x].level)]
end

struct Primal
  forw::IRCode
  perm::Vector{Int}
  wrt::Set{Any}
  varargs::Union{Int,Nothing}
end

Primal(ir::IRCode, xs, vs) = Primal(ir, reverse_order(ir.cfg), xs, vs)

function Primal(ir::IRCode; varargs = nothing)
  ir = merge_returns(ir)
  forw, xs = record!(record_branches!(record_globals!(ir)))
  Primal(forw, xs, varargs)
end

newblock(pr::Primal, b) = invperm(pr.perm)[b]
oldblock(pr::Primal, b) = pr.perm[b]

function blockinfo(pr::Primal)
  preds(b) = pr.forw.cfg.blocks[b].preds
  info = Dict(b => (phis=Dict(),partials=[],grads=[]) for b in 1:length(pr.forw.cfg.blocks))
  append!(info[1].grads, filter(x -> x isa Argument, pr.wrt))
  for b in 1:length(pr.forw.cfg.blocks), i in pr.forw.cfg.blocks[b].stmts
    ex = pr.forw[SSAValue(i)]
    if ex isa ReturnNode
      ex.val in pr.wrt && push!(info[b].partials, ex.val)
    elseif ex isa PiNode
      (SSAValue(i) in pr.wrt && ex.val in pr.wrt) || continue
      push!(info[b].grads, SSAValue(i))
      push!(info[b].partials, ex.val)
    elseif ex isa PhiNode
      any(x -> x in pr.wrt, ex.values) && push!(info[b].grads, SSAValue(i))
      for (c, x) in zip(ex.edges, ex.values)
        x in pr.wrt && push!(@get!(info[b].phis, c, []), x)
      end
    elseif iscall(ex, Zygote, :_forward)
      y = isassert(pr.forw, i) ? SSAValue(i+3) : SSAValue(i+1)
      push!(info[b].grads, y)
      for x in ex.args[3:end]
        x in pr.wrt && push!(info[b].partials, x)
      end
    end
  end
  worklist = collect(1:length(pr.forw.cfg.blocks))
  while !isempty(worklist)
    b = pop!(worklist)
    for c in preds(b)
      in = union(get(info[b].phis, c, []), setdiff(info[b].partials, info[b].grads))
      out = union(info[c].partials, info[c].grads)
      new = setdiff(in, out)
      if !isempty(new)
        append!(info[c].partials, new)
        c ∉ worklist && push!(worklist, c)
      end
    end
  end
  return info
end

function IRCode(ir::Primal)
  stmts = []
  blocks = []
  newidx(i) = invperm(ir.perm)[i]
  for block in ir.perm
    old = ir.forw.cfg.blocks[block]
    start = length(stmts)+1
    block == length(ir.perm) && push!(stmts, :(Δ()))
    preds, succs = newidx.(old.succs), newidx.(sort(old.preds))
    if isempty(succs)
      push!(stmts, nothing)
    else
      for (i, b) in enumerate(succs[1:end-1])
        push!(stmts, xcall(Base, :(!==), Alpha(range(old)[1]), Int8(i)))
        push!(stmts, GotoIfNot(SSAValue(length(stmts)), b))
      end
      push!(stmts, GotoNode(succs[end]))
    end
    push!(blocks, BasicBlock(StmtRange(start,length(stmts)), preds, succs))
  end
  ir = IRCode(ir.forw, stmts, Any[Any for _ in stmts], Int32[0 for _ in stmts],
              [0x00 for _ in stmts], CFG(blocks), NewNode[])
end

function reverse_ir(pr::Primal)
  ir = IRCode(pr)
  grads = Dict()
  partials = Dict(x => [] for x in pr.wrt)
  phis = Dict()
  for b in pr.perm
    j = ir.cfg.blocks[newblock(pr, b)].stmts[1]
    j = max(j, 2)
    for i in reverse(pr.forw.cfg.blocks[b].stmts)
      ex = pr.forw[SSAValue(i)]
      if ex isa ReturnNode
        ex.val in pr.wrt && push!(partials[ex.val], SSAValue(1))
      elseif ex isa PiNode
        (SSAValue(i) in pr.wrt && ex.val in pr.wrt) || continue
        Δ = insert_node!(ir, j, Any, xcall(Zygote, :accum))
        ir.lines[j] = pr.forw.lines[i]
        grads[SSAValue(i)] = Δ
        push!(partials[ex.val], Δ)
      elseif ex isa PhiNode
        any(x -> x in pr.wrt, ex.values) || continue
        Δ = insert_node!(ir, j, Any, xcall(Zygote, :accum))
        ir.lines[j] = pr.forw.lines[i]
        grads[SSAValue(i)] = Δ
        for (c, x) in zip(ex.edges, ex.values)
          x in pr.wrt || continue
          @assert !haskey(phis, (newblock(pr, b), newblock(pr, c), x)) "not implemented"
          phis[(newblock(pr, b), newblock(pr, c), x)] = Δ
        end
      elseif iscall(ex, Zygote, :_forward)
        # TODO remove with type hacks above
        y = isassert(pr.forw, i) ? SSAValue(i+3) : SSAValue(i+1)
        J = Alpha(i+2)
        dy = insert_node!(ir, j, Any, xcall(Zygote, :accum))
        ir.lines[j] = pr.forw.lines[i]
        dxs = insert_node!(ir, j, Any, Expr(:call, J, dy))
        ir.lines[j] = pr.forw.lines[i]
        grads[y] = dy
        for (i, x) in enumerate(ex.args[3:end])
          x in pr.wrt || continue
          dx = insert_node!(ir, j, Any, xgradindex(dxs, i))
          ir.lines[j] = pr.forw.lines[i]
          push!(partials[x], dx)
        end
      end
    end
    if b == 1
      gs = []
      for i = 3:length(pr.forw.argtypes)
        Argument(i) in pr.wrt || (push!(gs, nothing); continue)
        dx = insert_node!(ir, j, Any, xcall(Zygote, :accum))
        grads[Argument(i)] = dx
        push!(gs, dx)
      end
      if pr.varargs == nothing
        Δ = insert_node!(ir, j, Any, xcall(Zygote, :tuple, gs...))
      else
        Δ = insert_node!(ir, j, Any, xcall(Zygote, :tuple_va, Val(pr.varargs), gs...))
      end
      insert_node!(ir, j, Any, ReturnNode(Δ))
    end
  end
  ir, m = _compact!(ir)
  return ir, rename(grads, m), rename(partials, m), rename(phis, m)
end

function simplify!(ir)
  ir = IncrementalCompact(ir)
  for (i, x) in ir
    iscall(x, Zygote, :accum) || continue
    filter!(x -> x != nothing, x.args)
    nargs = length(x.args)-1
    ir[i] = nargs == 0 ? nothing : nargs == 1 ? x.args[end] : x
  end
  return finish(ir)
end

function accumulators!(pr::Primal, ir::IRCode, grads, partials, phis)
  blockpartials(b, x) = filter(x -> x.id in ir.cfg.blocks[b].stmts, get(partials, x, []))
  accums = Dict()
  info = blockinfo(pr)
  for b = 1:length(ir.cfg.blocks), x in setdiff(info[b].partials, info[b].grads)
    ps = blockpartials(newblock(pr, b), x)
    p = insert_blockend!(ir, newblock(pr, b), Any, xcall(Zygote, :accum, ps...))
    setdiff!(partials[x], ps)
    push!(partials[x], p)
    accums[(newblock(pr, b),x)] = p
  end

  function predpartial(b, x)
    function blockpartial(b, c, x)
      if haskey(accums, (b, x))
        @assert !haskey(phis, (b, c, x)) "not implemented"
        return accums[(b, x)]
      elseif haskey(phis, (b, c, x))
        return phis[b, c, x]
      end
    end
    preds = ir.cfg.blocks[b].preds
    isempty(preds) && return
    ps = map(c -> blockpartial(c, b, x), preds)
    all(==(nothing), ps) && return
    length(ps) == 1 ? ps[1] : insert_blockstart!(ir, b, Any, PhiNode(preds, ps))
  end

  for ((b, x), p) in accums
    push!(ir[p].args, predpartial(b, x))
  end
  for (x, dx) in grads
    b = blockidx(ir, dx)
    append!(ir[dx].args, blockpartials(b, x))
    push!(ir[dx].args, predpartial(b, x))
  end
  return simplify!(ir)
end

struct Adjoint
  forw::IRCode
  back::IRCode
  perm::Vector{Int}
end

function Adjoint(pr::Primal)
  back = accumulators!(pr, reverse_ir(pr)...)
  Adjoint(pr.forw, compact!(compact!(back)), pr.perm)
end

Adjoint(ir::IRCode; varargs = nothing) = Adjoint(Primal(ir, varargs = varargs))

using InteractiveUtils: @which

macro code_adjoint(ex)
  :(Adjoint($(code_irm(ex)), varargs = varargs($(esc(:(@which $ex))), length(($(esc.(ex.args)...),)))))
end


text:
function setref!(x::EXPR, binding::Binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
    binding.refs !== nothing && push!(binding.refs, x)
end

function setref!(x::EXPR, binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
end

# Main function to be called
function _resolve_ref(x, state)
    if !(parentof(x) isa EXPR && typof(parentof(x)) == CSTParser.Quotenode)
        resolved = resolve_ref(x, state.scope, state)
        if !resolved && (state.delayed || isglobal(valof(x), state.scope))
            push!(state.urefs, x)
        end
    end
end

function resolve_ref(x::EXPR, m, state::State, visited_scopes = 0)
    return hasref(x)
end

function resolve_ref(x1::EXPR, m::SymbolServer.ModuleStore, state::State, visited_scopes = 0)
    hasref(x1) && return true
    if isidentifier(x1)
        x = x1
        if valof(x) == m.name
            setref!(x, m)
            return true
        elseif valof(x) in m.exported && haskey(m.vals, valof(x))
            val = m.vals[valof(x)]
            if val isa SymbolServer.PackageRef
                val1 = SymbolServer._lookup(val, getsymbolserver(state.server))
                if val1 !== nothing
                    setref!(x, val1)
                    return true
                else
                    return false
                end
            else
                setref!(x, val)
                return true
            end
        end
    elseif typof(x1) === MacroName
        x = x1.args[2]
        mn = string("@", valof(x))
        if mn in m.exported
            setref!(x, m.vals[mn])
            return true
        end
    elseif typof(x1) === x_Str
        mac = x1
        mn = string("@", valof(mac.args[1]), "_str")
        if mn in m.exported && haskey(m.vals, mn)
            setref!(mac.args[1], m.vals[mn])
            return true
        end
    end
    return false
end

function resolve_ref(x::EXPR, scope::Scope, state::State, visited_scopes = [])
    hasref(x) && return true

    resolved = false
    if (typof(scope.expr) === CSTParser.ModuleH || typof(scope.expr) === CSTParser.BareModule) && CSTParser.length(scope.expr.args) > 1 && CSTParser.typof(scope.expr.args[2]) === IDENTIFIER
        s_m_name = scope.expr.args[2].val isa String ? scope.expr.args[2].val : ""
        if s_m_name in visited_scopes
            return resolved
        else
            push!(visited_scopes, s_m_name)
        end
    end

    if typof(x) === BinaryOpCall && kindof(x.args[2]) === CSTParser.Tokens.DOT
        return resolve_getindex(x, scope, state)
    elseif isidentifier(x)
        if typof(x) === IDENTIFIER
            mn = valof(x)
            x1 = x
        else
            # NONSTDIDENTIFIER, e.g. var"name"
            mn = valof(x.args[2])
            x1 = x
        end
        if (mn == "__source__" || mn == "__module__") && _in_macro_def(x)
            setref!(x, Binding(noname, nothing, nothing, [], nothing, nothing))
            return true
        end
    elseif resolvable_macroname(x)
        x1 = x.args[2]
        mn = string("@", valof(x1))
    elseif typof(x) === x_Str
        if typof(x.args[1]) === IDENTIFIER
            x1 = x.args[1]
            mn = string("@", valof(x1), "_str")
        else
            return false
        end
    elseif typof(x) === CSTParser.Kw
        if typof(x.args[1]) === IDENTIFIER
            setref!(x.args[1], Binding(noname, nothing, nothing, [], nothing, nothing))
        elseif typof(x.args[1]) === BinaryOpCall && kindof(x.args[1].args[2]) === CSTParser.Tokens.DECLARATION && typof(x.args[1].args[1]) === IDENTIFIER
            setref!(x.args[1].args[1], Binding(noname, nothing, nothing, [], nothing, nothing))
        end
        return true
    else
        return true
    end

    if haskey(scope.names, mn)
        setref!(x1, scope.names[mn])
        resolved = true
    elseif scope.modules isa Dict && length(scope.modules) > 0
        for m in scope.modules
            resolved = resolve_ref(x, m[2], state, visited_scopes)
            resolved && return true
        end
    end
    if !resolved && !scope.ismodule && parentof(scope) isa Scope
        return resolve_ref(x, parentof(scope), state, visited_scopes)
    end
    return resolved
end

function resolve_getindex(x::EXPR, scope::Scope, state::State)
    hasref(x) && return true
    resolved = false
    if typof(x.args[1]) === IDENTIFIER
        resolved = resolve_ref(x.args[1], scope, state)
        if resolved && typof(x.args[3]) === Quotenode && typof(x.args[3].args[1]) === IDENTIFIER
            resolved = resolve_getindex(x.args[3].args[1], refof(x.args[1]), state)
        end
    elseif typof(x.args[1]) === BinaryOpCall && kindof(x.args[1].args[2]) === CSTParser.Tokens.DOT
        resolved = resolve_ref(x.args[1], scope, state)
        if resolved && typof(x.args[3]) === Quotenode && typof(x.args[3].args[1]) === IDENTIFIER
            resolved = resolve_getindex(x.args[3].args[1], refof(x.args[1].args[3].args[1]), state)
        end
    end
    return resolved
end

function resolve_getindex(x::EXPR, b::Binding, state::State)
    hasref(x) && return true
    resolved = false
    if b.type isa Binding
        resolved = resolve_getindex(x, b.type.val, state)
    elseif b.val isa SymbolServer.ModuleStore
        resolved = resolve_getindex(x, b.val, state)
    elseif b.val isa EXPR && (typof(b.val) === ModuleH || typof(b.val) === BareModule)
        resolved = resolve_getindex(x, b.val, state)
    elseif b.val isa Binding && b.val.val isa EXPR && (typof(b.val.val) === ModuleH || typof(b.val.val) === BareModule)
        resolved = resolve_getindex(x, b.val.val, state)
    end
    return resolved
end

function resolve_getindex(x::EXPR, parent_type, state::State)
    hasref(x) && return true
    return false
end
function resolve_getindex(x::EXPR, parent_type::EXPR, state::State)
    hasref(x) && return true
    resolved = false
    if CSTParser.isidentifier(x)
        if (typof(parent_type) === ModuleH || typof(parent_type) === BareModule) && scopeof(parent_type) isa Scope
            resolved = resolve_ref(x, scopeof(parent_type), state)
        elseif CSTParser.defines_struct(parent_type)
            if haskey(scopeof(parent_type).names, valof(x))
                setref!(x, scopeof(parent_type).names[valof(x)])
                resolved = true
            end
        end
    end
    return resolved
end

function resolve_getindex(x::EXPR, parent::SymbolServer.SymStore, state::State)
    hasref(x) && return true
    resolved = false
    if CSTParser.isidentifier(x)
        if parent isa SymbolServer.ModuleStore && haskey(parent.vals, valof(x))
            val = parent.vals[valof(x)]
            if val isa SymbolServer.PackageRef
                val = SymbolServer._lookup(val, getsymbolserver(state.server))
                !(val isa SymbolServer.SymStore) && return false
            end
            setref!(x, val)
            resolved = true
        elseif parent isa SymbolServer.DataTypeStore && valof(x) in parent.fields
            fi = findfirst(f->valof(x) == f, parent.fields)
            val = resolve_typeref(parent.ts[fi], state)
            if val !== nothing
                setref!(x, Binding(noname, nothing, val, [], nothing, nothing))
                resolved = true
            end
        end
    end
    return resolved
end

function resolve_typeref(tr::SymbolServer.TypeRef{N}, state::State) where N
    if haskey(getsymbolserver(state.server), tr.mod.name[1])
        m = getsymbolserver(state.server)[tr.mod.name[1]]
    else
        return
    end
    for i = 2:N
        if haskey(m.vals, tr.mod.name[i])
            m = m.vals[tr.mod.name[i]]
        else
            return
        end
    end
    return get(m.vals, tr.name, nothing)
end

resolvable_macroname(x::EXPR) = typof(x) === MacroName && length(x.args) == 2 && isidentifier(x.args[2]) && refof(x.args[2]) === nothing

# function hasref(x::EXPR)
#     refof(x) !== nothing && refof(x) !== NoReference
# end

function _in_macro_def(x::EXPR)
    if typof(x) === CSTParser.Macro
        return true
    elseif parentof(x) isa EXPR
        return _in_macro_def(parentof(x))
    else
        return false
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

function get_test_model()
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 0.5],
        4 => [0.0, 0.5],
        5 => [0.0, 0.6],
        6 => [1.0, 0.6],
        7 => [1.0, 1.1],
        8 => [0.0, 1.1])
    el1 = Element(Quad4, [1, 2, 3, 4])
    el2 = Element(Quad4, [5, 6, 7, 8])
    el3 = Element(Seg2, [1, 2])
    el4 = Element(Seg2, [7, 8])
    el5 = Element(Seg2, [4, 3])
    el6 = Element(Seg2, [5, 6])
    update!([el1, el2, el3, el4, el5, el6], "geometry", X)
    update!([el1, el2], "youngs modulus", 96.0)
    update!([el1, el2], "poissons ratio", 1/3)
    update!([el3], "displacement 1", 0.0)
    update!([el3], "displacement 2", 0.0)
    update!([el4], "displacement 1", 0.0)
    update!([el4], "displacement 2", 0.0)
    update!(el6, "master elements", [el5])
    p1 = Problem(Elasticity, "body1", 2)
    p2 = Problem(Elasticity, "body2", 2)
    p3 = Problem(Dirichlet, "fixed", 2, "displacement")
    p4 = Problem(Mortar, "interface", 2, "displacement")
    push!(p1, el1)
    push!(p2, el2)
    push!(p3, el3, el4)
    push!(p4, el5, el6)
    return p1, p2, p3, p4
end

@testset "test adjust setting in 2d tie contact" begin
    p1, p2, p3, p4 = get_test_model()
    p1.properties.formulation = :plane_stress
    p2.properties.formulation = :plane_stress
    p4.properties.adjust = true
    p4.properties.rotate_normals = false
    solver = Solver(Linear)
    push!(solver, p1, p2, p3, p4)
    solver()
    el5 = p4.elements[1]
    u = el5("displacement", [0.0], 0.0)
    info("u = $u")
    @test isapprox(u, [0.0, 0.05])
end

@testset "test that interface transfers constant field without error" begin
    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d.med"
    mesh = aster_read_mesh(meshfile)

    upper = Problem(Heat, "upper", 1)
    upper.properties.formulation = "2D"
    upper.elements = create_elements(mesh, "UPPER")
    update!(upper.elements, "temperature thermal conductivity", 1.0)

    lower = Problem(Heat, "lower", 1)
    lower.properties.formulation = "2D"
    lower.elements = create_elements(mesh, "LOWER")
    update!(lower.elements, "temperature thermal conductivity", 1.0)

    bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
    update!(bc_upper.elements, "temperature 1", 0.0)

    bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
    update!(bc_lower.elements, "temperature 1", 1.0)

    interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
    update!(interface_slave_elements, "master elements", interface_master_elements)
    interface.elements = [interface_master_elements; interface_slave_elements]

    solver = Solver(Linear)
    push!(solver, upper, lower, bc_upper, bc_lower, interface)
    solver()

    #interface_norm = norm(interface.assembly)
    # for bi-orthogonal:
    #interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.44870723441585775, 0.44870723441585775, 0.0, 0.0, 0.0]
    #interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.39361633468943247, 0.39361633468943247, 0.0, 0.0, 0.0]
    #info("Interface norm:          $interface_norm")
    #info("Interface norm expected: $interface_norm_expected")
    #@test isapprox(interface_norm, interface_norm_expected)

    T_upper = first(bc_upper.elements)("temperature", [0.0], 0.0)
    T_lower = first(bc_lower.elements)("temperature", [0.0], 0.0)
    T_middle = first(interface.elements)("temperature", [0.0], 0.0)
    info("T upper: $T_upper, T lower: $T_lower, T interface: $T_middle")

    node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
    T = [t[1] for t in temperature]
    minT = minimum(T)
    maxT = maximum(T)
    info("minT = $minT, maxT = $maxT")
    @test isapprox(minT, 0.5)
    @test isapprox(maxT, 0.5)
end

#=
# TODO: if one forget plane_stress solver gives singular exception and it's
    hard to trace to the source of problem
@testset "expect clear error when trying to solve 2d model in 3d setting" begin
    p1, p2, p3, p4 = get_test_model()
#   p1.properties.formulation = :plane_stress
#   p2.properties.formulation = :plane_stress
    p4.properties.adjust = true
    p4.properties.rotate_normals = false
    solver = Solver(Nonlinear)
    solver.properties.linear_system_solver = :DirectLinearSolver_UMFPACK
    push!(solver, p1, p2, p3, p4)
    solver()
    el5 = p4.elements[1]
    u = el5("displacement", [0.0], 0.0)
    info("u = $u")
    @test isapprox(u, [0.0, 0.05])
end
=#

function JuliaFEM.get_mesh(::Type{Val{Symbol("1x1 block splitted to upper and lower")}})
    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d.med"
    mesh = aster_read_mesh(meshfile)
end

function JuliaFEM.get_model(::Type{Val{Symbol("splitted block, plane stress elasticity and mesh tie")}})
    mesh = get_mesh("1x1 block splitted to upper and lower")

    upper = Problem(Elasticity, "upper", 2)
    upper.properties.formulation = :plane_stress
    upper.elements = create_elements(mesh, "UPPER")
    update!(upper.elements, "youngs modulus", 100.0)
    update!(upper.elements, "poissons ratio", 1/3)

    lower = Problem(Elasticity, "lower", 2)
    lower.properties.formulation = :plane_stress
    lower.elements = create_elements(mesh, "LOWER")
    update!(lower.elements, "youngs modulus", 100.0)
    update!(lower.elements, "poissons ratio", 1/3)

    bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
#   update!(bc_upper.elements, "displacement 1", 0.1)
    update!(bc_upper.elements, "displacement 2", -0.1)

    bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
#   update!(bc_lower.elements, "displacement 1", 0.0)
    update!(bc_lower.elements, "displacement 2", 0.0)

    bc_corner = Problem(Dirichlet, "fix model from lower left corner to prevent singularity", 2, "displacement")
    node_ids = find_nearest_nodes(mesh, [0.0, 0.0])
    bc_corner.elements = [Element(Poi1, node_ids)]
    update!(bc_corner.elements, "geometry", mesh.nodes)
    update!(bc_corner.elements, "displacement 1", 0.0)

    interface = Problem(Mortar, "interface between upper and lower block", 2, "displacement")
    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
    update!(interface_slave_elements, "master elements", interface_master_elements)
    interface.elements = [interface_master_elements; interface_slave_elements]

    solver = Solver(Linear)
    push!(solver, upper, lower, bc_upper, bc_lower, interface, bc_corner)

    return solver

end

@testset "test mesh tie with splitted block and plane stress elasticity" begin
    solver = get_model("splitted block, plane stress elasticity and mesh tie")
    upper, lower, bc_upper, bc_lower, interface = solver.problems
    solver()
    slave_elements = get_slave_elements(interface)
    node_ids, la = get_nodal_vector(slave_elements, "reaction force", 0.0)
    for lai in la
        @test isapprox(lai, [0.0, -10.0])
    end
    # FIXME
end

function JuliaFEM.get_mesh(::Type{Val{Symbol("curved 2d block splitted to upper and lower")}})
    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med"
    mesh = aster_read_mesh(meshfile)
end

function JuliaFEM.get_model(::Type{Val{Symbol("mesh tie with curved 2d block")}};
            dy=0.0, adjust=false, tolerance=0.0, rotate_normals=false, swap=false,
            dual_basis=false, use_forwarddiff=false)

    mesh = get_mesh("curved 2d block splitted to upper and lower")

    upper = Problem(Elasticity, "upper", 2)
    upper.properties.formulation = :plane_stress
    upper.elements = create_elements(mesh, "UPPER")
    update!(upper.elements, "youngs modulus", 96.0)
    update!(upper.elements, "poissons ratio", 1/3)

    lower = Problem(Elasticity, "lower", 2)
    lower.properties.formulation = :plane_stress
    lower.elements = create_elements(mesh, "LOWER")
    update!(lower.elements, "youngs modulus", 96.0)
    update!(lower.elements, "poissons ratio", 1/3)

    bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
    update!(bc_upper.elements, "displacement 1", 0.0)
    update!(bc_upper.elements, "displacement 2", dy)

    bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
    update!(bc_lower.elements, "displacement 1", 0.0)
    update!(bc_lower.elements, "displacement 2", 0.0)

    interface = Problem(Mortar, "interface between upper and lower block", 2, "displacement")
    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
    if swap
        interface_slave_elements, interface_master_elements = interface_master_elements, interface_slave_elements
    end
    update!(interface_slave_elements, "master elements", interface_master_elements)
    interface.elements = [interface_master_elements; interface_slave_elements]
    interface.properties.adjust = adjust
    interface.properties.distval = tolerance
    interface.properties.rotate_normals = rotate_normals
    interface.properties.dual_basis = dual_basis
    interface.properties.use_forwarddiff = use_forwarddiff

    solver = Solver(Linear)
    push!(solver, upper, lower, bc_upper, bc_lower, interface)

    return solver

end

@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=0.0" begin
    # TODO: analytical solution now known, verify using other fem software
    solver = get_model("mesh tie with curved 2d block";
        adjust=true, tolerance=10, dy=0.0, rotate_normals=true,
        dual_basis=false)
    solver()
    interface = solver["interface between upper and lower block"]
    @test isapprox(norm(interface.assembly.u), 0.11339715157447851)
end

@testset "curved surface with adjust=true, dual lagrange, slave=lower surface, dy=0.0" begin
    # TODO: analytical solution now known, verify using other fem software
    solver = get_model("mesh tie with curved 2d block";
        adjust=true, tolerance=10, dy=0.0, rotate_normals=true,
        dual_basis=true)
    solver()
    interface = solver["interface between upper and lower block"]
    @test isapprox(norm(interface.assembly.u), 0.11660422877751599)
end

@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=-0.1" begin
    # TODO: analytical solution now known, verify using other fem software
    solver = get_model("mesh tie with curved 2d block";
        adjust=true, tolerance=10, dy=-0.1, rotate_normals=true,
        dual_basis=false)
    solver()
    interface = solver["interface between upper and lower block"]
    @test isapprox(norm(interface.assembly.u), 0.34230262165505887)
end

@testset "curved surface, adjust=true, dual basis, slave=lower surface, dy=-0.1" begin
    # TODO: analytical solution now known, verify using other fem software
    solver = get_model("mesh tie with curved 2d block";
        adjust=true, tolerance=10, dy=-0.1, rotate_normals=true,
        dual_basis=true)
    solver()
    interface = solver["interface between upper and lower block"]
    @test isapprox(norm(interface.assembly.u), 0.34318800698017704)
end


text:
###
# Sampler is currently implemented equally for  Unordered_Design and Ordered_Design base on index.
# Can be improved, special for Unordered_Design but for small space it should be oke
###

####
# implementaiton base on the StatsBase version [StatsBase](https://juliastats.org/StatsBase.jl/stable/sampling/)
####
"""
    StatsBase.sample!(rng::AbstractRNG, space::Eff_Space,x::AbstractArray; with_index::Bool = false , replace::Bool=false, ordered::Bool=false)

Extend the StatsBase.sample! function to sample for a`Eff_Space`.
Using the default settings, It fills `x` with unique samples from the `Eff_Space`. It returns a vector with the corresponding samples.
If with_index = true, corresponding indexes returned together with the sampled constructs.
For more information see [StatsBase](https://juliastats.org/StatsBase.jl/stable/sampling/)
"""

function StatsBase.sample!(rng::AbstractRNG, space::Eff_Space,x::AbstractArray; with_index::Bool = false , replace::Bool=false, ordered::Bool=false)
    index = sample!(rng,1:length(space),x;replace=replace,ordered=ordered)
    if with_index
        return [[space.space[i] for i in index]  index]
    else
        return ([space.space[i] for i in index])
    end
end

# if no RNG object is given
StatsBase.sample!(a::Eff_Space, x::AbstractArray; with_index::Bool = false, replace::Bool=false, ordered::Bool=false) =
    sample!(Random.GLOBAL_RNG, a, x; whit_index = with_index, replace=replace, ordered=ordered)

"""
    StatsBase.sample(rng::AbstractRNG, a::Eff_Space, n::Integer; with_index::Bool = false,replace::Bool=false, ordered::Bool=false)

Extend the StatsBase.Sample function to sample for a`Eff_Space`.
Using the default settings,  It draws `n` unique constructs form the `Eff_Space`.
It returns a vector with the corresponding samples.
If with_index = true, corresponding indexes are returned together with the sampled constructs.
For more information see [StatsBase](https://juliastats.org/StatsBase.jl/stable/sampling/)
"""
function StatsBase.sample(rng::AbstractRNG, a::Eff_Space, n::Integer; with_index::Bool = false,
                replace::Bool=false, ordered::Bool=false)
    sample!(rng, a, Vector{Int}(undef, n); with_index = with_index, replace=replace, ordered=ordered)
end

# if no RNG object is given
StatsBase.sample(a::Eff_Space, n::Integer;with_index::Bool = false, replace::Bool=false, ordered::Bool=false) =
    sample(Random.GLOBAL_RNG, a, n; whit_index = with_index, replace=replace, ordered=ordered)

"""
    sample(rng::AbstractRNG,space::Frame_Space,n::Int)

Return Array with dimension (n x 2).  The function samples `n` constructs from the input `space`, these are stored in the first collum.
The second column returns the corresponding indices of the construct in the closest efficient design space.
These are not the correct indexes in the constrained space.
Currently not sure if returning them is useful. The with_index = true is to be similar to the other sampler, has no influence on the output

The sampling is based on a recursive function with a reject the unwanted constructs.
These are the constructs which are not allowed based on the given constraints.
For design space with many constraints, it is not efficient, and explicit calculation of the design space may be a better alternative using
`getspace(Frame Space, full = true) `

"""
function StatsBase.sample(rng::AbstractRNG,space::Frame_Space,n::Int; with_index = true)
    sample_reject!(rng,space.space,n,space.con,Array{Int}(undef, 0, 1),  Array{eltype(space)}(undef, 0, 2))
end

"""
    sample_reject!(rng::AbstractRNG,space::Eff_Space,n::Int,con,save_index::Array,check_constructs::Array)

Return Array with dimension (n x 2).  The function samples `n` constructs from the input `space` and corresponding indexes.
These constructs are added to the  `check_constructs` array containing the previously sampled constructs.
The `save_index` contains al previous evaluated constructs and prevents resampling  and the revaluation of an unwanted constructs.

The sampling is based on a recursive function which rejects the unwanted constructs.
These are the constructs which are not allowed based on the given constraints.
For design space with many constraints, it is not efficient, and explicit calculation of the design space may be a better alternative using
`getspace(Frame Space, full = true) `

[`sample(rng::AbstractRNG,space::Frame_Space,n::Int)`]@ref

"""
# build to make maximum use of the StatsBase implementaiton to sample multiple unique samples, sequential sampler is avoided
function sample_reject!(rng::AbstractRNG,space::Eff_Space,n::Int,con,save_index::Array,check_constructs::Array)

            # number of started constructs
            n_start = size(check_constructs)[1]
            #make new samples
            new_sample = sample(rng,space,n,with_index = true)
            #check if draw before
            new_sample = _filter_sample!(new_sample,save_index)
            #updated check indexes
            save_index = [save_index;new_sample[:,2]]
            #check constrain
            new_sample = _filter_sample!(new_sample,con)
            #updata approved constructs
            check_constructs = [check_constructs;new_sample]
            #number of constructs ad the end
            n_stop = size(check_constructs)[1]
            # evaluated if their where rejections, if not return else resample for the remaining open positions
            if  n_stop-n_start == n
                return (check_constructs)

            else
                delta = n - (n_stop - n_start)
                sample_reject!(rng,space,delta,con,save_index,check_constructs)
            end
end

"""
    _filter_sample!(new_sample,con::Construct_Constrains)

Internal function to sample_reject!. Evaluates all samples based on the given constraints and removes the unallowed constructs.

"""
_filter_sample!(new_sample,con::Construct_Constrains) = map(y -> !filter_constrain(y,con),new_sample[:,1]) |> (y -> new_sample[y,:])

"""
    _filter_sample!(new_sample,saved_index::Array)

 Evaluates is all samples were sampled or evaluated before and remove them if this is true

"""
_filter_sample!(new_sample,saved_index::Array) = map((y -> y ∉ saved_index),new_sample[:,2]) |> (y -> new_sample[y,:])

"""
    StatsBase.sampleStatsBase.sample(rng::AbstractRNG, space::Multi_Space, n::Integer; with_index::Bool = true)

Return Array with dimension (n x 2).  The function samples `n`  unique constructs from the input `space` and corresponding indexes.
The function sampless from the underlying space that constructs the `Multi_space`.

**note** For `Multi_Space{Frame_Space}` these indexes are form the whole design space with considering the constrains.
The sampling is based on a recursive function which rejects the unwanted constructs.
These are the constructs which are not allowed based on the given constraints.
For design space with many constraints, it is not efficient, and explicit calculation of the design space may be a better alternative using
`getspace(Frame Space, full = true) `

[`sample(rng::AbstractRNG,space::Frame_Space,n::Int)`](@ref)
[sample(rng::AbstractRNG,space::Eff_Space,n::Int)](@ref)

"""

function StatsBase.sample(rng::AbstractRNG, space::Multi_Space, n::Integer; with_index::Bool = true)
    #prelocated
    samples = Array{T where T,2}(undef,n,2)
    #set weigths of the sampler

    w,max = [length(s) for s in space.space] |> x ->(FrequencyWeights(x),sum(x))
    println(max)
    @assert n < max "can only sample $max points"
    # one point
    samples[1,:] = sample(rng,collect(1:_nspace(space)), w,1)  |> x -> space.space[x[1]] |> y -> (sample(rng,y,1,with_index = true))
    #sampels form whole the design space n times
    #sampler one selects with space is chosen, this is weighted, sample step 2 samples form that space
    if n == 1
         return samples
    end
    for i in 2:n
        new_point = sample(rng,collect(1:_nspace(space)), w,1) |> x -> space.space[x[1]] |> y -> (sample(rng,y,1,with_index = true))
        # select n unique sample, evaluated if draw before

        while  new_point[2] in samples[1:i-1,2]
            new_point = sample(rng,collect(1:size(space)[2]), w,1) |> x -> space.space[x[1]] |> y -> (sample(rng,y,1,with_index = true))
        end

        samples[i,:] = new_point
    end

    return samples
end


text:
@testset "MLJ" begin
  Z = [ones(10,10) 2ones(10,10); 3ones(10,10) 4ones(10,10)]
  𝒮 = georef((Z=Z,))

  # non-probabilistic model
  kmeans = @load KMeans pkg=Clustering verbosity=0
  C = cluster(𝒮, kmeans(k=4))
  @test Set(C.cluster) == Set([1,2,3,4])

  # probabilistic model
  gmm = @load GMMClusterer pkg=BetaML verbosity=0
  C = cluster(𝒮, gmm(K=4))
  # see https://github.com/alan-turing-institute/MLJ.jl/issues/846
  @test_broken Set(C.cluster) == Set([1,2,3,4])
end


text:
using LowRankModels

# loss types to test
real_loss_types = [QuadLoss, HuberLoss]
bool_loss_types = [HingeLoss]
ordinal_loss_types = [OrdinalHingeLoss, BvSLoss]
categorical_loss_types = [MultinomialLoss, OvALoss]

#instantiate losses
ncat = 4 # maximum categorical levels
nord = 5 # maximum ordinal levels
real_losses = [l() for l in real_loss_types]
bool_losses = [l() for l in bool_loss_types]
ordinal_losses = [l(rand(3:nord)) for l in ordinal_loss_types]
categorical_losses = [l(rand(3:ncat)) for l in categorical_loss_types]
losses = [real_losses..., bool_losses..., ordinal_losses..., categorical_losses...]

# scale losses for different columns
for loss in losses
  scale!(loss, rand())
end

# regularizers to test
regularizers = [QuadReg(), OneReg(5), NonNegConstraint(), KSparseConstraint(2)]

m,n = length(regularizers), length(losses)

A_real = rand(m, length(real_losses))
A_bool = rand(Bool, m, length(bool_losses))
A_ord = rand(1:5, m, length(ordinal_losses))
A_cat = rand(1:3, m, length(categorical_losses))

# without saying "Any", upconverts to array of Floats
A = Any[A_real A_bool A_ord A_cat] # XXX upconverts to Array{Float64,2}

glrm = GLRM(A, losses, regularizers, QuadReg(), 2)
fit!(glrm)
println("successfully fit matrix")

### now fit data frame

using DataFrames
A = DataFrame(A)
for i=1:10
  A[rand(1:m), rand(1:n)] = NA
end
obs = observations(A)
glrm = GLRM(A, losses, QuadReg(), QuadReg(), 2, obs=obs)
fit!(glrm)
println("successfully fit dataframe")

# imputation and sampling
impute(glrm)
println("successfully imputed entries")
sample(glrm)
sample_missing(glrm)
println("successfully sampled from model")

### now fit sparse matrix

m, n = 100, 200
A = sprandn(m, n, .5)
glrm = GLRM(A, QuadLoss(), QuadReg(), QuadReg(), 5)
fit!(glrm)
println("successfully fit sparse GLRM")


text:
config = MOI.Test.TestConfig()
_optimize!(mock) = MOIU.mock_optimize!(mock, [4.0, 1.0])
for mock in mocks(_optimize!)
    Tests.dsos_quartic_constant_test(mock, config)
    Tests.sdsos_quartic_constant_test(mock, config)
    Tests.sos_quartic_constant_test(mock, config)
end


text:
"""
    compute_frames!(elements::Vector{UA}; parent=nothing)
        where UA<:Union{AbstractObject,AbstractAction}

Set elem.frames.frames to the computed frames for each elem in elements.
"""
function compute_frames!(
    elements::Vector{UA};
    parent = nothing,
    parent_counter = 0,
) where {UA<:Union{AbstractObject,AbstractAction}}
    available_subframes = typemin(Int):typemax(Int)
    if parent !== nothing
        last_frames = get_frames(parent)
        available_subframes = 1:length(get_frames(parent))
    else
        last_frames = nothing
    end
    is_first = true
    counter = 1
    for elem in elements
        # update (CURRENT/PREVIOUS)_(OBJECT/ACTION) and CURRENT_OBJECT_ACTION_TYPE
        if elem isa Object
            set_current_object(elem)
            empty!(CURRENT_ACTION)
        elseif elem isa Action
            set_current_action(elem)
        end
        set_current_action_type(elem)

        if last_frames === nothing && get_frames(elem) === nothing
            throw(ArgumentError("Frames need to be defined explicitly in the initial
                Object/Background or Action."))
        end
        if get_frames(elem) === nothing
            set_frames!(parent, elem, last_frames; is_first = is_first)
        end
        last_frames = get_frames(elem)
        if !(get_frames(elem) ⊆ available_subframes)
            @warn(
                "Action defined outside the frame range of the parent object.
          Action #$counter for Object #$parent_counter is defined for frames
          $(get_frames(elem)) but Object #$parent_counter exists only for $(available_subframes).
          (Info: Background is counted as Object #1)"
            )
        end
        is_first = false
        counter += 1
    end
end

"""
    get_current_setting()

Return the current setting of the current object
"""
function get_current_setting()
    object = CURRENT_OBJECT[1]
    return object.current_setting
end

"""
    get_interpolation(frames::UnitRange, frame)

Return a value between 0 and 1 which represents the relative `frame` inside `frames`.
"""
function get_interpolation(frames::UnitRange, frame)
    frame == last(frames) && return 1.0
    t = (frame - first(frames)) / (length(frames) - 1)
    # makes sense to only allow 0 ≤ t ≤ 1
    t = min(1.0, t)
end

"""
    get_interpolation(action::AbstractAction, frame)

Return the value of the `action.anim` Animation based on the relative frame given by
`get_interpolation(get_frames(action), frame)`
"""
function get_interpolation(action::AbstractAction, frame)
    t = get_interpolation(get_frames(action), frame)
    if !(action.anim.frames[end].t ≈ 1)
        @warn "Animations should be defined from 0.0 to 1.0"
    end
    return interpolation_to_transition_val(at(action.anim, t), action.transition)
end

"""
    interpolation_to_transition_val(interpolation_val, Transition)

Returns the transition value for the given `interpolation_val`.
If the interpolation value is already of the correct form it just gets returned.
Otherwise the Transition function like `get_position` is called and the interpolated value
is calculated.
"""
interpolation_to_transition_val(t, ::Nothing) = t
interpolation_to_transition_val(t::Point, trans::Translation) = t
interpolation_to_transition_val(t::Float64, trans::Rotation) = t
interpolation_to_transition_val(t::Scale, trans::Scaling) = t

function interpolation_to_transition_val(t, trans::Translation)
    # does interpolate between `to` and `from` and assumes we are at `from` already
    from = get_position(trans.from)
    to = get_position(trans.to)
    return t * (to - from)
end

function interpolation_to_transition_val(t, trans::Scaling)
    from = get_scale(trans.from)
    to = get_scale(trans.to)
    return from + t * (to - from)
end

function isapprox_discrete(val; atol = 1e-4)
    return isapprox(val, round(val); atol = atol)
end

function polywh(polygon::Vector{Vector{Point}})
    T = typeof(polygon[1][1].x)
    min_x = typemax(T)
    min_y = typemax(T)
    max_x = typemin(T)
    max_y = typemin(T)
    for poly in polygon
        for p in poly
            min_x = min(min_x, p.x)
            min_y = min(min_y, p.y)
            max_x = max(max_x, p.x)
            max_y = max(max_y, p.y)
        end
    end
    return max_x - min_x, max_y - min_y
end

function get_polypoint_at(points, t; pdist = polydistances(points))
    if t ≈ 0
        return points[1]
    end
    ind, surplus = nearestindex(pdist, t * pdist[end])

    nextind = mod1(ind + 1, length(points))
    overshootpoint = between(
        points[ind],
        points[nextind],
        surplus / distance(points[ind], points[nextind]),
    )
    return overshootpoint
end

"""
    set_previous_object(object::Object)

Set the `object` as `PREVIOUS_OBJECT`
"""
function set_previous_object(object::Object)
    if isempty(PREVIOUS_OBJECT)
        push!(PREVIOUS_OBJECT, object)
    else
        PREVIOUS_OBJECT[1] = object
    end
end

"""
    set_previous_action(action::Action)

Set the `action` as `PREVIOUS_ACTION`
"""
function set_previous_action(action::Action)
    if isempty(PREVIOUS_ACTION)
        push!(PREVIOUS_ACTION, action)
    else
        PREVIOUS_ACTION[1] = action
    end
end

"""
    set_current_object(object::Object)

Set the `object` as `CURRENT_OBJECT` and update `PREVIOUS_OBJECT`/`PREVIOUS_ACTION`
"""
function set_current_object(object::Object)
    update_previous_object_or_action()

    if isempty(CURRENT_OBJECT)
        push!(CURRENT_OBJECT, object)
    else
        CURRENT_OBJECT[1] = object
    end
end

"""
    set_current_action(action::Action)

Set the `action` as `CURRENT_ACTION` and update `PREVIOUS_OBJECT`/`PREVIOUS_ACTION`
"""
function set_current_action(action::Action)
    update_previous_object_or_action()

    if isempty(CURRENT_ACTION)
        push!(CURRENT_ACTION, action)
    else
        CURRENT_ACTION[1] = action
    end
end

"""
    update_previous_object_or_action()

Update the `PREVIOUS_OBJECT` or `PREVIOUS_ACTION` depending on whether the
last element was an object or an action. This is still saved in `CURRENT_OBJECT_ACTION_TYPE`.
"""
function update_previous_object_or_action()
    if !isempty(CURRENT_OBJECT_ACTION_TYPE)
        if CURRENT_OBJECT_ACTION_TYPE[1] == :Object
            !isempty(CURRENT_OBJECT) && set_previous_object(CURRENT_OBJECT[1])
        else
            !isempty(CURRENT_ACTION) && set_previous_action(CURRENT_ACTION[1])
        end
    end
end

"""
    set_current_action_type(t)

Set `CURRENT_OBJECT_ACTION_TYPE` to `:Object` or `:Action` depending
on the type of `t`.
"""
function set_current_action_type(t)
    type = :Object
    if t isa AbstractAction
        type = :Action
    end
    if isempty(CURRENT_OBJECT_ACTION_TYPE)
        push!(CURRENT_OBJECT_ACTION_TYPE, type)
    else
        CURRENT_OBJECT_ACTION_TYPE[1] = type
    end
end

"""
    empty_CURRENT_constants()

empty all `CURRENT_` constants like `CURRENT_OBJECT`
"""
function empty_CURRENT_constants()
    empty!(CURRENT_VIDEO)
    empty!(CURRENT_OBJECT)
    empty!(CURRENT_ACTION)
    empty!(PREVIOUS_OBJECT)
    empty!(PREVIOUS_ACTION)
    empty!(CURRENT_OBJECT_ACTION_TYPE)
end


text:
const τ_e = 5
const τ_i = 10
const g_spike = 1

# at the expense of using BLAS + slicing to perform layer updates, every synapse is evaluated separately

function SpikingSynapse(stype::SynapseType, s::Int64, t::Int64, w::Float64, d::Float64)
	tau = (stype == EXC) ? τ_e : τ_i
	g = 0
	_ql = int(d/dt)+1
	incoming = falses(_ql)
	head = 1
	SpikingSynapse(stype,s,t,w,d,g,tau,incoming,head,_ql)
end

function update!(syn::SpikingSynapse, source::AbstractSpikingNeuron)
	syn.g += dt*(-syn.g/syn.tau)
	if syn.incoming[syn._head]
		syn.g += g_spike
	end
	# cycle forward

	# in a 0-indexed language:
	# end = (head + cl-1 ) % cl
	# therefore, in a 1-indexed language (convert to 0-index language, then back):
	# end = ((head-1) + ql - 1) % ql + 1
	e = ((syn._head-1)+syn._ql-1) % syn._ql + 1
	syn.incoming[e] = source.didSpike
	syn._head = syn._head%syn._ql+1
end


text:
module TestMorrisLecar
include("preamble.jl")

using StaticArrays: SVector

using Bifurcations: resolved_points
using Bifurcations.BifurcationsBase: count_special_points
using Bifurcations.Codim2LimitCycle: FoldLimitCycleProblem
using Bifurcations.Examples.MorrisLecar

# To test start_from_nearest_root, let's start from an arbitrary
# place:
prob = MorrisLecar.make_prob(
    MorrisLecar.MorrisLecarParam(y=0.1, z=0.1);
    u0 = SVector(0.0, 0.0),
)
solver1 = init(
    prob;
    start_from_nearest_root = true,
)
@time solve!(solver1)

all_specials = special_intervals(solver1)
all_hopf = special_intervals(solver1, Codim1.PointTypes.hopf)
all_sn = special_intervals(solver1, Codim1.PointTypes.saddle_node)
@test length(all_specials) == 3
@test length(all_sn) == 2
@test length(all_hopf) == 1

hopf_point, = all_hopf
hopf_point

codim2_prob = BifurcationProblem(
    hopf_point,
    solver1,
    (@lens _.z),
    (-1.0, 1.0),
)
hopf_solver = init(
    codim2_prob;
    nominal_angle_rad = 0.01,
)
solve!(hopf_solver)

all_bautin = special_intervals(hopf_solver, Codim2.PointTypes.bautin)
@assert length(all_bautin) == 1
bautin_point, = all_bautin
bautin_point

using Bifurcations.Codim2LimitCycle: FoldLimitCycleProblem

@time flc_prob = FoldLimitCycleProblem(
    bautin_point,
    hopf_solver;
    num_mesh = 30,
    degree = 4,
);
flc_solver = init(
    flc_prob;
    start_from_nearest_root = true,
    max_branches = 0,
    max_samples = 20,
    bidirectional_first_sweep = false,
    nominal_angle_rad = 2π * (5 / 360),
    verbose = true,
)
@time solve!(flc_solver)

sn_prob = BifurcationProblem(
    special_intervals(solver1, Codim1.PointTypes.saddle_node)[1],
    solver1,
    (@lens _.z),
    (-1.0, 1.0),
)
sn_solver = init(
    sn_prob;
    nominal_angle_rad = 2π * (5 / 360),
    max_samples = 1000,
    start_from_nearest_root = true,
)
@time solve!(sn_solver)

@test count_special_points(sn_solver) == Dict(
    Codim2.PointTypes.bogdanov_takens => 1,
    Codim2.PointTypes.cusp => 1,
)

hopf_prob2 = BifurcationProblem(
    special_intervals(sn_solver, Codim2.PointTypes.bogdanov_takens)[1],
    sn_solver,
)
hopf_solver2 = init(hopf_prob2)
@time solve!(hopf_solver2)

count_hopf_solver2 = count_special_points(hopf_solver2)

@test Set(keys(count_hopf_solver2)) <= Set([Codim2.PointTypes.bogdanov_takens])
if Codim2.PointTypes.bogdanov_takens in keys(count_hopf_solver2)
    @test count_hopf_solver2 == Dict(
        Codim2.PointTypes.bogdanov_takens => 1,
    )
end

end  # module


text:
using DataFrames

concentrations = sort([[0:8]//8, 1//16, 15//16, 13//16, 3//16])
target = 800 # MicroLitre

mixtures = map(concentrations) do c
  float(c*target), float(target - c*target)
end

cs = [zip(mixtures...)...]
df = DataFrame(CAcid = [cs[1]...], Water = [cs[2]...])
writetable("mixtures.csv", df)

max_stars_repo_name: vchuravy/TrackingCBoat

max_issues_repo_name: vchuravy/TrackingCBoat

max_forks_repo_name: vchuravy/TrackingCBoat


text:
# This file was generated, do not modify it. # hide
ensemble_model = EnsembleModel(atom=knn_model, n=20);


text:
import GeneticJulia
include("fitnessFunctions.jl")

GeneticJulia.setRandomSeed(5)
GeneticJulia.setIndividualType(GeneticJulia.BinaryGenotype)
GeneticJulia.setStopCondition(maxIterations=500)
GeneticJulia.setEvaluator(GeneticJulia.FitnessFunction(countBinaryOnes))
GeneticJulia.setGenerator(GeneticJulia.randomBinaryGenerator, 50, popSize = 100)
GeneticJulia.setSelector(GeneticJulia.tournamentSelector, 20, samplingWithRep=true, nSelected=2.)
GeneticJulia.setCrossoverOperator(GeneticJulia.uniformCross, nChildren=2)
GeneticJulia.setMutationOperator(GeneticJulia.standardMutation, 0.1)
GeneticJulia.setReplacementOperator(GeneticJulia.replaceAllPopulation)

@time GeneticJulia.runGenJ(verbose=false)

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
type MatrixDataSet{N} <: DataSet
    X::Array{N, 2}
    y::Array{N, 1}
    indepvars::Vector{Symbol}
    depvar::Symbol
end

# A matrix data set is given the X matrix and the y vector explicitly and can (optionally)
# be given names of the independent vars as well as the dependent var.
function MatrixDataSet{N <: Number}(X::Array{N, 2}, y::Array{N, 1}; indepvars = nothing, depvar = :y)
    @assert length(y) == size(X, 2)
    if indepvars == nothing
        indepvars = map(i -> Symbol("x$i"), 1:size(X, 1))
    end
    MatrixDataSet{eltype(X)}(X, y, indepvars, depvar)
end

function MatrixDataSet{N <: Number}(X::Array{N, 2}, y::Array{N, 2}, names = nothing)
    @assert size(y, 1) == 1 || size(y, 2) == 1
    MatrixDataSet(X, y[:], names)
end

indepvars(m::MatrixDataSet) = m.indepvars
depvar(m::MatrixDataSet) = m.depvar
as_matrices(m::MatrixDataSet) = return m.X, m.y
num_cases(m::MatrixDataSet) = size(m.X, 2)
indepvars_as_row_major_matrix(m::MatrixDataSet) = m.X'
depvar_as_vector(m::MatrixDataSet) = m.y


text:
using Pkg
Pkg.activate("..")
using Documenter, ThreadPools

push!(LOAD_PATH,"../src/")

makedocs(sitename="ThreadPool Documentation")
cp("src/img", "build/img", force=true)


text:
using SimpleDirectMediaLayer.LibSDL2

SDL_GL_SetAttribute(SDL_GL_MULTISAMPLEBUFFERS, 16)
SDL_GL_SetAttribute(SDL_GL_MULTISAMPLESAMPLES, 16)

@assert SDL_Init(SDL_INIT_EVERYTHING) == 0 "error initializing SDL: $(unsafe_string(SDL_GetError()))"

win = SDL_CreateWindow("Game", SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED, 1000, 1000, SDL_WINDOW_SHOWN)
SDL_SetWindowResizable(win, SDL_TRUE)

renderer = SDL_CreateRenderer(win, -1, SDL_RENDERER_ACCELERATED | SDL_RENDERER_PRESENTVSYNC)

surface = IMG_Load(joinpath(@__DIR__, "..", "assets", "cat.png"))
tex = SDL_CreateTextureFromSurface(renderer, surface)
SDL_FreeSurface(surface)

w_ref, h_ref = Ref{Cint}(0), Ref{Cint}(0)
SDL_QueryTexture(tex, C_NULL, C_NULL, w_ref, h_ref)

try
    w, h = w_ref[], h_ref[]
    x = (1000 - w) ÷ 2
    y = (1000 - h) ÷ 2
    dest_ref = Ref(SDL_Rect(x, y, w, h))
    close = false
    speed = 300
    while !close
        event_ref = Ref{SDL_Event}()
        while Bool(SDL_PollEvent(event_ref))
            evt = event_ref[]
            evt_ty = evt.type
            if evt_ty == SDL_QUIT
                close = true
                break
            elseif evt_ty == SDL_KEYDOWN
                scan_code = evt.key.keysym.scancode
                if scan_code == SDL_SCANCODE_W || scan_code == SDL_SCANCODE_UP
                    y -= speed / 30
                    break
                elseif scan_code == SDL_SCANCODE_A || scan_code == SDL_SCANCODE_LEFT
                    x -= speed / 30
                    break
                elseif scan_code == SDL_SCANCODE_S || scan_code == SDL_SCANCODE_DOWN
                    y += speed / 30
                    break
                elseif scan_code == SDL_SCANCODE_D || scan_code == SDL_SCANCODE_RIGHT
                    x += speed / 30
                    break
                else
                    break
                end
            end
        end

        x + w > 1000 && (x = 1000 - w;)
        x < 0 && (x = 0;)
        y + h > 1000 && (y = 1000 - h;)
        y < 0 && (y = 0;)

        dest_ref[] = SDL_Rect(x, y, w, h)
        SDL_RenderClear(renderer)
        SDL_RenderCopy(renderer, tex, C_NULL, dest_ref)
        dest = dest_ref[]
        x, y, w, h = dest.x, dest.y, dest.w, dest.h
        SDL_RenderPresent(renderer)

        SDL_Delay(1000 ÷ 60)
    end
finally
    SDL_DestroyTexture(tex)
    SDL_DestroyRenderer(renderer)
    SDL_DestroyWindow(win)
    SDL_Quit()
end

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
module Workflow
    export Launcher;
    export get_ts_s;
    export get_units_by_kind;
    export get_units_by_bus;
    export get_bus;
    export sc_opf;
    export print_nz;

    using ..AmplTxt;
    using Dates: Date, DateTime;
    using JuMP;
    using Printf;
    using Parameters;
    U_NAME=1;
    U_SCENARIO=2;
    U_H=3;
    U_ECH=4;
    K_LIMITABLE="Limitable";
    K_IMPOSABLE="Imposable";

    @with_kw    mutable struct Launcher
        NO_IMPOSABLE::Bool;
        NO_LIMITABLE::Bool
        NO_LIMITATION::Bool;

        dirpath::String;
        # amplTxt description of the network
        ampltxt;
        # uncertainties, name-scenario-h-ech->value
        uncertainties::Dict{Tuple{String,String,DateTime,DateTime},Float64};
        # previsionnal planning
        previsions::Dict{Tuple{String,DateTime,DateTime},Float64};
        # name->minP-maxP-startCost-propCost
        units::Dict{String, Vector{Float64}};
        # gen->type-bus
        gen_type_bus::Dict{String, Vector{String}};
        # line-bus->ptdf
        ptdf::Dict{Tuple{String,String}, Float64};

        limits::Dict{String, Float64};

    end

    @with_kw    mutable struct ImposableModeler
        p_imposable = Dict{Tuple{String,DateTime,String},VariableRef}();
        p_is_imp = Dict{Tuple{String,DateTime},VariableRef}();
        p_is_imp_and_on = Dict{Tuple{String,DateTime},VariableRef}();
        p_imp = Dict{Tuple{String,DateTime},VariableRef}();
        p_start = Dict{Tuple{String,DateTime},VariableRef}();
        p_on = Dict{Tuple{String,DateTime},VariableRef}();
        c_imp_pos = Dict{Tuple{String,DateTime,String},VariableRef}();
        c_imp_neg = Dict{Tuple{String,DateTime,String},VariableRef}();
    end

    @with_kw mutable struct LimitableModeler
        p_lim = Dict{Tuple{String,DateTime},VariableRef}();

        is_limited = Dict{Tuple{String,DateTime,String},VariableRef}();
        p_enr = Dict{Tuple{String,DateTime,String},VariableRef}();
        is_limited_x_p_lim = Dict{Tuple{String,DateTime,String},VariableRef}();
        c_lim = Dict{Tuple{String,DateTime,String},VariableRef}();
    end

    @with_kw mutable struct ReserveModeler
        p_res_pos = Dict{Tuple{DateTime,String},VariableRef}();
        p_res_neg = Dict{Tuple{DateTime,String},VariableRef}();
    end

    function read_ptdf(dir_path::String)
        result = Dict{Tuple{String,String}, Float64}();
        open(joinpath(dir_path, "pscopf_ptdf.txt"), "r") do file
            for ln in eachline(file)
                # don't read commentted line
                if ln[1] != '#'
                    buffer = AmplTxt.split_with_space(ln);
                    result[buffer[1], buffer[2]] = parse(Float64, buffer[3]);
                end
            end
        end
        return result;
    end
    function read_limits(dir_path::String)
        result = Dict{String, Float64}();
        open(joinpath(dir_path, "pscopf_limits.txt"), "r") do file
            for ln in eachline(file)
                # don't read commentted line
                if ln[1] != '#'
                    buffer = AmplTxt.split_with_space(ln);
                    result[buffer[1]] = parse(Float64, buffer[2]);
                end
            end
        end
        return result;
    end

    function read_gen_type_bus(dir_path::String)
        result = Dict{String, Vector{String}}();
        open(joinpath(dir_path, "pscopf_gen_type_bus.txt"), "r") do file
            for ln in eachline(file)
                # don't read commentted line
                if ln[1] != '#'
                    buffer = AmplTxt.split_with_space(ln);
                    result[buffer[1]] = buffer[2:3];
                end
            end
        end
        return result
    end

    function read_units(dir_path::String)
        result = Dict{String, Vector{Float64}}();
        open(joinpath(dir_path, "pscopf_units.txt"), "r") do file
            for ln in eachline(file)
                # don't read commentted line
                if ln[1] != '#'
                    buffer = AmplTxt.split_with_space(ln);
                    result[buffer[1]] = (parse.(Float64, buffer[3:6]));
                end
            end
        end
        return result
    end

    function read_uncertainties(dir_path::String)
        result = Dict{Tuple{String,String,DateTime,DateTime},Float64}();
        open(joinpath(dir_path, "pscopf_uncertainties.txt"), "r") do file
            for ln in eachline(file)
                # don't read commentted line
                if ln[1] != '#'
                    # println(ln)
                    buffer = AmplTxt.split_with_space(ln);
                    name = buffer[1];
                    h = DateTime(buffer[2]);
                    ech = DateTime(buffer[3]);
                    scenario = buffer[4];
                    value = parse(Float64, buffer[5]);
                    result[name, scenario, h, ech] = value;
                end
            end
        end
        return result;
    end
    function read_previsions(dir_path::String)
        result = Dict{Tuple{String,DateTime,DateTime},Float64}();
        open(joinpath(dir_path, "pscopf_previsions.txt"), "r") do file
            for ln in eachline(file)
                # don't read commentted line
                if ln[1] != '#'
                    # println(ln)
                    buffer = AmplTxt.split_with_space(ln);
                    name = buffer[1];
                    h = DateTime(buffer[2]);
                    ech = DateTime(buffer[3]);
                    value = parse(Float64, buffer[4]);
                    result[name, h, ech] = value;
                end
            end
        end
        return result;
    end

    function Launcher(dir_path::String)
        return Launcher(false, false, false, dir_path, AmplTxt.read(dir_path),read_uncertainties(dir_path), read_previsions(dir_path), read_units(dir_path), read_gen_type_bus(dir_path), read_ptdf(dir_path), read_limits(dir_path))
    end

    function read_ampl_txt(launcher::Launcher, dir_path::String)
        launcher.ampltxt = AmplTxt.read(dir_path);
    end

    function add_uncertainties(launcher::Launcher, name::String, bus_name::String, ech::DateTime, ts::DateTime, value)
        launcher.uncertainties[name, bus_name, ech, ts] = value;
    end

    function get_bus(launcher::Launcher, names::Set{String})
        result = Set{String}();
        for name in names
            if !haskey(launcher.gen_type_bus, name)
                push!(result, name);
            end
        end
        return result;
    end

    function get_ech_ts_s_name(launcher::Launcher)
        ech_set = Set{DateTime}();
        ts_set = Set{DateTime}();
        s_set =  Set{String}();
        name_set = Set{String}();
        for uncertainty in launcher.uncertainties
            ts = uncertainty[1][U_H];
            ech = uncertainty[1][U_ECH];
            s = uncertainty[1][U_SCENARIO];
            name = uncertainty[1][U_NAME];
            push!(ech_set, ech);
            push!(ts_set, ts);
            push!(s_set, s);
            push!(name_set, name);
        end
        return sort(collect(ech_set)), sort(collect(ts_set)), sort(collect(s_set)), name_set;
    end

    function get_units_by_kind(launcher::Launcher)
        result = Dict{String,Dict{String,String}}();
        result[K_LIMITABLE] = Dict{String,String}();
        result[K_IMPOSABLE] = Dict{String,String}();
        for gen in launcher.gen_type_bus
            kind = gen[2][1];
            push!(result[kind],  gen[1] => gen[2][2]);
        end
        return result;
    end

    function get_units_by_bus(launcher::Launcher, buses::Set{String})
        result = Dict{String,Dict{String,Vector{String}}}();
        result[K_LIMITABLE] = Dict{String,Vector{String}}([bus => Vector{String}() for bus in buses]);
        result[K_IMPOSABLE] = Dict{String,Vector{String}}([bus => Vector{String}() for bus in buses]);
        for gen in launcher.gen_type_bus
            name = gen[1];
            kind = gen[2][1];
            bus  = gen[2][2];
            tmp = get(result[kind], bus, Vector{String}())
            push!(tmp, name);
            result[kind][bus] = tmp;
        end
        return result;
    end
    include("WorkflowModeler.jl");
    include("WorkflowWorstCase.jl");
end

max_stars_repo_name: azized/pscopf

max_issues_repo_name: azized/pscopf

max_forks_repo_name: azized/pscopf


text:
module Types

#abstract type AgentStates{T <: AbstractFloat} end

#abstract type Agent{S <: AgentStates} end

# abstract type Donor end # T::generated by agent, to be processed for sending. should change as "to be put to channel"
# abstract type Acceptor end # T::processed as what agent accepts. should be change as "taken from channel"

end


text:
################################
######    Polymorphism    ######
################################
# Expected number of polymorphism above frequency x from the standard diffusion theory
# f(x) = ∫(s) θs* (1/x(1-x)) * ( (ℯ^(Ns)*(1-ℯ^(-4Ns(1-x))) ) / (ℯ^(4Ns-1)))
# Convolote with the binomial to obtain the downsampled sfs
# E[P(x)] = ∑_x=x*→x=1 fB(x*)
############Neutral#############

"""

	DiscSFSNeutDown()

Expected rate of neutral allele frequency reduce by backgrou	nd selection. The spectrum depends on the number of individual []

```math
\\mathbb{E}[Ps_{(x)}] = \\sum{x^{*}=x}{x^{*}=1}f_{B}(x)
```

# Return:
 - `Array{Float64}`: expected rate of neutral alleles frequencies.
"""
function DiscSFSNeutDown(param::parameters,binom::SparseMatrixCSC{Float64,Int64})

	NN2 = convert(Int64,ceil(param.NN*param.B))
	# Allocating variables

	neutralSfs(i::Int64) = 1.0/(i)

	x = collect(0:NN2)
	solvedNeutralSfs = neutralSfs.(x)
	replace!(solvedNeutralSfs, Inf => 0.0)

	# subsetDict = get(param.bn,param.B,1)
	# subsetDict = binom
	out::Array{Float64,1} = param.B * (param.thetaMidNeutral) * 0.25 * (binom*solvedNeutralSfs)
	#=out::Array{Float64,1} = param.B .* (param.θᵣ) .*0.25 .* (binom*solvedNeutralSfs)=#
	return out
end

############Positive############
# Variable gamma in function changed to gammaValue to avoid problem with exported SpecialFunctions.gamma
"""

	DiscSFSSelPosDown(gammaValue,ppos)

Expected rate of positive selected allele frequency reduce by background selection. The spectrum depends on the number of individuals.

# Arguments
 - `gammaValue::Int64`: selection strength.
 - `ppos::Float64`: positive selected alleles probabilty.
# Return:
 - `Array{Float64}`: expected positive selected alleles frequencies.
"""
function DiscSFSSelPosDown(param::parameters,gammaValue::Int64,ppos::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if ppos == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else

		redPlus = phiReduction(param,gammaValue)

		# Solving sfs
		NN2 = convert(Int64,ceil(param.NN*param.B))
		xa1  = collect(0:NN2)
		xa2  = xa1/(NN2)

		# Solving float precision performance using exponential rule. Only one BigFloat estimation.
		gammaCorrected = gammaValue*param.B

		gammaExp1 = exp(gammaCorrected*2)
		gammaExp2 = exp(gammaCorrected*-2)

		positiveSfs(i::Float64,g1::Float64=gammaExp1,g2::Float64=gammaExp2,ppos::Float64=ppos) = Float64(ppos*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))

		# Original
		# ppos*0.5*(ℯ^(2*gammaCorrected)*(1-ℯ^(-2.0*gammaCorrected*(1.0-i)))/((ℯ^(2*gammaCorrected)-1.0)*i*(1.0-i)))

		# Allocating outputs
		solvedPositiveSfs::Array{Float64,1} = (1.0/(NN2)) * (positiveSfs.(xa2))
		replace!(solvedPositiveSfs, NaN => 0.0)

		# subsetDict = get(param.bn,param.B,1)
		# out               = (param.thetaMidNeutral)*redPlus*0.75*(subsetDict*solvedPositiveSfs)
		out::Array{Float64,1} = param.thetaMidNeutral * redPlus * 0.75 * (binom*solvedPositiveSfs)
		#=out::Array{Float64,1} = param.θᵣ .* redPlus .* 0.75 .* (binom*solvedPositiveSfs)=#

	end

	return out
end

function DiscSFSSelPosDownArb(param::parameters,gammaValue::Int64,ppos::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if ppos == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else
		redPlus = phiReduction(param,gammaValue)

		# Solving sfs
		NN2 = convert(Int64,ceil(param.NN*param.B))
		xa1  = collect(0:NN2)
		xa2  = xa1/(NN2)

		# Solving float precision performance using exponential rule. Only one BigFloat estimation.
		gammaCorrected = gammaValue*param.B
		gammaExp1::Quadmath.Float128 = exp(Quadmath.Float128(gammaCorrected*2))
		gammaExp2::Quadmath.Float128 = exp(Quadmath.Float128(gammaCorrected*-2))

		positiveSfs(i::Float64,g1::Quadmath.Float128=gammaExp1,g2::Quadmath.Float128=gammaExp2,ppos::Float64=ppos) = Float64(ppos*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))
		# Allocating outputs
		solvedPositiveSfs::Array{Float64,1} = (1.0/(NN2)) * (positiveSfs.(xa2))
		replace!(solvedPositiveSfs, NaN => 0.0)
		out::Array{Float64,1} = param.thetaMidNeutral * redPlus * 0.75 * (binom*solvedPositiveSfs)
		#=out::Array{Float64,1} = param.θᵣ .* redPlus .* 0.75 .* (binom*solvedPositiveSfs)=#

	end

	return out
end

# function DiscSFSSelPosDown(gammaValue::Int64,ppos::Float64)

# 	if ppos == 0.0
# 		out = zeros(Float64,adap.nn + 1)
# 	else

# 		redPlus = phiReduction(gammaValue)

# 		# Solving sfs
# 		NN2 = convert(Int64,ceil(adap.NN*adap.B))
# 		xa  = collect(0:NN2)
# 		xa  = xa/(NN2)

		# function positiveSfs(i,gammaCorrected=gammaValue*adap.B,ppos=ppos)
		# 	if i > 0 && i < 1.0
		# 		return ppos*0.5*(
		# 			ℯ^(2*gammaCorrected)*(1-ℯ^(-2.0*gammaCorrected*(1.0-i)))/((ℯ^(2*gammaCorrected)-1.0)*i*(1.0-i)))
		# 	end
		# 	return 0.0
		# end

# 		# Allocating outputs
# 		solvedNeutralSfs = Array{Float64}(undef,NN2 + 1)
# 		out              = Array{Float64}(undef,NN2 + 1)

# 		solvedPositiveSfs = (1.0/(NN2)) * (xa .|> positiveSfs)
# 		replace!(solvedPositiveSfs, NaN => 0.0)
# 		out               = (adap.thetaMidNeutral)*redPlus*0.75*(adap.bn[adap.B]*solvedPositiveSfs)
# 	end

# 	return view(out,2:lastindex(out)-1,:)
# end

######Slightly deleterious######
"""

	DiscSFSSelNegDown(param,ppos)

Expected rate of positive selected allele frequency reduce by background selection. Spectrum drawn on a gamma DFE. It depends on the number of individuals.

# Arguments
 - `ppos::Float64`: positive selected alleles probabilty.
# Return:
 - `Array{Float64}`: expected negative selected alleles frequencies.
"""
function DiscSFSSelNegDown(param::parameters,ppos::Float64,binom::SparseMatrixCSC{Float64,Int64})
	# subsetDict = get(param.bn,param.B,1)
	solvedNegative = DiscSFSSelNeg(param,ppos)
	out = param.B .* (param.thetaMidNeutral) .* 0.75 .* (binom*solvedNegative)
	#=out = param.B .* (param.θᵣ) .* 0.75 .* (binom*solvedNegative)=#

	return out
end

function DiscSFSSelNeg(param::parameters,ppos::Float64)

	beta     = param.be/(1.0*param.B)
	NN2      = convert(Int64, ceil(param.NN*param.B))
	xa       = collect(0:NN2)/NN2

	solveZ   = similar(xa)

	z(x::Float64,p::Float64=ppos) = (1.0-p)*(2.0^-param.al)*(beta^param.al)*(-zeta(param.al,x+beta/2.0) + zeta(param.al,(2+beta)/2.0))/((-1.0+x)*x)

	solveZ   = z.(xa)

	if (solveZ[1] == Inf || isnan(solveZ[1]))
		solveZ[1] = 0.0
	end
	if (solveZ[lastindex(solveZ)] == Inf || isnan(solveZ[lastindex(solveZ)]))
		solveZ[lastindex(solveZ)] = 0.0
	end

	return 1.0/(NN2+0.0).*solveZ
end

"""
	cumulativeSfs(sfsTemp)

Changing SFS considering all values above a frequency *x*. The original asymptotic-MK approach takes Pn(x) and Ps(x) as the number of polymorphic sites at frequency *x* rather than above *x*, but this approach scales poorly as sample size increases. We define the polymorphic spectrum as stated above since these quantities trivially have the same asymptote but are less affected by changing sample size.
"""
function cumulativeSfs(sfsTemp::Array,freqs::Bool=true)

	out      = Array{Float64}(undef, size(sfsTemp,1),size(sfsTemp,2))

	if freqs
		idx = 2
	else
		idx = 1
	end

	out[1,idx:end] = sum(sfsTemp[:,idx:end],dims=1)

	@simd for i in 2:(size(sfsTemp)[1])

		#=app = view(out,i-1,:) .- view(sfsTemp,i-1,:)=#
		app = out[i-1,idx:end] .- sfsTemp[i-1,idx:end]

		if sum(app) > 0.0
			out[i,idx:end] = app
		else
			out[i,idx:end] = zeros(length(app))
		end
	end

	if freqs
		out[:,1] = sfsTemp[:,1]
	end

	return out
end

"""
	reduceSfs(sfsTemp,bins)

Function to reduce the SFS into N bins.
"""
function reduceSfs(sfsTemp::Array,bins::Int64)

	freq  = collect(0:(size(sfsTemp,1)-1))/size(sfsTemp,1)
	h1    = StatsBase.fit(StatsBase.Histogram,freq,0:(1/(bins-1)):1)
	xmap1 = StatsBase.binindex.(Ref(h1), freq)

	tmp = hcat(sfsTemp,xmap1)
	out = Array{Float64}(undef,bins-1 ,size(sfsTemp,2))
	vIter =  convert(Array,unique(xmap1)')
	@simd for i = eachindex(vIter)
		@inbounds out[i,2:end] = sum(tmp[tmp[:,end] .== i,2:end-1],dims=1)
	end

	out[:,1] = collect(1:(bins-1)) ./ bins

	return (out)
end


text:
using Revise
using statmech_tm_solver
using TensorKit
using Zygote
using LinearAlgebra
using Optim

# MPO
#T = mpo_triangular_AF_ising()
T = mpo_triangular_AF_ising_alternative()
T2 = mpo_triangular_AF_ising()
T_adapter, _ = mpo_triangular_AF_ising_adapter()
d = 4

psi = TensorMap(rand, ComplexF64, ℂ^2*ℂ^d, ℂ^2)
_, psi = left_canonical(psi)

io = open("result_iTEBD.txt", "w")
close(io)
for chi in [2; 4; 8; 16; 32; 64]
    global psi
    convergence_flag = -10
    io = open("result_iTEBD.txt", "a+")
    while convergence_flag <= 0
        _, Tpsi = left_canonical(act(T, psi))
        _, psi = iTEBD_truncate(Tpsi, chi)
        #psi = variational_truncate(Tpsi, chi)

        F_value = free_energy(T, psi)
        F_value2 = free_energy(T2, act(T_adapter, psi))
        costfunc_value = nonherm_cost_func(T, psi)
        costfunc_value2 = nonherm_cost_func(T2, act(T_adapter, psi))
        fidelity_value = ln_fidelity(psi, Tpsi)

        if isapprox(costfunc_value, -fidelity_value; rtol=0.02, atol=1e-13)
            convergence_flag += 1
        end

        msg = "$chi $F_value $costfunc_value $fidelity_value $F_value2 $costfunc_value2"

        println(msg)
        println(io, msg)
    end
    close(io)
    quicksave("ckpt_iTEBD_chi$(chi)", psi)
end