text:
import LLVM

nfuncs = Ref(0)
struct Scope
    values::Dict{Symbol, LLVM.User}
    parent::Scope

    Scope() = new(Dict{Symbol, LLVM.AllocaInst}())
    Scope(parent::Scope) = new(Dict{Symbol, LLVM.User}(), parent)
end

# Look for variable in `scope` and if not found, look recursively in parent
function _get(scope::Scope, var::Symbol, default)
    v = get(scope.values, var, nothing)
    if v !== nothing
        return v
    else
        return isdefined(scope, :parent) ? _get(scope.parent, var, default) : default
    end
end

mutable struct CurrentScope
    scope::Scope
end

CurrentScope() = CurrentScope(Scope())
isglobalscope(c::CurrentScope) = !isdefined(c.scope, :parent)
open_scope!(currentscope::CurrentScope) = (currentscope.scope = Scope(currentscope.scope); currentscope)
Base.pop!(currentscope::CurrentScope) = (currentscope.scope = currentscope.scope.parent)
Base.setindex!(currentscope::CurrentScope, v::LLVM.User, var::Symbol) = currentscope.scope.values[var] = v
Base.get(currentscope::CurrentScope, var::Symbol, default) = _get(currentscope.scope, var, default)

mutable struct CodeGen
    ctx::LLVM.Context
    builder::LLVM.Builder
    current_scope::CurrentScope
    mod::LLVM.Module

    CodeGen(ctx::LLVM.Context) =
        new(
            ctx,
            LLVM.Builder(ctx),
            CurrentScope(),
            LLVM.Module("MSScheme"; ctx),
        )
end

current_scope(cg::CodeGen) = cg.current_scope
function new_scope(f, cg::CodeGen)
    open_scope!(current_scope(cg))
    f()
    pop!(current_scope(cg))
end

parse_expr(val::Float64, root = false) = NumberExprAST(val)
parse_expr(val::Symbol, root = false) = VarExprAST(val)
function parse_expr(expr::MyExpr, root = false)
    call = expr.call

    if call == :define
        root || error("Syntax Error: Nested definitions not allowed \n $expr")
        val = expr.args[1]
        if val isa MyExpr
            fname = val.call
            args = val.args
            all(isa.(args, Symbol)) || error("Syntax Error: Only identifiers can be used as parameters")
            body = expr.args[end]
            newexpr = parse_expr(body)
            if (newexpr isa IfExprAST) && (newexpr.elsee isa CallExprAST)
                if newexpr.elsee.call == fname
                    return BodyTailExprAST(HeadFunExprAST(fname, Symbol[args...]),newexpr)
                end
            end
            return BodyFunExprAST(HeadFunExprAST(fname, Symbol[args...]), newexpr)
        end
        if val isa Symbol
            varname = val
            length(expr.args) != 2 && error("Syntax Error: variable definition can only have one expression \n $expr")
            body = expr.args[end]
            newexpr =  parse_expr(body)
            return DefVarExprAST(varname, newexpr)
        end

    elseif any(call .== (:+, :*, :/, :-))  #reduce functions
        return ReduceExprAST(call, parse_expr.(expr.args))

    elseif call == :if
        length(expr.args) != 3 && error("Syntax Error: if block requires the (if cond then else) Syntax \n $expr")
        cond = expr.args[1]
        then = expr.args[2]
        els = expr.args[3]
        return IfExprAST(parse_expr(cond), parse_expr(then), parse_expr(els))

    elseif any(call .== (:<, :>, :(==), :!=, :(<=), :(>=)))
        length(expr.args) != 2 && error("Syntax Error: Boolean Expression must have only two arguments \n $expr")
        op = call
        lhs = expr.args[1]
        rhs = expr.args[2]
        return BoolOpExprAST(op, parse_expr(rhs), parse_expr(lhs))

    else
        return CallExprAST(call, parse_expr.(expr.args))
    end
end

function create_entry_block_allocation(cg::CodeGen, fn::LLVM.Function, varname::String)
    local alloc
    LLVM.Builder(cg.ctx) do builder
        # Set the builder at the start of the function
        entry_block = LLVM.entry(fn)
        if isempty(LLVM.instructions(entry_block))
            LLVM.position!(builder, entry_block)
        else
            LLVM.position!(builder, first(LLVM.instructions(entry_block)))
        end
        alloc = LLVM.alloca!(builder, LLVM.DoubleType(cg.ctx), varname)
    end
    return alloc
end

function codegen(cg::CodeGen, expr::NumberExprAST)
    return LLVM.ConstantFP(LLVM.DoubleType(cg.ctx), expr.val)
end

function codegen(cg::CodeGen, expr::VarExprAST)
    V = get(current_scope(cg), expr.name, nothing)
    String(expr.name)
    V == nothing && error("Undefined Variable $(expr.name)")
    return LLVM.load!(cg.builder, V, String(expr.name))
end

function codegen(cg::CodeGen, expr::DefVarExprAST)
    local initval
    initval = codegen(cg, expr.init)
    varname = expr.varname
    local V
    if isglobalscope(current_scope(cg))
        V = LLVM.GlobalVariable(cg.mod, LLVM.DoubleType(cg.ctx), String(varname))
        LLVM.initializer!(V, initval)
    else
        func = LLVM.parent(LLVM.position(cg.builder))
        V = create_entry_block_allocation(cg, func, String(varname))
        LLVM.store!(cg.builder, initval, V)
    end
    current_scope(cg)[varname] = V
    return initval
end

function codegen(cg::CodeGen, expr::BoolOpExprAST)
    if isglobalscope(current_scope(cg))
        fname = Symbol(nfuncs[])
        nfuncs[] += 1
        head = HeadFunExprAST(fname,Symbol[])
        return codegen(cg, BodyFunExprAST(head, expr))
    end
    L = codegen(cg, expr.lhs)
    R = codegen(cg, expr.rhs)
    if expr.op == :(==)
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealOEQ, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(cg.ctx), "booltmp")
    elseif expr.op == :(!=)
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealONE, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(cg.ctx), "booltmp")
    elseif expr.op == :(<)
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealOLT, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(cg.ctx), "booltmp")
    elseif expr.op == :(>)
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealOGT, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(cg.ctx), "booltmp")
    elseif expr.op == :(<=)
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealOLE, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(cg.ctx), "booltmp")
    elseif expr.op == :(>=)
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealOGE, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(cg.ctx), "booltmp")
    else
        error("Operação booleana invalida $(expr.op)")
    end
end

function codegen(cg::CodeGen, expr::ReduceExprAST)
    if isglobalscope(current_scope(cg))
        fname = Symbol(nfuncs[])
        nfuncs[] += 1
        head = HeadFunExprAST(fname,Symbol[])
        return codegen(cg, BodyFunExprAST(head, expr))
    end
    op = expr.fun
    args = LLVM.Value[]
    for v in expr.args
        push!(args, codegen(cg, v))
    end
    tot = args[1]
    for arg in args[2:end]
        if op == :+
            tot = LLVM.fadd!(cg.builder, tot, arg, "addtmp")
        elseif op == :-
            tot = LLVM.fsub!(cg.builder, tot, arg, "addtmp")
        elseif op == :*
            tot = LLVM.fmul!(cg.builder, tot, arg, "addtmp")
        elseif op == :/
            tot = LLVM.fdiv!(cg.builder, tot, arg, "addtmp")
        else
            error("Não devia chegar aqui")
        end
    end
    return tot
end

function codegen(cg::CodeGen, expr::CallExprAST)
    if isglobalscope(current_scope(cg))
        fname = Symbol(nfuncs[])
        nfuncs[] += 1
        head = HeadFunExprAST(fname,Symbol[])
        return codegen(cg, BodyFunExprAST(head, expr))
    end
    if !haskey(LLVM.functions(cg.mod), String(expr.call))
        error("Undefined Function $(expr.call)")
    end
    func =  LLVM.functions(cg.mod)[String(expr.call)]

    if length(LLVM.parameters(func)) != length(expr.args)
        error("Incorrect number of arguments $fname requires $(length(expr.args)) arguments")
    end

    args = LLVM.Value[]
    for v in expr.args
        push!(args, codegen(cg, v))
    end

    return LLVM.call!(cg.builder, func, args, "calltmp")
end

function codegen(cg::CodeGen, expr::HeadFunExprAST)
    if haskey(LLVM.functions(cg.mod), String(expr.name))
            error("Can't redefine function")
    end
    args = [LLVM.DoubleType(cg.ctx) for i in 1:length(expr.args)]
    func_type = LLVM.FunctionType(LLVM.DoubleType(cg.ctx), args)
    func = LLVM.Function(cg.mod, String(expr.name), func_type)
    LLVM.linkage!(func, LLVM.API.LLVMExternalLinkage)

    for (i, param) in enumerate(LLVM.parameters(func))
        LLVM.name!(param, String(expr.args[i]))
    end
    return func
end

function codegen(cg::CodeGen, expr::BodyFunExprAST)
    # create new function...
    the_function = codegen(cg, expr.head)

    entry = LLVM.BasicBlock(the_function, "entry"; cg.ctx)
    LLVM.position!(cg.builder, entry)

    new_scope(cg) do
        for (i, param) in enumerate(LLVM.parameters(the_function))
            argname = expr.head.args[i]
            alloc = create_entry_block_allocation(cg, the_function, String(argname))
            LLVM.store!(cg.builder, param, alloc)
            current_scope(cg)[argname] = alloc
        end

        body = codegen(cg, expr.body)
        LLVM.ret!(cg.builder, body)
        LLVM.verify(the_function)
    end
    return the_function
end

function codegen(cg::CodeGen, expr::BodyTailExprAST)
    # create new function...
    the_function = codegen(cg, expr.head)

    entry = LLVM.BasicBlock(the_function, "entry"; cg.ctx)
    exit = LLVM.BasicBlock(the_function, "exit"; cg.ctx)
    tail = LLVM.BasicBlock(the_function, "tail"; cg.ctx)
    merge = LLVM.BasicBlock(the_function, "ifcont"; cg.ctx)
    body = expr.body
    local phi
    LLVM.position!(cg.builder, entry)
    new_scope(cg) do
        for (i, param) in enumerate(LLVM.parameters(the_function))
            argname = expr.head.args[i]
            alloc = create_entry_block_allocation(cg, the_function, String(argname))
            LLVM.store!(cg.builder, param, alloc)
            current_scope(cg)[argname] = alloc
        end

        #cond
        cond = codegen(cg, body.cond)
        zero = LLVM.ConstantFP(LLVM.DoubleType(cg.ctx), 0.0)
        condv = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealONE, cond, zero, "ifcond")
        LLVM.br!(cg.builder, condv, exit, tail)
        #exit
        LLVM.position!(cg.builder, exit)
        exitcg = codegen(cg, body.then)
        LLVM.br!(cg.builder, merge)
        exit_block = position(cg.builder)

        # tail
        LLVM.position!(cg.builder, tail)
        tailcg =codegen(cg, body.elsee)
        LLVM.tailcall!(tailcg, true)
        LLVM.ret!(cg.builder, tailcg)
        # LLVM.br!(cg.builder, merge)
        # tail_block = position(cg.builder)

        #merge
        LLVM.position!(cg.builder, merge)
        phi = LLVM.phi!(cg.builder, LLVM.DoubleType(cg.ctx), "iftmp")
        append!(LLVM.incoming(phi), [(exitcg, exit_block)])
        LLVM.ret!(cg.builder, phi)
        # LLVM.verify(the_function)
    end
    return the_function
end

function codegen(cg::CodeGen, expr::IfExprAST)
    if isglobalscope(current_scope(cg))
        fname = Symbol(nfuncs[])
        nfuncs[] += 1
        head = HeadFunExprAST(fname,Symbol[])
        return codegen(cg, BodyFunExprAST(head, expr))
    end

    func = LLVM.parent(LLVM.position(cg.builder))
    then = LLVM.BasicBlock(func, "then"; cg.ctx)
    elsee = LLVM.BasicBlock(func, "else"; cg.ctx)
    merge = LLVM.BasicBlock(func, "ifcont"; cg.ctx)

    local phi
    new_scope(cg) do
        # if
        cond = codegen(cg, expr.cond)
        zero = LLVM.ConstantFP(LLVM.DoubleType(cg.ctx), 0.0)
        condv = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealONE, cond, zero, "ifcond")
        LLVM.br!(cg.builder, condv, then, elsee)
        # then
        LLVM.position!(cg.builder, then)
        thencg = codegen(cg, expr.then)
        LLVM.br!(cg.builder, merge)
        then_block = position(cg.builder)

        # else
        LLVM.position!(cg.builder, elsee)
        elsecg = codegen(cg, expr.elsee)
        LLVM.br!(cg.builder, merge)
        else_block = position(cg.builder)

        # merge
        LLVM.position!(cg.builder, merge)
        phi = LLVM.phi!(cg.builder, LLVM.DoubleType(cg.ctx), "iftmp")
        append!(LLVM.incoming(phi), [(thencg, then_block), (elsecg, else_block)])
    end

    return phi
end

function generate_IR(code::String, ctx::LLVM.Context)
    nfuncs[] = 0
    exprs =global_machine(code)
    cg = CodeGen(ctx)
    for expr in exprs
        codegen(cg, parse_expr(expr, true))
    end
    LLVM.verify(cg.mod)
    LLVM.dispose(cg.builder)
    return cg.mod
end

function run_module(mod::LLVM.Module, entry::String, optimize = false)
    res_jl = 0.0
    if optimize
        mod = optimize!(mod)
    end
    LLVM.JIT(mod) do engine
        if !haskey(LLVM.functions(engine), entry)
            error("Did not find entry function '$entry' in module")
        end
        f = LLVM.functions(engine)[entry]
        res = LLVM.run(engine, f)
        res_jl = convert(Float64, res, LLVM.DoubleType(LLVM.context(mod)))
        LLVM.dispose(res)
    end
    return res_jl
end

function write_objectfile(mod::LLVM.Module, path::String)
    host_triple = Sys.MACHINE # LLVM.triple() might be wrong (see LLVM.jl#108)
    host_t = LLVM.Target(triple=host_triple)
    LLVM.TargetMachine(host_t, host_triple) do tm
        LLVM.emit(tm, mod, LLVM.API.LLVMObjectFile, path)
    end
end

function optimize!(mod::LLVM.Module)
    LLVM.ModulePassManager() do pass_manager
        LLVM.instruction_combining!(pass_manager)
        LLVM.reassociate!(pass_manager)
        LLVM.gvn!(pass_manager)
        LLVM.cfgsimplification!(pass_manager)
        LLVM.promote_memory_to_register!(pass_manager)
        LLVM.run!(pass_manager, mod)
    end
    return mod
end


text:
# This file constains the Clock tools for time synchronization of DsSimulator.

import Base: iterate, take!, length

Generator(t0, dt, tf) =
    Channel{promote_type(typeof(t0),typeof(dt),typeof(tf))}(channel -> foreach(t -> put!(channel, t), t0:dt:tf))

"""
    Clock(t::Real, dt::Real, tf::Real)

Constructs a `Clock` with starting time `t`, final time `tf` and sampling inteval `dt`. When iterated, the `Clock` returns its current time.

!!! warning
    When constructed, `Clock` is not running. To take clock ticks from `Clock`, the `Clock` must be setted. See [`take!(clk::Clock)`](@ref) and [`set!`](@ref)
"""
mutable struct Clock{T, CB}
    t::T
    ti::T
    dt::T
    tf::T
    generator::Channel{T}
    paused::Bool
    callbacks::CB
    name::Symbol
    uuid::UUID
    Clock{T}(ti::T, dt::T, tf::T; callbacks::CB=nothing, name=Symbol()) where {T, CB} =
        new{T, CB}(ti, ti, dt, tf, Channel{T}(0), false, callbacks, name, uuid4())
end
Clock(t::T, dt::T, tf::T; kwargs...) where T = Clock{T}(t, dt, tf; kwargs...)
Clock(t=0., dt=0.01, tf=1.; kwargs...) = Clock(promote(t, dt, tf)...; kwargs...)

show(io::IO, clk::Clock) = print(io,
    "Clock(t:$(clk.t), dt:$(clk.dt), tf:$(clk.tf), paused:$(clk.paused), isrunning:$(isrunning(clk)))")

##### Reading from clock
"""
    take!(clk::Clock)

Takes a values from `clk`.

# Example
```jldoctest
julia> clk = Clock(0., 0.1, 0.5)
Clock(t:0.0, dt:0.1, tf:0.5, paused:false, isrunning:false)

julia> set!(clk)
Clock(t:0.0, dt:0.1, tf:0.5, paused:false, isrunning:true)

julia> for i in 0 : 5
       @show take!(clk)
       end
take!(clk) = 0.0
take!(clk) = 0.1
take!(clk) = 0.2
take!(clk) = 0.3
take!(clk) = 0.4
take!(clk) = 0.5
```
"""
function take!(clk::Clock)
    if ispaused(clk)
        @warn "Clock is paused."
        return clk.t
    end
    if isoutoftime(clk)
        @warn "Clock is out of time."
        return clk.t
    end
    ##
    ##  NOTE: If this code block is uncommented, a bug occurs
    ##  The same clock times is sent multiple time.
    ##
    # if !isrunning(clk)
    #     @warn "Clock is not running."
    #     return clk.t
    # end
    clk.t = take!(clk.generator)
    applycallbacks(clk)
    clk.t
end

##### Clock state check
"""
    isrunning(clk::Clock)

Returns `true` if `clk` if `clk` is running.
"""
isrunning(clk::Clock) = isready(clk.generator)

"""
    ispaused(clk::Clock)

Returns `true` if `clk` is paused. When paused, the currnent time of `clk` is not advanced. See also [`pause!(clk::Clock)`](@ref)
"""
ispaused(clk::Clock) = clk.paused

"""
    isoutoftime(clk::Clock)

Returns `true` if `clk` is out of time, i.e., the current time of `clk` exceeds its final time.
"""
isoutoftime(clk::Clock) = clk.t >= clk.tf

##### Controlling clock.
"""
    set(clk::Clock, t::Real, dt::Real, tf::Real)

Sets `clk` for current clock time `t`, sampling time `dt` and final time `tf`. After the set,  it is possible to take clock tick from `clk`. See also [`take!(clk::Clock)`](@ref)

# Example
```jldoctest
julia> clk = Clock(0., 0.1, 0.5)
Clock(t:0.0, dt:0.1, tf:0.5, paused:false, isrunning:false)

julia> set!(clk)
Clock(t:0.0, dt:0.1, tf:0.5, paused:false, isrunning:true)

julia> take!(clk)
0.0
```
"""
function set!(clk::Clock, t::Real=clk.ti, dt::Real=clk.dt, tf::Real=clk.tf)
    # set!(clk, Generator(t, dt, tf))
    clk.generator = Generator(t, dt, tf)
    clk.t = t
    clk.dt = dt
    clk.tf = tf
    clk.paused = false
    clk
end
# function set!(clk::Clock, generator::Channel=Generator(clk.ti, clk.dt, clk.tf))
#     clk.generator = generator
#     clk.paused = false
#     clk
# end

"""
    stop!(clk::Clock)

Unsets `clk`. After the stpp, it is possible to take clock ticks from `clk`. See also [`take!(clk::Clock)`](@ref)
"""
function stop!(clk::Clock)
    set!(clk, Channel{typeof(clk.t)}(0))
    clk
end

"""
    pause!(clk::Clock)

Pauses `clk`. When paused, the current time of `clk` does not advance.

# Example
```julia
julia> clk = Clock(0., 0.1, 0.5);

julia> set!(clk);

julia> for i = 1 : 5
              i > 3 && pause!(clk)
              @show take!(clk)
              end
take!(clk) = 0.0
take!(clk) = 0.1
take!(clk) = 0.2
┌ Warning: Clock is paused.
└ @ Causal ~/.julia/dev/Causal/src/components/sources/clock.jl:61
take!(clk) = 0.2
┌ Warning: Clock is paused.
└ @ Causal ~/.julia/dev/Causal/src/components/sources/clock.jl:61
take!(clk) = 0.2
```
"""
pause!(clk::Clock) = (clk.paused = true; clk)

##### Iterating clock.
"""
    iterate(clk::Clock[, t=clk.t)

Iterationk interface for `clk`. `clk` can be iterated in a loop.

# Example
```jldoctest
julia> clk = Clock(0., 0.1, 0.3);

julia> set!(clk)
Clock(t:0.0, dt:0.1, tf:0.3, paused:false, isrunning:true)

julia> for t in clk
       @show t
       end
t = 0.0
t = 0.1
t = 0.2
t = 0.3
```
"""
iterate(clk::Clock, t=clk.t) = isready(clk.generator) ? (take!(clk), clk.t) : nothing

##### ProgressMeter interface.
### This `length` method is implemented for [ProgressMeter](https://github.com/timholy/ProgressMeter.jl)
length(clk::Clock) = length(clk.t:clk.dt:clk.tf)

ext: jl


text:
using Documenter, Cassette

makedocs(modules=[Cassette],
         doctest = false,
         format = :html,
         sitename = "Cassette",
         pages = ["Introduction" => "index.md",
                  "Disclaimers" => "disclaimers.md",
                  "Why Cassette?" => "whycassette.md",
                  "The Overdubbing Mechanism" => "overdub.md",
                  "Contextual Dispatch" => "contextualdispatch.md",
                  "Contextual Compiler Pass Injection" => "contextualpass.md",
                  "Contextual Tagging of Values" => "contextualtagging.md",
                  "Cassette API Documentation" => "api.md",
                  "Related Work" => "relatedwork.md"])

deploydocs(repo = "github.com/jrevels/Cassette.jl.git",
           osname = "linux",
           julia = "0.7",
           target = "build",
           deps = nothing,
           make = nothing)


text:
typealias WorkInt Union{Int64, Int32}

ClCl{R<:Real}(x::Rvl{ClCl,R}) = x
ClOp{R<:Real}(x::Rvl{ClCl,R}) = Rvl{ClOp,R}(x.lo,nextFloat(x.hi))
OpCl{R<:Real}(x::Rvl{ClCl,R}) = Rvl{OpCl,R}(prevFloat(x.lo),x.hi)
OpOp{R<:Real}(x::Rvl{ClCl,R}) = Rvl{OpOp,R}(prevFloat(x.lo),nextFloat(x.hi))

ClCl{R<:Real}(x::Rvl{ClOp,R}) = Rvl{ClCl,R}(x.lo,x.hi)
ClOp{R<:Real}(x::Rvl{ClOp,R}) = x
OpCl{R<:Real}(x::Rvl{ClOp,R}) = Rvl{OpCl,R}(prevFloat(x.lo),x.hi)
OpOp{R<:Real}(x::Rvl{ClOp,R}) = Rvl{OpOp,R}(prevFloat(x.lo),x.hi)

ClCl{R<:Real}(x::Rvl{OpCl,R}) = Rvl{ClCl,R}(x.lo,x.hi)
ClOp{R<:Real}(x::Rvl{OpCl,R}) = Rvl{ClOp,R}(x.lo,nextFloat(x.hi))
OpCl{R<:Real}(x::Rvl{OpCl,R}) = x
OpOp{R<:Real}(x::Rvl{OpCl,R}) = Rvl{OpOp,R}(x.lo,nextFloat(x.hi))

ClCl{R<:Real}(x::Rvl{OpOp,R}) = Rvl{ClCl,R}(x.lo,x.hi)
ClOp{R<:Real}(x::Rvl{OpOp,R}) = Rvl{ClOp,R}(x.lo,x.hi)
OpCl{R<:Real}(x::Rvl{OpOp,R}) = Rvl{OpCl,R}(x.lo,x.hi)
OpOp{R<:Real}(x::Rvl{OpOp,R}) = x

ClCl!{G<:Grasp,R<:Real}(x::Rvl{G,R}) = ClCl(x.lo,x.hi)
ClOp!{G<:Grasp,R<:Real}(x::Rvl{G,R}) = ClOp(x.lo,x.hi)
OpCl!{G<:Grasp,R<:Real}(x::Rvl{G,R}) = OpCl(x.lo,x.hi)
OpOp!{G<:Grasp,R<:Real}(x::Rvl{G,R}) = OpOp(x.lo,x.hi)


text:
# Standard rectlin layer: f(A,x,b)=rectlin(A*x+b)
# Note that the bias (a column vector) gets expanded here to match the dimension of A*x

FrectlinAXplusBias(A,X,b)=begin; a=A*X+b*ones(1,size(X,2)); return (max(a,0.),(a.>0)); end

function FrectlinAXplusBias_inplace(value,aux,A,X,b)
    a=A*X+b*ones(1,size(X,2))
    copy!(value,max(a,0.0))
    copy!(aux,(a.>0))
end

function DrectlinAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X,b)
    if derivativeIDX==1
        axpy!(1.0,(grad_c.*faux_c)*X',grad_n)
    elseif derivativeIDX==2
        axpy!(1.0,A'*(grad_c.*faux_c),grad_n)
    elseif derivativeIDX==3
        axpy!(1.0,sum(grad_c.*faux_c,2),grad_n)
    end
end

if PROC=="GPU"
    function FrectlinAXplusBias(A::CudaArray,X::CudaArray,b::CudaArray)
        return(rectlin(FAXplusBias(A,X,b)[1]),FAXplusBias(A,X,b)[1])
        # TODO: FIX THIS UNNECESSARY DOUBLE COMPUTATION
    end

    function FrectlinAXplusBias_inplace(value,aux,A::CudaArray,X::CudaArray,b::CudaArray)
        FAXplusBias_inplace(value,aux,A,X,b)
        aux=CudaArray(Float64,size(value))
        copy!(aux,value)
        rectlin!(value,value)
    end

    function DrectlinAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray,b::CudaArray)
        tmp=CudaArray(Float64,size(grad_c)); fill!(tmp,0.0)
        A_emult_Bg0!(grad_c,faux_c,tmp)
        if derivativeIDX==1
            gemm!('N','T',1.0,tmp,X,1.0,grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',1.0,A,tmp,1.0,grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float64,(size(X,2),1)); fill!(ons,1.0)
            gemm!('N','N',1.0,tmp,ons,1.0,grad_n)
            free(ons)
        end
        free(tmp)
    end

end

Derivative[FrectlinAXplusBias]=DrectlinAXplusBias
Inplace[FrectlinAXplusBias]=FrectlinAXplusBias_inplace

rectlinAXplusBias(A,X,b)=ADnode(FrectlinAXplusBias,[A X b])
export rectlinAXplusBias

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
export postprocess
function postprocess(Events::Array{Event,1})
    set_context("poster")
    set_context("poster")
    close("all")
    pathlengths = zeros(0)
    distance = 0.0

    distancetofirst = true
    for i=1:length(Events)
        distance += Events[i].distance
        if Events[i].collision
            if distancetofirst
                distancetofirst = false # only store distances between collision
            else
                push!(pathlengths, distance)
            end
            distance = 0.0
        end
    end
    return postprocess(pathlengths)
end

function postprocess(pathlengths::Array{Float64,1},latticeflag = false)
    set_context("poster")
    set_context("poster")
    close("all")

    sort!(pathlengths)

    npl = length(pathlengths)

    fig, ax = plt.subplots(1,1,figsize=(40,20))

    logpathlengths = log10.(pathlengths)

    h = fit(Histogram, pathlengths, nbins = 200)

    centers = (h.edges[1][1:end-1] .+ h.edges[1][2:end])/2
    weights = h.weights
    weights = weights/sum(weights)/(centers[2]-centers[1])
    ax.plot(centers,weights,lw = 5,"-o",color=tuple([0,150,130]/255...))
    ax.set_xlim([0,1.05*maximum(pathlengths)])
    plt.yscale("log")

    d = fit(Exponential, pathlengths)
	if !latticeflag
		ax.plot(centers,pdf.(d,centers),lw=5,"--",color=tuple([163,16,124]/255...))
	end

    m(t, p) = p[1] .* t.^(-2.0)#p[2]
    p0 = [3.0]#, -2.0]
    lsqfit = curve_fit(m, centers, weights, p0)

    m(t) = lsqfit.param[1].*t.(-2.0)#^lsqfit.param[2]
    if latticeflag
		ax.plot(centers,m.(centers),lw=2,"--")
	end
	ax.set_title("Path-length distribution for n=$npl samples.", pad = 20)
    ax.set_xlabel("Distance s")
    ax.set_ylabel("Rel. frequency")
   # ax.set_ylim([1e-16,1])
    theta = params(d)[1]
    lambda = round(1.0/theta,digits = 3)
    sigma = round(lsqfit.param[1],digits = 3)
    #alpha = round(lsqfit.param[2],digits = 3)
	#if !latticeflag
#		ax.legend(["Data","Best fit, \$ \\lambda \\cdot \\exp( -\\lambda s)\$, \$\\lambda=\$ $lambda"])#
#	else
#		ax.legend(["Best fit, \$ \\sigma \\cdot s^\\alpha \$, \$\\sigma=\$ $sigma, \$\\alpha=\$ $alpha"])
#	end
    plt.tight_layout()
    plt.savefig("Histogram.png",dpi=200)

    return pathlengths
end


text:
# Print the floating-point number with precision high enough that rounding never occurs
function float_string(val::AbstractFloat)
    # Note: This function could simply be `@sprintf("%.325f", val)` once the following
    # issue is solved: https://github.com/JuliaLang/julia/issues/22137
    dp = 325  # Number of decimal places to print `nextfloat(0.0)` without rounding
    int = trunc(BigInt, val)
    if abs(int) > 0
        dp -= ndigits(int)
    end
    @eval @sprintf($("%.$(dp)f"), $val)
end

# FixedDecimal methods which perform an alternative method of trunc, floor, and ceil which
# primarily use the string representation of the floating-point. These `*_alt` methods exist
# to ensure the accuracy of the packages mathematical methods for trunc, floor, and ceil.

function integer_alt(::Type{T}, dp::Integer, val::AbstractFloat) where {T<:Integer}
    # Note: Use a precision larger than the value can represent so that `sprintf` doesn't
    # perform any rounding.
    str = float_string(val)
    sign = T(first(str) == '-' ? -1 : 1)
    decimal = findfirst(equalto('.'), str)
    int_start = sign < 0 ? 2 : 1
    int_end = decimal + dp
    v = parse(T, str[int_start:(decimal - 1)] * str[(decimal + 1):int_end])
    r = T(any(d -> d != '0', str[(int_end + 1):end]))  # Remaining digits
    return (sign, v, r)
end

function trunc_alt(::Type{FD{T,f}}, val::AbstractFloat) where {T<:Integer, f}
    s, v, r = integer_alt(T, f, val)
    reinterpret(FD{T,f}, copysign(v, s))
end

function floor_alt(::Type{FD{T,f}}, val::AbstractFloat) where {T<:Integer, f}
    s, v, r = integer_alt(T, f, val)
    reinterpret(FD{T,f}, copysign(v + (s < 0 ? r : zero(T)), s))
end

function ceil_alt(::Type{FD{T,f}}, val::AbstractFloat) where {T<:Integer, f}
    s, v, r = integer_alt(T, f, val)
    reinterpret(FD{T,f}, copysign(v + (s > 0 ? r : zero(T)), s))
end


text:
"""
    fresnel2(Uin, d, λ, lx, ly)

calculate the propagation light field based on the Fresnel diffraction with double Fourier transform.

## Arguments

- `Uin::AbstractArray{<:Number,2}`: Complex array of input complex amplitude.
- `d::Real`: Distance to propagate in metres.
- `λ::Real`: Wavelength of light to propagate.
- `lx::Real`: The physical size of the input data along the x-axis.
- `ly::Real`: The physical size of the input data along the y-axis.

## Returns

- `::Array{<:Number,2}`: Complex amplitude data after propagation.
"""
function fresnel2(Uin::AbstractArray{<:Number,2}, d::Real, λ::Real, lx::Real, ly::Real)
    (n, m) = size(Uin)
    ua = (-n/2:n/2-1)/lx*λ
    va = (-m/2:m/2-1)/ly*λ
    ifft(fft(Uin).*ifftshift([exp(2im*pi*d/λ)*exp(-im*pi*d/λ*(u^2+v^2)) for u in ua, v in va]))
end

"""
    fresnel2(Uin::MonoLightField2D, d::Unitful.Length) -> MonoLightField2D

calculate the propagation light field for [`MonoLightField2D`](@ref) based on the Fresnel diffraction with double Fourier transform.
"""
function fresnel2(Uin::MonoLightField2D, d::Unitful.Length)
    Uout = fresnel2(Uin.data, auval(d), auval(Uin.wavelength), auval(Uin.size[1]), auval(Uin.size[2]))
    MonoLightField2D(Uin, data=Uout)
end

"""
    fresnel2(d)

represent a [`fresnel2`](@ref) propagation functor with the specified distance.
"""
fresnel2(d) = x->fresnel2(x, d)

"""
    fresnel1(Uin, d, λ, lx, ly)

calculate the propagation light field based on the Fresnel diffraction with single Fourier transform.

## Arguments

- `Uin::AbstractArray{<:Number,2}`: Complex array of input complex amplitude.
- `d::Real`: Distance to propagate in metres.
- `λ::Real`: Wavelength of light to propagate.
- `lx::Real`: The physical size of the input data along the x-axis.
- `ly::Real`: The physical size of the input data along the y-axis.

## Returns

    (Uout, (lxo, lyo))

- `Uout::Array{<:Number,2}`: Complex amplitude data after propagation.
- `lxo::Real`: The physical size of the diffraction data along the x-axis.
- `lyo::Real`: The physical size of the diffraction data along the y-axis.
"""
function fresnel1(Uin::AbstractArray{<:Number,2}, d::Real, λ::Real, lx::Real, ly::Real)
    (n, m) = size(Uin)
    dx = lx/n
    dy = ly/m
    lxo = λ*d/dx
    lyo = λ*d/dy
    xa = (-n/2:n/2-1)/lx*λ
    ya = (-m/2:m/2-1)/ly*λ
    amp = fftshift(fft(Uin.*[exp(im*pi/d/λ*(x^2+y^2)) for x in xa, y in ya]))
    xa *= lx/lxo
    ya *= ly/lyo
    Uout = amp.*[exp(im*pi/d/λ*(x^2+y^2)) for x in xa, y in ya].*exp(2im*pi*d/λ).*dx.*dy./(im*λ*d)
    (Uout, (lxo, lyo))
end

"""
    fresnel1(Uin::MonoLightField2D, d::Unitful.Length) -> MonoLightField2D

calculate the propagation light field for [`MonoLightField2D`](@ref) based on the Fresnel diffraction with single Fourier transform.
"""
function fresnel1(Uin::MonoLightField2D, d::Unitful.Length)
    (Uout, Lout) = fresnel1(Uin.data, auval(d), auval(Uin.wavelength), auval(Uin.size[1]), auval(Uin.size[2]))
    MonoLightField2D(Uin, data=Uout, size=Lout)
end

"""
    fresnel1(d)

represent a [`fresnel1`](@ref) propagation functor with the specified distance.
"""
fresnel1(d) = x->fresnel1(x, d)


text:
##
#   File:    qo1.jl
#
#   Purpose: Demonstrate how to solve a quadratic
#            optimization problem using the MOSEK Python API.
##

using Mosek

printstream(msg::String) = print(msg)

# Define a stream printer to grab output from MOSEK

bkc   = [ MSK_BK_LO ]
blc   = [ 1.0 ]
buc   = [ Inf ]

bkx   = [ MSK_BK_LO, MSK_BK_LO, MSK_BK_LO ]
blx   = [ 0.0,  0.0, 0.0 ]
bux   = [ Inf,  Inf, Inf ]
numvar = length(bkx)
numcon = length(bkc)

c     = [ 0.0, -1.0, 0.0 ]
A     = sparse( [ 1, 1, 1 ],
                [ 1, 2, 3 ],
                [ 1.0, 1.0, 1.0 ],
                numcon, numvar )

maketask() do task
    putstreamfunc(task,MSK_STREAM_LOG,printstream)

    # Append 'numcon' empty constraints.
    # The constraints will initially have no bounds.
    appendcons(task,numcon)

    # Append 'numvar' variables.
    # The variables will initially be fixed at zero (x=0).
    appendvars(task,numvar)

    # Set the linear term c_j in the objective.
    putclist(task,[1:numvar;],c)

    # Set the bounds on variable j
    # blx[j] <= x_j <= bux[j]
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)

    putacolslice(task,1,numvar+1,
                 A.colptr[1:numvar],A.colptr[2:numvar+1],
                 A.rowval,A.nzval)

    # Set up and input quadratic objective
    qsubi = [ 1,   2,    3,   3   ]
    qsubj = [ 1,   2,    1,   3   ]
    qval  = [ 2.0, 0.2, -1.0, 2.0 ]

    putqobj(task,qsubi,qsubj,qval)

    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)

    # Optimize
    optimize(task)
    # Print a summary containing information
    # about the solution for debugging purposes
    solutionsummary(task,MSK_STREAM_MSG)

    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    if solsta in     [ MSK_SOL_STA_OPTIMAL,
                      MSK_SOL_STA_NEAR_OPTIMAL ]
        xx = getxx(task,MSK_SOL_ITR)
        println("Optimal solution:")
        print(xx)
    elseif solsta in [ MSK_SOL_STA_DUAL_INFEAS_CER,
                      MSK_SOL_STA_PRIM_INFEAS_CER,
                      MSK_SOL_STA_NEAR_DUAL_INFEAS_CER,
                      MSK_SOL_STA_NEAR_PRIM_INFEAS_CER ]
        println("Primal or dual infeasibility certificate found.\n")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status")
    else
        @printf("Other solution status (%d)\n",solsta)
    end

end


text:
@with_kw struct Analysis{T}
    data::SelectedData{T}
    compute_error = nothing
    x = nothing
    y = nothing
    z = nothing
    xfunc = nothing
    yfunc = nothing
    zfunc = nothing
    axis_type = nothing
    smoother = nothing
    package = nothing
    plot = nothing
    plot_kwargs = []
end

struct StatPlotsRecipe; end
struct GroupedError; end

function analysistype(a)
    a.package != nothing && return a.package
    a.plot in [boxplot, violin, histogram2d, marginalhist] && return StatPlotsRecipe
    a.compute_error !== nothing && return GroupedError
    (a.y in colnames(a.data.table) || a.y === nothing) ? StatPlotsRecipe : GroupedError
end

process(a::Analysis) = process(analysistype(a), a)
splitby(a::Analysis) = splitby(a.data)
orderby(a::Analysis) = orderby(a.data)


text:
module Reversible
using TropicalYao, TropicalYao.Reversible
using NiLang, NiLang.AD
using Yao
using Compose
using ..TropicalTensors
using ..TropicalTensors: _init_reg, SpinglassOptConfig
using Viznet: isconnected

export opt_config, isolve, isolve_largemem, opt_config_largemem, stack4reg

include("square.jl")
include("chimera.jl")
include("second_neighbor.jl")

function opt_config(sg::Spinglass{LT,T}) where {LT,T}
    reg = _init_reg(Tropical{T}, sg.lattice, Val(false))
    A = stack4reg(reg, cachesize_A(sg.lattice))
    B = stack4reg(reg, cachesize_B(sg.lattice))
    eng, sg, reg, A, B = isolve(T(0.0), sg, reg, A, B)
    sgg = Spinglass(sg.lattice, GVar.(sg.Js, zero(sg.Js)), GVar.(sg.hs, zero(sg.hs)))
    gres = (~isolve)(GVar(eng, T(1)), sgg, GVar(reg), GVar(A), GVar(B))
    empty!(NiLang.GLOBAL_STACK)
    return SpinglassOptConfig(sg, eng, grad.(gres[2].Js), grad.(gres[2].hs))
end

function opt_config_largemem(sg::Spinglass{LT,T}) where {LT,T}
    reg = _init_reg(Tropical{T}, sg.lattice, Val(false))
    A = stack4reg(reg, cachesize_largemem(sg.lattice))
    eng, sg, reg, A = isolve_largemem(T(0.0), sg, reg, A)
    sgg = Spinglass(sg.lattice, GVar.(sg.Js, zero(sg.Js)), GVar.(sg.hs, zero(sg.hs)))
    gres = (~isolve_largemem)(GVar(eng, T(1)), sgg, GVar(reg), GVar(A))
    empty!(NiLang.GLOBAL_STACK)
    return SpinglassOptConfig(sg, eng, grad.(gres[2].Js), grad.(gres[2].hs))
end
end


text:
using COAST
using JuMP
using Test

@testset "COAST.jl" begin
    @test loaded_COAST()==true
end

@testset "he_forward.jl" begin
    ## input
    # rdaam params
    # eqs, fanning curvilinear, ketcham et al 07
    alpha = 0.04672
    c0 = 0.39528
    c1 = 0.01073
    c2 = -65.12969
    c3 =  -7.91715
    rmr0 = 0.79
    eta_q = 0.91
    L_dist = 8.1*1e-4 # cm (!)

    # general params
    U238 = 4*1e-6
    R = 1.9872*1e-3
    L = 60*1e-4 # m
    n_T_pts = 100
    n_iter = 60.0

    # preprocess concentrations
    (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol) = conc_to_atoms_per_volume(U238,0.0,density=3.20)

    """
    check if RDAAM reproduces vanilla Durango
    from farley et al, 00
    """
    L1 = 90*1e-4
    logD0_a2 = log10(10^1.5/(L1^2)) # cm^2/s
    Ea = 32.9

    T = collect(LinRange(120.0,0.01,ceil(Int,n_T_pts)).+273.15)
    times = collect(LinRange(120.0,0.01,ceil(Int,n_T_pts)+1).*3.1558e7*1e6)
    mass_he_t1 = rdaam_forward_diffusion(alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist,0.0,0.0,0.0,R,Ea,logD0_a2,n_iter,
                  U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,L1,times...,T...)
    print("mass t1 Farley, 00 is: $mass_he_t1 \n")
    pre_he_t1 = (8*(U238_mol*exp(60*1e6*sec_in_yrs/τ38)-U238_mol)+
                 7*(U235_mol*exp(60*1e6*sec_in_yrs/τ35)-U235_mol)) # Dodson value
    print("mass t1 Farley, 00 predicted is: $pre_he_t1 \n")
    @test isapprox(mass_he_t1/pre_he_t1,1.0; rtol = 1e-2)

    """
    calculate erho
    rdaam params in diffn
    test fig 2, 2009 RDAAM
    """
    Rjoules = 8.314
    psi = 1e-13
    omega = 1e-22
    Etrap = 34*1e3 # J/mol
    E_L = 122.3*1e3 # J/mol
    L2 = 60*1e-4
    log10D0L_a2_rdaam = log10(exp(9.733)*L^2/L2^2)
    T = collect(LinRange(120.0,0.0,ceil(Int,n_T_pts)).+273.15)
    times = collect(LinRange(120.0,0.01,ceil(Int,n_T_pts)+1).*3.1558e7*1e6)
    mass_he_t2 = rdaam_forward_diffusion(alpha,c0,c1,c2,c3,0.79,eta_q,L_dist,psi,omega,Etrap,Rjoules,E_L,log10D0L_a2_rdaam,n_iter,
                   U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,L,times...,T...)
    pre_he_t2 = (8*(U238_mol*exp(43.6*1e6*sec_in_yrs/τ38)-U238_mol)+
                 7*(U235_mol*exp(43.6*1e6*sec_in_yrs/τ35)-U235_mol))# Dodson value
    print("mass of he is: $(mass_he_t2) \n")
    print("predicted mass of he is: $(pre_he_t2) \n")
    # non-redistributed value
    @test isapprox(mass_he_t2/pre_he_t2,1.0; rtol = 1e-2)

    """
    wolf,98 fig. 5 - history 3
    """
    Ea = 36.2
    logD0_a2 = log10(10^(7.7))

    T1 = collect(LinRange(60.0,60.0,ceil(Int,n_T_pts/2)).+273.15)
    T2 = collect(LinRange(15.0,15.0,floor(Int,n_T_pts/2)).+273.15)
    T = vcat(T1,T2)
    times1 = collect(LinRange(100.0,20.0,ceil(Int,n_T_pts/2)+1).*3.1558e7*1e6)
    times2 = collect(LinRange(19.9,0.01,floor(Int,n_T_pts/2)).*3.1558e7*1e6)
    times = vcat(times1,times2)
    mass_he_t3 = rdaam_forward_diffusion(alpha,c0,c1,c2,c3,0.83,eta_q,L_dist,
                                         0.0,0.0,0.0,R,Ea,logD0_a2,n_iter,
                                         U238_mol,U238_V,U235_mol,U235_V,
                                         Th232_mol,Th232_V,L,times...,T...)
    pre_he_t3 = (8*(U238_mol*exp(41.4*1e6*sec_in_yrs/τ38)-U238_mol)+ # HeFTy output
                 7*(U235_mol*exp(41.4*1e6*sec_in_yrs/τ35)-U235_mol)) # HeFTy output
    print("Wolf test measured is $(mass_he_t3)! \n")
    print("Wolf test predicted is $(pre_he_t3)! \n")
    @test isapprox(mass_he_t3/pre_he_t3,1.0; atol = 1e-2) # <1 Ma error vs
end

@testset "he_vars_constraints.jl" begin
     # basic setup
     model1 = initialize_JuMP_model("mumps")
     @variable(model1, T)
     @test num_variables(model1)==1

     # add bounds
     model2 = initialize_JuMP_model("mumps")
     @variable(model2, 273.0<=T<=400.0)
     @test num_variables(model2)==1
     @test has_lower_bound(T)
     @test has_upper_bound(T)
     print("Basic JuMP tests passed!\n")
end

# basic model with custom constraints
@testset "custom constraints" begin
    # general params
    U238 = 4*1e-6
    R = 1.9872*1e-3
    L = 60*1e-4 # m
    n_T_pts = 100
    n_iter = 60.0
     model3 = initialize_JuMP_model("mumps",print_level=1)
     time_segs = 10
     (T,set_val) = define_variables!(time_segs-1,1,model3,0.1*ones(time_segs-1))
     @test num_variables(model3)==time_segs-1
     @test register_forward_model!(time_segs,model3)==true
     register_objective_function!(time_segs-1,model3)

     ## input
     # rdaam params
     # eqs, fanning curvilinear, ketcham et al 07
     alpha = 0.04672
     c0 = 0.39528
     c1 = 0.01073
     c2 = -65.12969
     c3 =  -7.91715
     rmr0 = 0.79
     eta_q = 0.91
     L_dist = 8.1*1e-4 # cm (!)

     # general params
     U238 = 28*1e-6
     R = 1.9872*1e-3
     L = 60*1e-4 # m
     times = collect(LinRange(120.0,0.01,time_segs).*3.1558e7*1e6)

     # preprocess concentrations
     (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol) = conc_to_atoms_per_volume(U238,0.0,density=3.20)

     L1 = 90*1e-4
     logD0_a2 = log10(10^1.5/(L1^2)) # cm^2/s
     Ea = 32.9

     constr1=@NLconstraint(model3,[i = 1:1],rdaam_forward_diffusion(alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist,0.0,0.0,0.0,R,Ea,logD0_a2,n_iter,
                   U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,L,times...,T...)==3e-9)
     @NLobjective(model3,Min,mod_constraints(T...))
     optimize!(model3)
     @test termination_status(model3) == MOI.LOCALLY_SOLVED
end

@testset "one specified constraint" begin
    model4 = initialize_JuMP_model("mumps",print_level=1)
    time_segs = 10
    (T,set_dev) = define_variables!(time_segs-1,1,model4,0.1*ones(time_segs-1))
    register_objective_function!(time_segs-1,model4)
    register_forward_model!(time_segs,model4)
    ## input
    # rdaam params
    # eqs, fanning curvilinear, ketcham et al 07
    alpha = 0.04672
    c0 = 0.39528
    c1 = 0.01073
    c2 = -65.12969
    c3 =  -7.91715
    rmr0 = 0.79
    eta_q = 0.91
    L_dist = 8.1*1e-4 # cm (!)

    # general params
    U238 = 28*1e-6
    R = 1.9872*1e-3
    L = 60*1e-4 # m
    times = collect(LinRange(120.0,0.01,time_segs).*3.1558e7*1e6)

    # preprocess concentrations
    (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol) = ppm_to_atoms_per_volume(U238,0.0,density=3.20)

    L1 = 90*1e-4
    logD0_a2 = log10(10^1.5/(L1^2)) # cm^2/s
    Ea = 32.9
    density = 3.2
    He_conc = [3e-9]
    rdaam_define_constraints(model4,density,set_dev,He_conc,L,times,T,U238_mol,U235_mol,Th232_mol,U238_V,U235_V,Th232_V)
    optimize!(model4)
    @test termination_status(model4) == MOI.LOCALLY_SOLVED
end

@testset "FT_setup.jl" begin
     # Ketcham 07, 5_5M, Dpar case
    (l0m_K07_5_5Dpar,l0cm_K07_5_5Dpar)= init_track_len(calc_type="Dpar",
                                                Dpar=1.65,Ketcham07_5_5m=true)
    @test isapprox(l0m_K07_5_5Dpar,16.1,rtol=1e-3)
    @test isapprox(l0cm_K07_5_5Dpar,16.3,rtol=1e-3)

    (l0m_K07_5_5rmr0,l0cm_K07_5_5rmr0)= init_track_len(calc_type="rmr0",
                                                rmr0=1.65,Ketcham07_5_5m=true)
    @test isapprox(l0m_K07_5_5rmr0,16.14,rtol=1e-3)
    @test isapprox(l0cm_K07_5_5rmr0,16.35,rtol=1e-3)

    (l0m_K07_5_5cl_a,l0cm_K07_5_5cl_a)= init_track_len(calc_type="Cl_apfu",
                                                Cl=1.65,Ketcham07_5_5m=true)
    @test isapprox(l0m_K07_5_5cl_a,17.80,rtol=1e-3)
    @test isapprox(l0cm_K07_5_5cl_a,18.33,rtol=1e-3)

    (l0m_K07_5_5a_vol,l0cm_K07_5_5a_vol)= init_track_len(calc_type="a_len",
                                                a_vol=10.65,Ketcham07_5_5m=true)
    @test isapprox(l0m_K07_5_5a_vol,23.18,rtol=1e-3)
    @test isapprox(l0cm_K07_5_5a_vol,24.53,rtol=1e-3)

    # Donelick99 case
   (l0m_D99Dpar,l0cm_D99Dpar)= init_track_len(calc_type="Dpar",
                                               Dpar=1.65,Ketcham07_5_5m=false,
                                               Donelick99=true)
   @test isapprox(l0m_D99Dpar,16.1,rtol=1e-3)
   @test isapprox(l0cm_D99Dpar,16.3,rtol=1e-3)

   (l0m_D99rmr0,l0cm_D99rmr0)= init_track_len(calc_type="rmr0",
                                               rmr0=1.65,Donelick99=true)
   @test isapprox(l0m_D99rmr0,16.14,rtol=1e-3)
   @test isapprox(l0cm_D99rmr0,16.35,rtol=1e-3)

   (l0m_D99cl_a,l0cm_D99cl_a)= init_track_len(calc_type="Cl_apfu",
                                               Cl=1.65,Donelick99=true)
   @test isapprox(l0m_D99cl_a,17.80,rtol=1e-3)
   @test isapprox(l0cm_D99cl_a,18.33,rtol=1e-3)

   # Ketcham 07, 5_0M, Dpar case
  (l0m_K07_5_0Dpar,l0cm_K07_5_0Dpar)= init_track_len(calc_type="Dpar",
                                              Dpar=1.65,Ketcham07_5_5m=false,
                                              Ketcham07_5_0m=true)
  @test isapprox(l0m_K07_5_0Dpar,15.82,rtol=1e-3)
  @test isapprox(l0cm_K07_5_0Dpar,16.06,rtol=1e-3)

  (l0m_K07_5_0rmr0,l0cm_K07_5_0rmr0)= init_track_len(calc_type="rmr0",
                                              rmr0=1.65,Ketcham07_5_0m=true)
  @test isapprox(l0m_K07_5_0rmr0,15.94,rtol=1e-3)
  @test isapprox(l0cm_K07_5_0rmr0,16.19,rtol=1e-3)

  (l0m_K07_5_0cl_a,l0cm_K07_5_0cl_a)= init_track_len(calc_type="Cl_apfu",
                                              Cl=1.65,Ketcham07_5_0m=true)
  @test isapprox(l0m_K07_5_0cl_a,16.82,rtol=1e-3)
  @test isapprox(l0cm_K07_5_0cl_a,17.18,rtol=1e-3)

  (l0m_K07_5_0a_vol,l0cm_K07_5_0a_vol)= init_track_len(calc_type="a_len",
                                              a_vol=10.65,Ketcham07_5_0m=true)
  @test isapprox(l0m_K07_5_0a_vol,24.85,rtol=1e-3)
  @test isapprox(l0cm_K07_5_0a_vol,25.22,rtol=1e-3)

  # test custom functionality
  (l0m_cust,l0cm_cust)= init_track_len(custom=true,custom_param=1.65,
  custom_intercept=3.0,custom_slope=0.5,custom_c_intercept=2.0,
                                       custom_c_slope=1.0)
   @test isapprox(l0m_cust,(1.65*0.5+3.0),rtol=1e-3)
   @test isapprox(l0cm_cust,(1.65*1.0+2.0),rtol=1e-3)
end

# @testset "FT_forward.jl" begin
#     """
#     RDAAM type cooling
#     """
#     T = collect(LinRange(120.0,0.01,ceil(Int,n_T_pts)).+273.15)
#     times = collect(LinRange(120.0,0.01,ceil(Int,n_T_pts)+1).*3.1558e7*1e6)
#     alpha = 0.04672
#     beta = 0.0
#     c0 = 0.39528
#     c1 = 0.01073
#     c2 = -65.12969
#     c3 =  -7.91715

#     rcb2_test = FT_forward(alpha,beta,c0,c1,c2,c3,T,times)
# end


text:
# This file illustrates the simulation of Q-PSK signalling. The simulation is performed using
# waveform AWGN channel.

using Plots
using FFTW
using LinearAlgebra
using Statistics
using DigitalCommunications

# Settings
k = 2                                   # Bits per symbol
M = 2^k                                 # Constellation size
T = 1.                                  # Symbol duration
f = 1 / T                               # Symbols per second
fs = 10                                 # Sampling frequency
ts = 1 / fs                             # Sampling period
t = 0 : ts : T - ts                     # Time vector for symbol duration
l = length(t)                           # Number of samples per symbol
# g(t) = 0 ≤ t ≤ T ? sqrt(Eg / T) : 0     # Modulation pulse: rectangular pulse
# Eg = 1                                  # Modulation pulse energy: rectangular pulse
A = 1                                     # Amplitude for raised cosine pulse
g(t) = 0 ≤ t ≤ T ?  A/2 * (1 - cos(2π/T * t)) : 0       # Raised cosine modulation pulse
Eg = A^2 * 3 / 8                                        # Energy for raised cosine modulation pulse
nsymbols = Int(1e6)                     # Number of symbols

# Define basis
basis = [
    t -> g(t) / sqrt(Eg)
]
alphabet = [
    [1 + 0im],
    [0 + 1im],
    [-1 + 0im],
    [0 - 1im],
]
basisvals = map(base -> base.(t), basis)
references = map(1 : M) do mi
    Ai = alphabet[mi] * sqrt(Eg)
    sum(Ai .* basisvals)
end

# Simulate the system
ebno = collect(0 : 1 : 10)
esno = ebno .+ 10 * log10(k)
m = rand(1 : M, nsymbols)
tx = map(m) do mi
    Ai = alphabet[mi] * sqrt(Eg)
    sum(Ai .* basisvals)
end
Es = mean(txi -> sum(abs.(txi).^2) * ts, tx)
symerr = zeros(length(esno))
for i = 1 : length(symerr)
    # Generate noise
    σ = sqrt(Es * fs / 2 / (10^(esno[i] / 10)))
    n = [σ * (randn(l) + 1im * randn(l)) for i = 1 : length(tx)]
    # Corrupt transmitted signal
    rx = tx + n
    # Detect symbols
    mr = map(rx) do rxi
        argmax(map(reference -> real(rxi ⋅ reference), references))
    end
    # Calculate symbol
    ber = sum(m .!= mr) / length(m)
    symerr[i] = ber
end

# Plots
plot(ebno, symerr, marker=:circle, yscale=:log10, label="montecarlo")
plot!(ebno, berpsk.(esno, M), marker=:circle, yscale=:log10, label="theoretical")
title!("$M-PSK")
xlabel!("EbNo")
ylabel!("Pe")


text:
"""
    transform(::YeoJohnson, 𝐱)

Transform an array using Yeo-Johnson method.  The power parameter λ is derived
from maximizing a log-likelihood estimator.
"""
function transform(method::YeoJohnson, 𝐱; optim_args...)
    λ, details = lambda(method, 𝐱; optim_args...)
    @debug "YeoJohnson: estimated lambda = $λ"
    return transform(method, 𝐱, λ)
end

"""
    transform(::YeoJohnson, 𝐱, λ)

Transform an array using Yeo-Johnson method with the provided power parameter λ.
"""
function transform(::YeoJohnson, 𝐱, λ)
    𝐱′ = similar(𝐱, Float64)
    for (i, x) in enumerate(𝐱)
        if x >= 0
            𝐱′[i] = λ ≈ 0 ? log(x + 1) : ((x + 1)^λ - 1)/λ
        else
            𝐱′[i] = λ ≈ 2 ? -log(-x + 1) : -((-x + 1)^(2 - λ) - 1) / (2 - λ)
        end
    end
    return 𝐱′
end

"""
    lambda(::YeoJohnson, 𝐱; interval = (-2.0, 2.0), optim_args...)

Calculate lambda from an array using a log-likelihood estimator.
Keyword arguments:
- interval: search interval
- optim_args: keyword arguments accepted by Optim.optimize function
See also: [`log_likelihood`](@ref)
"""
function lambda(method::YeoJohnson, 𝐱; interval = (-2.0, 2.0), optim_args...)
    i1, i2 = interval
    res = optimize(λ -> -log_likelihood(method, 𝐱, λ), i1, i2; optim_args...)
    return (value=minimizer(res), details=res)
end

"""
    log_likelihood(::YeoJohnson, 𝐱, λ)

Return log-likelihood for the given array and lambda.
"""
function log_likelihood(method::YeoJohnson, 𝐱, λ)
    N = length(𝐱)
    𝐲 = transform(method, float.(𝐱), λ)
    σ² = var(𝐲, corrected = false)
    c = sum(sign.(𝐱) .* log.(abs.(𝐱) .+ 1))
    llf = -N / 2.0 * log(σ²) + (λ - 1) * c
    @debug "YeoJohnson: λ = $λ => σ²=$σ², c=$c, llf=$llf"
    return llf
end


text:
include("../imports.jl")

@testset ExtendedTestSet "Affine" begin
    P = CoordinateTransformations.Translation(10, 10)
    tfm = Project(P)

    keypoints = Keypoints([SVector(1., 1)], (50, 50))
    image = Image(rand(RGB, 50, 50))

    bounds = getbounds(keypoints)

    @testset ExtendedTestSet "Affine Keypoints" begin
        @test getprojection(tfm, getbounds(keypoints)) == P
        @test_nowarn project(P, keypoints, bounds)
        tkeypoints = apply(tfm, keypoints)
        @test itemdata(tkeypoints) == [SVector(11, 11)]
        testprojective(tfm)
    end

    @testset ExtendedTestSet "Affine Image" begin
        @test_nowarn project(P, image, bounds)
        timage = apply(tfm, image)
        @test itemdata(timage)[11, 11] ≈ itemdata(image)[1, 1]
    end

    tfmcropped = tfm |> CenterCrop((50, 50))

    @testset ExtendedTestSet "Cropped, Image" begin
        tfmcropped = tfm |> CenterCrop((50, 50))
        @test_nowarn apply(tfmcropped, keypoints)
        tkeypoints = apply(tfmcropped, keypoints)
        @test getbounds(tkeypoints).rs[1] |> first == 11
        @test itemdata(tkeypoints)[1] == SVector(11, 11)
    end

    @testset ExtendedTestSet "Cropped with PinOrigin, Image" begin
        tfmcropped = tfm |> PinOrigin() |> CenterCrop((50, 50))
        @test_nowarn apply(tfmcropped, keypoints)
        tkeypoints = apply(tfmcropped, keypoints)
        @test getbounds(tkeypoints).rs[1] |> first == 1
    end

    @testset ExtendedTestSet "Composed" begin
        tfms = Project(Translation(10, 10)) |> Project(Translation(-10, -10))
        randstate = getrandstate(tfms)
        @test getprojection(tfms, getbounds(image); randstate) == Translation(0, 0)
        @test itemdata(apply(tfms, image)) == itemdata(image)
        @test itemdata(apply(tfms, keypoints)) == itemdata(keypoints)
    end

    @testset ExtendedTestSet "ScaleKeepAspect" begin
        image = Image(zeros(RGB, 100, 50))
        keypoints = Keypoints([SVector(0., 0)], (50, 50))

        imbounds = getbounds(image)
        kpbounds = getbounds(keypoints)
        tfm1 = ScaleKeepAspect((50, 50))
        tfm2 = ScaleKeepAspect((25, 25))
        testprojective(tfm1)
        testprojective(tfm2)
    end

    @testset ExtendedTestSet "Scale" begin
        keypoints = Keypoints([SVector(0., 0)], (50, 50))
        image = Image(rand(RGB, 50, 50))

        @testset ExtendedTestSet "ScaleFixed" begin
            tfm = ScaleFixed((25, 25))
            @test_nowarn apply(tfm, image)
            @test_nowarn apply(tfm, keypoints)
            timage = apply(tfm, image)
            tkeypoints = apply(tfm, keypoints)
            @test length.(getbounds(timage).rs) == (25, 25)
            @test getbounds(timage) == getbounds(tkeypoints)
            testprojective(tfm)

        end
        @testset ExtendedTestSet "ScaleRatio" begin
            tfm = ScaleRatio((1/2, 1/2))
            @test_nowarn apply(tfm, image)
            @test_nowarn apply(tfm, keypoints)
            timage = apply(tfm, image)
            tkeypoints = apply(tfm, keypoints)
            @test getbounds(timage).rs == (0:25, 0:25)
            @test getbounds(timage) == getbounds(tkeypoints)
            testprojective(tfm)
        end

        @testset ExtendedTestSet "ScaleKeepAspect" begin
            tfm = ScaleKeepAspect((32, 32))
            img = rand(RGB{N0f8}, 64, 96)
            @test apply(tfm, Image(img)) |> itemdata |> size == (32, 48)

        end
    end

    @testset ExtendedTestSet "Rotate" begin
        tfm = Rotate(10)
        image = Image(rand(RGB, 50, 50))
        @test_nowarn apply(tfm, image)
        P = DataAugmentation.getprojection(tfm, getbounds(image))
        @test P isa AffineMap
        @test P.linear.mat[1] isa Float32
    end

    @testset ExtendedTestSet "Rotate" begin
        tfm = Zoom((0.1, 2.))
        image = Image(rand(RGB, 50, 50))
        @test_nowarn apply(tfm, image)
    end

    @testset ExtendedTestSet "Reflect" begin
        tfm = Reflect(10)
        testprojective(tfm)
    end

    @testset ExtendedTestSet "CenterResizeCrop" begin
        tfm = CenterResizeCrop((25, 40))
        keypoints = Keypoints([SVector(0., 0)], (50, 50))
        image = Image(rand(RGB, 50, 50))
        timage, tkeypoints = apply(tfm, (image, keypoints))
        @test getbounds(timage).rs == (1:25, 1:40)
        testprojective(Project(IdentityTransformation()) |> CenterCrop((12, 12)))
    end

    @testset ExtendedTestSet "CroppedAffine inplace" begin
        image1 = Image(rand(RGB, 100, 100))
        image2 = Image(rand(RGB, 100, 100))
        tfm = ScaleFixed((50, 50))
        buffer = makebuffer(tfm, image1)

        @test_nowarn apply!(buffer, tfm, image2)
    end

    @testset ExtendedTestSet "`RandomCrop` correct indices" begin
        # Flipping and cropping should be the same as reverse-indexing
        # the flipped dimension
        tfm = FlipX() |> RandomCrop((64, 64)) |> PinOrigin()
        img = rand(RGB, 64, 64)
        item = Image(img)
        titem = apply(tfm, item)
        timg = itemdata(titem)
        rimg = img[:, end:-1:1]
        @test titem.data == rimg
    end
end

@testset ExtendedTestSet "Big pipeline" begin
    @testset ExtendedTestSet "2D" begin
        tfms = compose(
            Rotate(10),
            FlipX(), FlipY(),
            ScaleRatio((.8, .8)),
            RandomCrop((10, 10)),
            WarpAffine(0.1),
            Zoom((1., 1.2))
        )
        testprojective(tfms)
    end

    @testset ExtendedTestSet "3D" begin
        sz = (50, 50, 50)
        items = (
            Image(rand(RGB, sz)),
            Keypoints(rand(SVector{3, Float32}, 50), sz),
            MaskBinary(rand(Bool, sz)),
            MaskMulti(rand(1:8, sz)),
        )

        tfms = compose(
            ScaleRatio((.8, .8, .8)),
            ScaleKeepAspect((12, 10, 10)),
            RandomCrop((10, 10, 10))
        )
        testprojective(tfms, items)

        buf = makebuffer(tfms, items)
        @test_nowarn apply!(buf, tfms, items)
    end
end


text:
#Example lowpass/bandpass real/quadrature modulator design
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import Printf: @sprintf
j=im

"""`result = dsexample1(dsm; ampdB::Float64=-3.0, ftest=nothing, LiveDemo=false)`

Example lowpass/bandpass real/quadrature modulator design.

# Output `result`: `Dict` of:
 - :nlev, :NTF, :ABCD, :umax, :amp, :snr, :peak_snr
 - :coeff (a NamedTuple of: a, g, b c)
"""
function dsexample1(dsm::AbstractDSM; ampdB::Float64=-3.0, ftest=nothing, LiveDemo::Bool=false)
	#Computed defaults
	if isnothing(ftest)
		ftest = default_ftest(dsm)
	end

	#Derived parameters
	nlev = dsm.M + 1
	dsmtypestr = str_modulatortype(dsm)
	println("*** ", dsmtypestr, " example... ")
	#TODO: support plot attributes?:
	sizes = Dict{Symbol, Any}(
		:lw => 1, #LineWidth
		:ms => 5, #MarkerSize
		:fs => 12, #FontSize
		:fw => "normal", #FontWeight
	)
	if LiveDemo #Overwrite settings:
		push!(sizes,
			:lw => 2,
			:ms => 6,
			:fw => "bold",
		)
	end

	LiveDemo && (@info("Realizing NTF..."); flush(stdout); flush(stderr))
#-------------------------------------------------------------------------------
	local NTF, ABCD
	if !isquadrature(dsm)
		NTF = synthesizeNTF(dsm)
		a,g,b,c = realizeNTF(NTF, dsm.form)
		z0 = exp(2π*j*dsm.f0)
		b = [abs(b[1]+b[2]*(1-z0)) zeros(1,length(b)-1)] #Use a single feed-in for the input
		ABCD = stuffABCD(a,g,b,c,dsm.form)
	else
		NTF = synthesizeQNTF(dsm)
		ABCD = realizeQNTF(NTF, dsm.form, 1)
	end
	LiveDemo && println("\tdone.")
	plot = documentNTF(dsm, ABCD, sizes=sizes, frespOnly=false)
	saveimage(:png, "dsexample1_NTF.png", plot, AR=2/1, width=900)
	displaygui(plot)

	LiveDemo && (@info("Performing time-domain simulations..."); flush(stdout); flush(stderr))
#-------------------------------------------------------------------------------
	N = 100
	t = 0:N-1
	#Do not use genTestTone() (N too small; ftest gets rounded to 0Hz).
	#(u, iftest) = genTestTone(dsm, ampdB, ftest, N=100)
	local simresult
	if !isquadrature(dsm)
		u = undbv(ampdB)*dsm.M*sin.( 2π*ftest*t )
		simresult = simulateDSM(u, ABCD, nlev=nlev)
	else
		u = undbv(ampdB)*dsm.M*exp.( 2π*j*ftest*t );
		simresult = simulateQDSM(u, ABCD, nlev=nlev)
	end
	LiveDemo && println("\tdone.")
	plot = plotModTransient(u, simresult.v, legend=false)
		ylim = (dsm.M+0.25)
		set(plot, xyaxes=set(ymin=-ylim, ymax=ylim))
	displaygui(plot)

	#Example spectrum
#-------------------------------------------------------------------------------
	plot = plotExampleSpectrum(dsm, NTF, ampdB=ampdB, ftest=ftest, sizes=sizes)
	displaygui(plot)

	#SQNR plot
#-------------------------------------------------------------------------------
	snrinfo = calcSNRInfo(dsm, NTF=NTF)
	plot = plotSNR(snrinfo, dsm)
		#TODO: Apply sizes
	displaygui(plot)

	if isquadrature(dsm) #Example I/Q mismatch
		error("Quadrature section not ported over.")
	end

	umax = nothing
	coeffs = nothing
	if !isquadrature(dsm) #Dynamic range scaling
		LiveDemo && (@info("Performing dynamic range scaling..."); flush(stdout); flush(stderr))
#-------------------------------------------------------------------------------
		ABCD0 = ABCD
		ABCD, umax = scaleABCD(ABCD0, nlev=nlev, f=dsm.f0)
		a,g,b,c = mapABCD(ABCD, dsm.form)
		coeffs = (a=a, g=g, b=b, c=c)
		LiveDemo && println("\tdone.")

		LiveDemo && (@info("Verifying dynamic range scaling..."); flush(stdout); flush(stderr))
#-------------------------------------------------------------------------------
		N = 10^4; N0 = 50
		u = range(0, stop=0.95*umax, length=30)
		test_tone = cos.(2π*dsm.f0 * (0:N-1))
		test_tone[1:N0] = test_tone[1:N0] .* (0.5 .- 0.5*cos.(2π/N0 * (0:N0-1)))

		plot = plotStateMaxima(u, ABCD, test_tone, nlev=nlev)
			set(plot, RSDeltaSigmaPort.linlin, xyaxes=set(ymin=0, ymax=1))
		LiveDemo && println("\tdone.")
		displaygui(plot)
	end

	results = Dict{Symbol, Any}(
		:nlev => nlev,
		:NTF => NTF,
		:ABCD => ABCD,
		:SNR_vs_amp_sim => snrinfo.vs_amp_sim,
		:peak_snr => snrinfo.peak[2],
		:umax => umax,
		:coeffs => coeffs
	)
	return results
end

#Last line


text:
#Methods operating on different special matrix types

#Interconversion between special matrix types
convert{T}(::Type{Bidiagonal}, A::Diagonal{T})=Bidiagonal(A.diag, zeros(T, size(A.diag,1)-1), true)
convert{T}(::Type{SymTridiagonal}, A::Diagonal{T})=SymTridiagonal(A.diag, zeros(T, size(A.diag,1)-1))
convert{T}(::Type{Tridiagonal}, A::Diagonal{T})=Tridiagonal(zeros(T, size(A.diag,1)-1), A.diag, zeros(T, size(A.diag,1)-1))
convert(::Type{Triangular}, A::Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal))=Triangular(full(A))
convert(::Type{Matrix}, D::Diagonal) = diagm(D.diag)

function convert(::Type{Diagonal}, A::Union(Bidiagonal, SymTridiagonal))
    all(A.ev .== 0) || throw(ArgumentError("Matrix cannot be represented as Diagonal"))
    Diagonal(A.dv)
end

function convert(::Type{SymTridiagonal}, A::Bidiagonal)
    all(A.ev .== 0) || throw(ArgumentError("Matrix cannot be represented as SymTridiagonal"))
    SymTridiagonal(A.dv, A.ev)
end

convert{T}(::Type{Tridiagonal}, A::Bidiagonal{T})=Tridiagonal(A.isupper?zeros(T, size(A.dv,1)-1):A.ev, A.dv, A.isupper?A.ev:zeros(T, size(A.dv,1)-1))

function convert(::Type{Bidiagonal}, A::SymTridiagonal)
    all(A.ev .== 0) || throw(ArgumentError("Matrix cannot be represented as Bidiagonal"))
    Bidiagonal(A.dv, A.ev, true)
end

function convert(::Type{Diagonal}, A::Tridiagonal)
    all(A.dl .== 0) && all(A.du .== 0) || throw(ArgumentError("Matrix cannot be represented as Diagonal"))
    Diagonal(A.d)
end

function convert(::Type{Bidiagonal}, A::Tridiagonal)
    if all(A.dl .== 0) return Bidiagonal(A.d, A.du, true)
    elseif all(A.du .== 0) return Bidiagonal(A.d, A.dl, true)
    else throw(ArgumentError("Matrix cannot be represented as Bidiagonal"))
    end
end

function convert(::Type{SymTridiagonal}, A::Tridiagonal)
    all(A.dl .== A.du) || throw(ArgumentError("Matrix cannot be represented as SymTridiagonal"))
    SymTridiagonal(A.d, A.dl)
end

function convert(::Type{Diagonal}, A::Triangular)
    full(A) == diagm(diag(A)) || throw(ArgumentError("Matrix cannot be represented as Diagonal"))
    Diagonal(diag(A))
end

function convert(::Type{Bidiagonal}, A::Triangular)
    fA = full(A)
    if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1)
        return Bidiagonal(diag(A), diag(fA,1), true)
    elseif fA == diagm(diag(A)) + diagm(diag(fA, -1), -1)
        return Bidiagonal(diag(A), diag(fA,-1), true)
    else
        throw(ArgumentError("Matrix cannot be represented as Bidiagonal"))
    end
end

convert(::Type{SymTridiagonal}, A::Triangular) = convert(SymTridiagonal, convert(Tridiagonal, A))

function convert(::Type{Tridiagonal}, A::Triangular)
    fA = full(A)
    if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1) + diagm(diag(fA, -1), -1)
        return Tridiagonal(diag(fA, -1), diag(A), diag(fA,1))
    else
        throw(ArgumentError("Matrix cannot be represented as Tridiagonal"))
    end
end

#Constructs two method definitions taking into account (assumed) commutativity
# e.g. @commutative f{S,T}(x::S, y::T) = x+y is the same is defining
#     f{S,T}(x::S, y::T) = x+y
#     f{S,T}(y::T, x::S) = f(x, y)
macro commutative(myexpr)
    @assert myexpr.head===:(=) || myexpr.head===:function #Make sure it is a function definition
    y = copy(myexpr.args[1].args[2:end])
    reverse!(y)
    reversed_call = Expr(:(=), Expr(:call,myexpr.args[1].args[1],y...), myexpr.args[1])
    esc(Expr(:block, myexpr, reversed_call))
end

for op in (:+, :-)
    SpecialMatrices = [:Diagonal, :Bidiagonal, :Tridiagonal, :Triangular, :Matrix]
    for (idx, matrixtype1) in enumerate(SpecialMatrices) #matrixtype1 is the sparser matrix type
        for matrixtype2 in SpecialMatrices[idx+1:end] #matrixtype2 is the denser matrix type
            @eval begin #TODO quite a few of these conversions are NOT defined...
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for  matrixtype1 in (:SymTridiagonal,)                      #matrixtype1 is the sparser matrix type
        for matrixtype2 in (:Tridiagonal, :Triangular, :Matrix) #matrixtype2 is the denser matrix type
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for matrixtype1 in (:Diagonal, :Bidiagonal) #matrixtype1 is the sparser matrix type
        for matrixtype2 in (:SymTridiagonal,)   #matrixtype2 is the denser matrix type
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
export VortexConfig, state_to_lagrange, lagrange_to_state

"""
    VortexConfig

An immutable structure representing a corner (`N=1`)/wedge (`N=2`)

## Fields
- `Nv::Int64`: number of vortices
- `U::ComplexF64`: freestream velocity
- `ss::Array{ComplexF64,1}`: sensor locations
- `Δt::Float64`: time step
- `δ::Float64` : blob radius

- `ϵX::Float64`: position noise std
- `ϵΓ::Float64`: circulation noise std

- `β::Float64`: multiplicative inflation

- `ϵY::Float64`: observation noise std
"""
struct VortexConfig

    "Number of vortices"
    Nv::Int64

    "Freestream velocity"
    U::ComplexF64

    "Sensor locations"
    ss::Array{ComplexF64,1}

    "Time step"
    Δt::Float64

    "Blob radius"
    δ::Float64

    "Position noise std"
    ϵX::Float64

    "Circulation noise std"
    ϵΓ::Float64

    "Multiplicative inflation"
    β::Float64

    "Observation noise std"
    ϵY::Float64
end

function state_to_lagrange(state::AbstractVector{Float64}, config::VortexConfig; isblob::Bool=true)
    Nv = config.Nv

    zv = (state[3i-2] + im*state[3i-1] for i in 1:Nv)
    Γv = (state[3i] for i in 1:Nv)

    if isblob == true #return collection of regularized point vortices
        δ = config.δ
        δv = (δ for i in 1:Nv)
        blobs₊ = Nv > 0 ? map(Vortex.Blob, zv, Γv, δv) : Vortex.Blob{Float64, Float64}[]
        blobs₋ = Nv > 0 ? map((z, Γ, δ) -> Vortex.Blob(conj(z), -Γ, δ), zv, Γv, δv) : Vortex.Blob{Float64, Float64}[]

        return (blobs₊, blobs₋)
    else #return collection of point vortices
        points₊ = Nv > 0 ? map(Vortex.Point, zv, Γv) : Vortex.Point{Float64, Float64}[]
        points₋ = Nv > 0 ? map((z, Γ) -> Vortex.Point(conj(z), -Γ), zv, Γv) : Vortex.Point{Float64, Float64}[]

        return (points₊, points₋)
    end
end

function lagrange_to_state(source, config::VortexConfig)
    Nv = length(source[1])

    @assert Nv == config.Nv

    states = Vector{Float64}(undef, 3*Nv)

    for i=1:Nv
        bi = source[1][i]
        states[3i-2] = real(bi.z)
        states[3i-1] = imag(bi.z)
        states[3i] =  circulation(bi)
    end

    return states
end


text:
using Documenter, JuliaCipher

#makedocs()

makedocs(
    ...,
    format = :html,
    sitename = "JuliaCipher",
)

deploydocs(
    repo   = "github.com/FLCN/JuliaCipher.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
@with_kw mutable struct TrainingParams
    loss
    optimizer = ADAM(3f-4)
    regularizer = (π) -> 0
    batch_size = 128
    epochs = 80
    update_every = 1
    early_stopping = (info) -> false
    name = ""
    max_batches=Inf
end

Flux.Optimise.train!(π::N, loss::Function, p::TrainingParams; info = Dict()) where {N <: Policy} = train!(Flux.params(π), (;info) -> loss(info = info) + p.regularizer(π), p, info=info)

function Flux.Optimise.train!(θ, loss::Function, p::TrainingParams; info = Dict())
    l, back = Flux.pullback(() -> loss(info = info), θ)
    typeof(l) == Float64 && @error "Float64 loss found: computation in double precision may be slow"
    grad = back(1f0)
    gnorm = norm(grad, p=2)
    isnan(gnorm) && error("NaN detected! Loss: $l")
    Flux.update!(p.optimizer, θ, grad)
    info[string(p.name, "loss")] = l
    info[string(p.name, "grad_norm")] = gnorm
    info
end

# Train with minibatches and epochs
function batch_train!(π, p::TrainingParams, 𝒫, 𝒟::ExperienceBuffer...; info=Dict(), π_loss=π)
    infos = [] # stores the aggregated info for each epoch
    total_batches = 0
    for epoch in 1:p.epochs
        minibatch_infos = [] # stores the info from each minibatch

        # Shuffle the experience buffers
        for D in 𝒟
            shuffle!(D)
        end

        # Call train for each minibatch
        partitions = [partition(1:length(D), p.batch_size) for D in 𝒟]
        for indices in zip(partitions...)
            mbs = [minibatch(D, i) for (D, i) in zip(𝒟, indices)]
            push!(minibatch_infos, train!(π, (;kwargs...)->p.loss(π_loss, 𝒫, mbs...; kwargs...), p, info=info))
            total_batches += 1
            if total_batches >= p.max_batches
                break
            end

            if p.early_stopping([infos...,  aggregate_info(minibatch_infos)])
                info[:early_stop_batch] = total_batches
                break
            end
        end
        push!(infos, aggregate_info(minibatch_infos))
        if p.early_stopping(infos)
            break
        end

    end
    merge!(info, aggregate_info(infos))
end


text:
# Food insecurity data from Census Household Pulse Survey =================================================
# NOTE: Percentage calculated from Table 2.b of each period (Week 1 pre-COVID, Week 12, and Week 20)
#       https://www.census.gov/programs-surveys/household-pulse-survey/data.html
#   Percentage calculation of total household with children
df_hps = DataFrame();
df_hps[:Time]   = ["Pre-pandemic", "Pre-pandemic", "Jul. 16 - Jul. 21", "Jul. 16 - Jul. 21", "Nov. 25 - Dec. 7", "Nov. 25 - Dec. 7"];
df_hps[:Income] = ["Less than \$25,000", "\$25,000 - \$34,999", "Less than \$25,000", "\$25,000 - \$34,999", "Less than \$25,000", "\$25,000 - \$34,999"];
df_hps[:value]  = [27.141, 19.324, 30.738, 20.408, 31.511, 22.957];

plt_hps = plot(df_hps, x="Time", y="value", color="Income", Geom.line,
    Guide.colorkey(title="Annual income"
    ),
    Guide.xlabel("Survey period (2020 Calendar year)"), Guide.ylabel("Household with food shortage, percent"),
)
pngout = SVG("plt_hps.svg", 5inch, 4inch)
draw(pngout, plt_hps)


text:
fails = 0
@context "foo is foo"
f = jtest( 1==1,
           2==2)
fails += f

@context "foo is not bar"
f = jtest( 3==3,
           4==4)
fails += f


text:
using Submodular
using Convex
using FactCheck

Tol = 1e-3

facts("Operators") do

  context("affine projection") do
    S = [3 2; 1 2]
    w1 = [2, 2]
    w2 = [3, 2]
    w3 = [5, 0]
    @fact aff_proj(w1, S) --> roughly([2, 2], Tol)
    @fact aff_proj(w2, S) --> roughly([2.5, 1.5], Tol)
    @fact aff_proj(w3, S) --> roughly([4.5, -0.5], Tol)
  end

  context("extreme point") do
    S = SetVariable(4)
    F = card([40, 70, 90, 100], S)
    @fact extreme_point(F, [3, 2, 4, 1]) --> roughly([10, 30, 40, 20], Tol)
  end

  context("gradient atom") do
    x = Variable(2)
    f = norm(x)
    g = grad(f)
    a = rand(2)
    @fact evaluate(g, a) * norm(a) --> roughly(a, Tol)
  end

  context("greedy algorithm") do
    S = SetVariable(4)
    F = card([40, 70, 90, 100], S)
    @fact greedy(F, [35, 25, 45, 15]) --> roughly([30, 20, 40, 10], Tol)
  end

  context("Lovasz extension") do
    # the Lovasz extension of a submodular function defined on a set of two elements
    S = SetVariable(4)
    F = modular([1, 2, 3, 4], S)
    x = Variable(4)
    f = lovasz(F, x)
    @fact vexity(f) --> AffineVexity()
    @fact evaluate(f, [1.5, 3.5, 2.5, 4.5]) --> roughly(34, Tol)
  end

  context("Lovasz extension of absolute value") do
    # the Lovasz extension of a submodular function defined on a set of two elements
    S = SetVariable(4)
    F = modular([1, 2, 3, 4], S)
    x = Variable(4)
    f = lovasz(F, abs(x))
    @fact vexity(f) --> ConvexVexity()
    @fact evaluate(f, [-1.5, 3.5, -2.5, 4.5]) --> roughly(34, Tol)
  end

end