text:
################################################################################
### FUNCTIONS FOR PREPACKAGED PLOTS
"""
    SSNetworkPlot(A,N,P,D,S,L,f; CaseDataDirectory=pwd(), PrintSpacer="*"^50, TestMode=false)

# Arguments

A - arc structure from ImportFunctions
N - node structure from ImportFunctions
P - product structure from ImportFunctions
D - demand structure from ImportFunctions
S - supply structure from ImportFunctions
L - technology site structure from ImportFunctions
f - transportation flows from JuMP model output structure; Output.f
CaseDataDirectory=pwd(): (optional keyword arg) file directory for supply chain case study data; defaults to current Julia directory if not specified
PrintSpacer="*"^50: (optional keyword arg): Change the characters used to space out the report
```
A simple network plot showing nodes, arcs, stakeholders, and product flows
"""
function SSNetworkPlot(A,N,P,D,S,L,f; CaseDataDirectory=pwd(), PrintSpacer="*"^50, TestMode=false)
    ################################################################################
    ### CHECK FOR OUTPUT DIRECTORY
    PlotFolder = "/plots"

    # Add records folder if not present
    if !isdir(CaseDataDirectory*PlotFolder)
        mkdir(CaseDataDirectory*PlotFolder)
    end

    ################################################################################
    ### SPECS
    fig_w = 600
    fig_h = 600
    dlb = false # default legend boolean
    dfs = 12 # default font size
    dlw = 5 # default line weight
    dlc = :black # default line color
    dms = 5 # default marker size
    dmc = :white # default marker color
    dmsw = 5 # default marker stroke weight
    dmsc = :black # default marker stroke color

    # Node custom settings
    nms = 13 # node marker size

    # Transport custom settings
    tlw = 2 # tranposrt line weight
    tlc = :white # transport line color

    # Supplier custom settings
    sms = 11 # supplier marker size
    smc = :red # supplier marker color
    smsw = 1 # supplier marker stroke weight
    smsc = :white # supplier marker stroke color

    # Technology custom settings
    lms = 8 # technology marker size
    lmc = :gold # technology marker color
    lmsw = 1 # technology marker stroke weight
    lmsc = :white # technology marker stroke color

    # Consumer custom settings
    cms = 5 # consumer marker size
    cmc = :blue # consumer marker color
    cmsw = 1 # consumer marker stroke weight
    cmsc = :white # consumer marker stroke color

    ################################################################################
    ### SETUP NODES AND ARCS ON BASE PLOT
    # Arrays for arc longitudes and lattitudes; these will be fed into plot()
    Alons = []
    Alats = []
    for a in A.ID
        push!(Alons,[N.lon[A.node_s[a]],N.lon[A.node_r[a]]])
        push!(Alats,[N.lat[A.node_s[a]],N.lat[A.node_r[a]]])
    end
    if !TestMode # skip plots during testing
        plt = plot(Alons,Alats,
            color = dlc,
            lw = dlw,
            legend = dlb)
    end

    # Arrays for node longitudes and lattitudes; these will be fed into plot()
    Nlons = []
    Nlats = []
    for n in N.ID
        push!(Nlons, N.lon[n])
        push!(Nlats, N.lat[n])
    end
    if !TestMode # skip plots during testing
        scatter!(Nlons,Nlats,
            ms = nms,
            mc = dmc,
            msw = dmsw,
            msc = dmsc,
            legend = dlb)
    end

    ################################################################################
    ### PLOT FLOWS (SO THEY APPEAR UNDER STAKEHOLDER NODES)
    flons = []
    flats = []
    for a in A.ID
        for p in P.ID
            if f[a,p] > 0
                push!(flons,[N.lon[A.node_s[a]],N.lon[A.node_r[a]]])
                push!(flats,[N.lat[A.node_s[a]],N.lat[A.node_r[a]]])
            end
        end
    end
    if !TestMode # skip plots during testing
        plot!(flons,flats,
            arrow=(:closed, 2.0),
            color = tlc,
            lw = tlw,
            legend = dlb)
    end

    ################################################################################
    ### PLOT SUPPLIER, TECHNOLOGY, CONSUMER LOCATIONS & FLOWS (IN THAT ORDER!)
    # Suppliers
    Slons = []
    Slats = []
    for i in S.ID
        push!(Slons, N.lon[S.node[i]])
        push!(Slats, N.lat[S.node[i]])
    end
    if !TestMode # skip plots during testing
        scatter!(Slons,Slats,
            ms = sms,
            mc = smc,
            msw = smsw,
            msc = smsc,
            legend = dlb)
    end

    # Technologies
    Llons = []
    Llats = []
    for l in L.ID
        push!(Llons, N.lon[L.node[l]])
        push!(Llats, N.lat[L.node[l]])
    end
    if !TestMode # skip plots during testing
        scatter!(Llons,Llats,
            ms = lms,
            mc = lmc,
            msw = lmsw,
            msc = lmsc,
            legend = dlb)
    end

    # Consumers
    Dlons = []
    Dlats = []
    for j in D.ID
        push!(Dlons, N.lon[D.node[j]])
        push!(Dlats, N.lat[D.node[j]])
    end
    if !TestMode # skip plots during testing
        scatter!(Dlons,Dlats,
            ms = cms,
            mc = cmc,
            msw = cmsw,
            msc = cmsc,
            legend = dlb)
    end

    # Update settings
    if !TestMode # skip plots during testing
        plot!(size=(fig_w,fig_h),lims=:round)
    end

    ################################################################################
    ### SAVE PLOT
    if !TestMode # skip plots during testing
        savefig(plt,CaseDataDirectory*PlotFolder*"/NetworkPlot.png")
    end

    ################################################################################
    ### END OF CODE ###
    println(PrintSpacer,"\nPlot Saved to "*CaseDataDirectory*PlotFolder*"!\n",PrintSpacer)
end


text:
export to_function
using Core: GlobalRef

const RESERVED = [:gate, :ctrl, :measure, :register]

function to_function(m::Module, ex)
    # parse macros first
    ex = eval_stmts(m, ex)
    ex = to_control(ex)
    ex = to_measure(ex)
    ex = to_gate_location(ex)

    quote
        $(Expr(:meta, :register, :new, gensym(:register)))
        $ex
        # force return nothing if no return declared
        return
    end
end

function eval_stmts(m, ex)
    ex isa Expr || return ex
    if ex.head === :$
        return Base.eval(m, ex.args[1])
    else
        return Expr(ex.head, map(x -> eval_stmts(m, x), ex.args)...)
    end
end

function is_gate_location(ex)
    ex isa Expr || return false
    return (ex.head === :call) && (ex.args[1] == :(=>))
end

function to_gate_location(ex)
    ex isa Expr || return ex
    if is_gate_location(ex)
        location = ex.args[2]
        gate = ex.args[3]
        return Expr(:call, GlobalRef(Compiler, :gate), gate, location)
    end

    if ex.head in [:block, :if, :for, :macrocall] #= make @inbounds etc. work =#
        return Expr(ex.head, map(to_gate_location, ex.args)...)
    end
    return ex
end

function to_control(ex)
    ex isa Expr || return ex
    if (ex.head === :macrocall) && (ex.args[1] == Symbol("@ctrl"))
        length(ex.args) == 4 || throw(ParseError("@ctrl expect 2 argument, got $(length(ex.args)-2)"))

        ctrl_location = ex.args[3]
        gate_ex = ex.args[4]

        if is_gate_location(gate_ex)
            location = gate_ex.args[2]
            gate = gate_ex.args[3]
            return Expr(:call, GlobalRef(Compiler, :ctrl), gate, location, ctrl_location)
        else
            throw(ParseError("@ctrl expect location=>gate at 2nd argument, got $gate_ex"))
        end
    end

    return Expr(ex.head, map(to_control, ex.args)...)
end

function to_measure(ex)
    ex isa Expr || return ex

    if (ex.head === :macrocall) && (ex.args[1] == Symbol("@measure"))
        length(ex.args) > 2 || throw(ParseError("@measure expect at least 1 argument, got $ex"))
        args = []
        parameters = []
        for each in ex.args[3:end]
            if is_measure_kwarg(each)
                push!(parameters, Expr(:kw, each...))
            else
                push!(args, each)
            end
        end

        length(parameters) <= 1 ||
            throw(ParseError("@measure takes only 1 keyword argument, got $(length(parameters))"))
        return Expr(:call, GlobalRef(Compiler, :measure), Expr(:parameters, parameters...), args...)
    end

    return Expr(ex.head, map(to_measure, ex.args)...)
end

function is_measure_kwarg(ex)
    ex isa Expr || return false
    ex.head == :(=) && ex.args[1] in [:reset_to, :remove]
end

ext: jl


text:
abstract type NotEqualConstraint <: Constraint end

"""
    NotEqualConstant(x::SeaPearl.IntVar, v::Int)

Inequality constraint, `x != v`
"""
struct NotEqualConstant <: NotEqualConstraint
    x       ::AbstractIntVar
    v       ::Int
    active  ::StateObject{Bool}
    function NotEqualConstant(x, v, trailer)
        constraint = new(x, v, StateObject(true, trailer))
        addOnDomainChange!(x, constraint)
        return constraint
    end
end

"""
    propagate!(constraint::NotEqualConstant, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`NotEqualConstant` propagation function. Basically remove the constant value from the domain of the variable.
"""
function propagate!(constraint::NotEqualConstant, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    setValue!(constraint.active, false)

    if constraint.v in constraint.x.domain
        addToPrunedDomains!(prunedDomains, constraint.x, remove!(constraint.x.domain, constraint.v))
        union!(toPropagate, constraint.x.onDomainChange)
        pop!(toPropagate, constraint)
        return !isempty(constraint.x.domain)
    end
    return true
end

function Base.show(io::IO, ::MIME"text/plain", con::NotEqualConstant)
    println(io, typeof(con), ": ", con.x.id, " != ", con.v, ", active = ", con.active)
    println(io, "   ", con.x)
end

function Base.show(io::IO, con::NotEqualConstant)
    print(io, typeof(con), ": ", con.x.id, " != ", con.v)
end

"""
    NotEqual(x::SeaPearl.IntVar, y::SeaPearl.IntVar)

Inequality constraint between two variables, stating that `x != y`.
"""
struct NotEqual <: NotEqualConstraint
    x       ::SeaPearl.AbstractIntVar
    y       ::SeaPearl.AbstractIntVar
    active  ::StateObject{Bool}

    function NotEqual(x::SeaPearl.AbstractIntVar, y::SeaPearl.AbstractIntVar, trailer::Trailer)
        constraint = new(x, y, StateObject(true, trailer))
        addOnDomainChange!(x, constraint)
        addOnDomainChange!(y, constraint)
        return constraint
    end
end

"""
    propagate!(constraint::NotEqual, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`NotEqual` propagation function.
"""
function propagate!(constraint::NotEqual, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    if !constraint.active.value
        return false
    end

    if isempty(constraint.x.domain) || isempty(constraint.y.domain)
        return false
    end

    if isbound(constraint.x)
        setValue!(constraint.active, false)
        pruned = remove!(constraint.y.domain, maximum(constraint.x.domain))
        if isempty(constraint.y.domain)
            return false
        end
        if !isempty(pruned)
            triggerDomainChange!(toPropagate, constraint.y)
            addToPrunedDomains!(prunedDomains, constraint.y, pruned)
        end
        return true
    end
    if isbound(constraint.y)
        setValue!(constraint.active, false)
        pruned = remove!(constraint.x.domain, maximum(constraint.y.domain))
        if isempty(constraint.x.domain)
            return false
        end
        if !isempty(pruned)
            triggerDomainChange!(toPropagate, constraint.x)
            addToPrunedDomains!(prunedDomains, constraint.x, pruned)
        end
        return true
    end
    return true
end

variablesArray(constraint::NotEqual) = [constraint.x, constraint.y]

function Base.show(io::IO, ::MIME"text/plain", con::NotEqual)
    println(io, typeof(con), ": ", con.x.id, " ≠ ", con.y.id, ", active = ", con.active)
    println(io, "   ", con.x)
    print(io, "   ", con.y)
end

function Base.show(io::IO, con::NotEqual)
    print(io, typeof(con), ": ", con.x.id, " ≠ ", con.y.id)
end


text:
module SimpleValueGraphs

using LightGraphs
using LightGraphs.SimpleGraphs: AbstractSimpleGraph, AbstractSimpleEdge, SimpleEdge

using Base: OneTo

import Base: eltype, show, reverse, iterate, length, replace_in_print_matrix, getindex, size, zero, Tuple
import SparseArrays: blockdiag
import LightGraphs:
    src, dst, edgetype, nv, ne, vertices, edges, is_directed,
    add_vertex!, add_edge!, rem_vertex!, rem_edge!,
    has_vertex, has_edge, inneighbors, outneighbors,
    indegree, outdegree, degree, has_self_loops, num_self_loops,
    add_vertex!, adjacency_matrix, weights,
    SimpleGraph, SimpleDiGraph, IsDirected,

    # operators
    complement

export AbstractSimpleValueGraph, SimpleValueGraph, SimpleValueDiGraph, SimpleValueEdge,
    get_edgeval, set_edgeval!,
    val,
    outedgevals, inedgevals, all_edgevals,
    default_edgeval, edgeval_type,
    map_edgevals! #, kruskal_mst_modified

# ===== AbstractSimpleValueGraph ==========

abstract type AbstractSimpleValueGraph{V<:Integer, E_VAL} <: AbstractGraph{V} end

const TupleOrNamedTuple = Union{Tuple, NamedTuple} # TODO maybe somewhere else?

const default_edgeval_type = Float64
default_edgeval(E_VAL)::E_VAL = oneunit(E_VAL)
default_edgeval(::Type{Nothing})::Nothing = nothing
default_edgeval(T::Type{<:TupleOrNamedTuple})::T = T( default_edgeval(U) for U in T.types )

default_zero_edgeval(::Type{Nothing}) = nothing
default_zero_edgeval(E_VAL) = zero(E_VAL)
default_zero_edgeval(E_VAL::Type{<:TupleOrNamedTuple}) = E_VAL( default_zero_edgeval(T) for T in E_VAL.types )

eltype(::AbstractSimpleValueGraph{V}) where {V} = V
edgetype(::AbstractSimpleValueGraph{V, E_VAL}) where {V, E_VAL} = SimpleValueEdge{V, E_VAL}
edgeval_type(g::AbstractSimpleValueGraph{V, E_VAL}) where {V, E_VAL} = E_VAL
edgeval_type(::Type{<:AbstractSimpleValueGraph{V, E_VAL}}) where {V, E_VAL} = E_VAL

# TODO edges with key
edges(g::AbstractSimpleValueGraph) = SimpleValueEdgeIter(g)

nv(g::AbstractSimpleValueGraph) = eltype(g)(length(g.fadjlist))
ne(g::AbstractSimpleValueGraph) = g.ne

zero(G::AbstractSimpleValueGraph{V}) where {V} = G(zero(V))

vertices(g::AbstractSimpleValueGraph) = Base.OneTo(nv(g))
has_vertex(g::AbstractSimpleValueGraph, v) = v ∈ vertices(g)

struct SimpleValueEdgeIter{G<:AbstractSimpleValueGraph} <: AbstractEdgeIter
    g::G
end

length(iter::SimpleValueEdgeIter) = ne(iter.g)
eltype(::Type{<:SimpleValueEdgeIter{<:AbstractSimpleValueGraph{V, E_VAL}}}) where {V, E_VAL} =
        SimpleValueEdge{V, E_VAL}
eltype(edges::SimpleValueEdgeIter) = eltype(typeof(edges))

function show(io::IO, g::AbstractSimpleValueGraph)
    dir = is_directed(g) ? "directed" : "undirected"
    println(io, "{$(nv(g)), $(ne(g))} $dir $(eltype(g)) SimpleValueGraph with edge values of type $(edgeval_type(g))")
end

# ==== Includes ===========================

include("utils.jl")
include("simplevalueedge.jl")
include("simplevaluegraph.jl")
include("simplevaluedigraph.jl")
include("matrices.jl")

include("operators.jl")
# include("modified_functions.jl")

include("export.jl")

end # module


text:
module BinaryOptimization

export
    BinaryOptimizationModel,
    HOBO,
    Ising

include("binary_optimization_model.jl")
include("variables.jl")
include("ising_model.jl")
include("hobo_model.jl")
include("conversions.jl")
include("special_qubo.jl")
include("special_hobo.jl")
include("simplifiers.jl")
include("hamiltonians.jl")

#n = 10
#x = Variables("TSP", [10,10], ["time", "city"])
#getindex(x, 1, 2)

end # module


text:
function CompilerContext(
            model::PricingModel,
            provider::MarketData.AbstractMarketDataProvider,
            attr::ContractAttributes,
            pricing_date::Date,
            input_riskfactors_variable::OptimizingIR.ImmutableVariable,
            target_pricer_type::DataType
        ) :: CompilerContext

    ctx = CompilerContext(
            model,
            provider,
            attr,
            pricing_date,
            OptimizingIR.BasicBlock(),
            ObservablesBuffer(),
            input_riskfactors_variable,
            OptimizingIR.LookupTable{RiskFactor}(),
            target_pricer_type,
            target_pricer_type <: AbstractCashflowPricer ? Dict{OptimizingIR.ImmutableVariable, CashflowType}() : nothing
        )

    # creates the pricing function argument for the vector of risk factor values
    OptimizingIR.addinput!(ctx.program, ctx.input_riskfactors_variable)

    return ctx
end

function clone_context_with_model(ctx::CompilerContext, new_model::PricingModel)
    return CompilerContext(
            new_model,
            ctx.provider,
            ctx.attr,
            ctx.pricing_date,
            ctx.program,
            ctx.observables_buffer,
            ctx.input_riskfactors_variable,
            ctx.input_riskfactors,
            ctx.target_pricer_type,
            ctx.output_var_to_cashflowtype
        )
end

function compile_pricer(
        model::PricingModel,
        provider::MarketData.AbstractMarketDataProvider,
        attributes::ContractAttributes,
        pricing_date::Date,
        contract::Contract,
        target_pricer_type::Type{T}
        ;
        compiler::Symbol=:interpreter) where {T<:Union{AbstractPricer, AbstractCashflowPricer}}

    compiler_result = compile_pricing_program(model, provider, attributes, pricing_date, contract, target_pricer_type)
    return create_pricer(compiler_result, compiler=compiler)
end

function _compile_pricing_function(
            program::OptimizingIR.CompiledBasicBlock,
            compiler::Symbol;
            memory_buffer::Union{Nothing, Vector{Any}}=nothing,
            input_values_buffer::Union{Nothing, Vector{Any}}=nothing,
            auto_resize_buffers::Bool=true
        )

    if compiler === :interpreter

        if memory_buffer === nothing
            memory_buffer = Vector{Any}()
        end

        if input_values_buffer === nothing
            input_values_buffer = Vector{Any}()
        end

        return OptimizingIR.BasicBlockInterpreter(program, memory_buffer, input_values_buffer, auto_resize_buffers=auto_resize_buffers)

    elseif compiler === :native
        return OptimizingIR.compile(OptimizingIR.Native, program)
    else
        error("Unknown compiler type: $compiler")
    end
end

function create_pricer(
            compiler_result::CompilerResult{OptimizingIR.CompiledBasicBlock, P};
            compiler::Symbol=:interpreter,
            memory_buffer::Union{Nothing, Vector{Any}}=nothing,
            input_values_buffer::Union{Nothing, Vector{Any}}=nothing,
            auto_resize_buffers::Bool=true
        ) where {P<:AbstractPricer}

    f = _compile_pricing_function(compiler_result.program, compiler, memory_buffer=memory_buffer, input_values_buffer=input_values_buffer, auto_resize_buffers=auto_resize_buffers)

    @assert length(compiler_result.program.outputs) == 1 "pricing function should have exactly one output"

    return Pricer(
            f,
            compiler_result.input_riskfactors,
            compiler_result.currency
        )
end

function create_pricer(
            compiler_result::CompilerResult{OptimizingIR.CompiledBasicBlock, P};
            compiler::Symbol=:interpreter,
            memory_buffer::Union{Nothing, Vector{Any}}=nothing,
            input_values_buffer::Union{Nothing, Vector{Any}}=nothing,
            auto_resize_buffers::Bool=true
        ) where {P<:AbstractCashflowPricer}

    f = _compile_pricing_function(compiler_result.program, compiler, memory_buffer=memory_buffer, input_values_buffer=input_values_buffer, auto_resize_buffers=auto_resize_buffers)

    @assert length(compiler_result.program.outputs) >= 1 "pricing function should have at least one output."

    return CashflowPricer(
        f,
        compiler_result.input_riskfactors,
        compiler_result.currency,
        compiler_result.price_output_index,
        compiler_result.output_index_to_cashflow_type
    )
end

function init_compiler_context(
            model::PricingModel,
            provider::MarketData.AbstractMarketDataProvider,
            attributes::ContractAttributes,
            pricing_date::Date,
            target_pricer_type::Type{T}
        ) where {T<:Union{AbstractPricer, AbstractCashflowPricer}}

    return CompilerContext(model, provider, attributes, pricing_date, OptimizingIR.ImmutableVariable(:risk_factors_values), target_pricer_type)
end

function compile_pricing_program(
            model::PricingModel,
            provider::MarketData.AbstractMarketDataProvider,
            attributes::ContractAttributes,
            pricing_date::Date,
            contract::Contract,
            target_pricer_type::Type{T}
        ) where {T<:Union{AbstractPricer, AbstractCashflowPricer}}

    ctx = init_compiler_context(model, provider, attributes, pricing_date, target_pricer_type)
    price_value = lower!(ctx, contract)
    bind_output!(ctx, :price, price_value)

    return CompilerResult(ctx)
end

function CompilerResult(ctx::CompilerContext{M, D, IR, P}) where {M, D, IR<:OptimizingIR.BasicBlock, P<:AbstractPricer}
    # identifies the index of the return value with the pricing result for the contract
    price_output_index = OptimizingIR.indexof(ctx.program.outputs, OptimizingIR.ImmutableVariable(:price))

    # an AbstractPricer has a single return value with the price of the contract
    output_index_to_cashflow_type = nothing

    return CompilerResult(
            OptimizingIR.CompiledBasicBlock(ctx.program),
            ctx.input_riskfactors,
            get_functional_currency(ctx.model),
            ctx.target_pricer_type,
            price_output_index,
            output_index_to_cashflow_type
        )
end

function CompilerResult(ctx::CompilerContext{M, D, IR, P}) where {M, D, IR<:OptimizingIR.BasicBlock, P<:AbstractCashflowPricer}
    # identifies the index of the return value with the pricing result for the contract
    price_output_index = OptimizingIR.indexof(ctx.program.outputs, OptimizingIR.ImmutableVariable(:price))

    # an AbstractCashflowPricer has a return value for each cashflow
    @assert ctx.output_var_to_cashflowtype !== nothing
    output_index_to_cashflow_type = Dict{Int, CashflowType}()
    for (output_variable, cftype) in ctx.output_var_to_cashflowtype
        output_index_to_cashflow_type[OptimizingIR.indexof(ctx.program.outputs, output_variable)] = cftype
    end

    # outputs: pricing result + one result for each cashflow
    @assert length(output_index_to_cashflow_type) + 1 == length(ctx.program.outputs) "Some output value was not considered in the pricing routine."

    return CompilerResult(
            OptimizingIR.CompiledBasicBlock(ctx.program),
            ctx.input_riskfactors,
            get_functional_currency(ctx.model),
            ctx.target_pricer_type,
            price_output_index,
            output_index_to_cashflow_type
        )
end

#
# Register inputs and outputs
#

function add_input_riskfactor!(ctx::CompilerContext, rf::RiskFactor) :: OptimizingIR.ImmutableValue
    rf_index = OptimizingIR.addentry!(ctx.input_riskfactors, rf)
    return OptimizingIR.addinstruction!(ctx.program, OptimizingIR.call(OP_GETINDEX, ctx.input_riskfactors_variable, OptimizingIR.constant(rf_index)))
end

function bind_local_immutable_variable!(ctx::CompilerContext, var_sym::Symbol, var_value::OptimizingIR.ImmutableValue) :: OptimizingIR.ImmutableVariable
    var = OptimizingIR.ImmutableVariable(var_sym)
    OptimizingIR.assign!(ctx.program, var, var_value)
    return var
end

function bind_output!(ctx::CompilerContext, var_sym::Symbol, var_value::OptimizingIR.ImmutableValue)
    var = OptimizingIR.ImmutableVariable(var_sym)
    bind_output!(ctx, var, var_value)
    nothing
end

function bind_output!(ctx::CompilerContext, var::OptimizingIR.ImmutableVariable, var_value::OptimizingIR.ImmutableValue)
    OptimizingIR.addoutput!(ctx.program, var)
    OptimizingIR.assign!(ctx.program, var, var_value)
    nothing
end

@inline function add_instruction!(ctx::CompilerContext, op::OptimizingIR.Op, arg...) :: OptimizingIR.ImmutableValue
    return OptimizingIR.addinstruction!(ctx.program, OptimizingIR.call(op, arg...))
end


text:
module Maps

using Test
using Gridap
using Gridap.Helpers

export Map
export evaluate!
export evaluate
export return_size
export test_map

"""
Abstract map that takes an `M`-dim array of `S` values and returns an `N`-dim
array of `T` values.
For efficiency reasons, `T` should be a concrete type. However, `S` can be also an abstract type or
a union of types
"""
abstract type Map{S,M,T,N} end

"""
Evaluate a `Map` on a set of points
"""
function evaluate!(
  this::Map{S,M,T,N},
  points::AbstractArray{<:S,M},
  v::AbstractArray{T,N}) where {S,M,T,N}
  @abstractmethod
end

"""
Return dimension of the output array
"""
function return_size(
  ::Map{S,M,T,N},::NTuple{M,Int})::NTuple{N,Int} where {S,M,T,N}
  @abstractmethod
end

"""
Same as `evaluate!` but allocates output
"""
function evaluate(
  this::Map{S,M,T,N},
  points::AbstractArray{<:S,M}) where {S,M,T,N}
  v_size = return_size(this, size(points))
  v = Array{T,N}(undef, v_size)
  evaluate!(this,points,v)
  return v
end

# Testers

function test_map(
  m::Map{S,M,T,N},x::AbstractArray{<:S,M},y::AbstractArray{T,N}) where {S,M,T,N}

  z = evaluate(m,x)

  @test z ≈ y

  s = return_size(m,size(x))

  @test s == size(y)

end

end # module Maps


text:
# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
using HDF5
using Plots
using Plots.PlotMeasures
using GR
using QuadGK
using Printf
using LaTeXStrings
using DelimitedFiles
ENV["GKSwstype"] = "100"
# %%
H = 0.5 # height from galactic plane in kpc
rlim = 10 # limit of radius
kmsinkpcyr = 3.1536e7/3.085677581e16 # 1 km/sec in kpc/yr
# timestep = np.linspace(0.01,0.8,80)
time = 100 # snapshot
binwidth = 100
binwidthtemp = 0.25
patchsize = 0.5
Range = Int(rlim/patchsize)
rootdirectory = "/home/oku/SimulationData/isogal"

models = [
    # ["Osaka2019_isogal", "Osaka2019 (Shimizu+19)"],
        #   , "geodome_model/geodome_original"\
    # ["geodome_model/ver_19.11.1", "Geodesic dome model & Cioffi+ 1988"],
    # ["geodome_model/OKU2020","Geodesic dome model & Athena fitting"],
    # ["centroid_model/ver07271_NoMomCeiling","Centroid model & Athena fitting (alpha = 0)"],
    # ["centroid_model/ver07272_nfb1","Centroid model & Athena fitting (nfb = 1)"],
    # ["centroid_model/ver07272_SFE001","Centroid model & Athena fitting (SFE = 0.01)"],
    # ["centroid_model/ver07311","Centroid model & Athena fitting (alpha = -1)"],
    # ["centroid_model/ver07311_fdens-2","Centroid model & Athena fitting (alpha = -2)"],
    # ["centroid_model/ver08041_alpha-1","Centroid model & Athena fitting (new, alpha = -1)"],
    # ["centroid_model/ver08041_alpha-2","Centroid model & Athena fitting (new, alpha = -2)"],
    # ["centroid_model/ver08041_NoThermal","Centroid model (alpha = -1)"],
    # ["centroid_model/ver08191","Centroid model (alpha = -0.24)"],
    # ["centroid_model/ver07272_CHEVALIER1974","Centroid model & Cioffi+ 1988"],
    # ["centroid_model/ver09011_n1f2a0","Centroid model"],
    # ["centroid_model/ver09011_n1f2a024","Centroid model (new, alpha = -0.24)"],
    # ["centroid_model/ver09011_n1f2a050","Centroid model (new, alpha = -0.50)"],
    # ["centroid_model/ver09011_n1f2a075","Centroid model (new, alpha = -0.75)"],
    # ["centroid_model/ver09011_n1f2a1","Centroid model (new, alpha = -1.00)"],
    # ["centroid_model/ver12151","Centroid"],
    # ["ss_model/ver01092", "Spherical superbubble"],
    # ["ss_model/ver01121", "Spherical superbubble"],
    # ["ss_model/ver01193", "Spherical superbubble"],
    # ["ss_model/ver03221", "Spherical superbubble"],
    # ["ss_model/ver12211", "Spherical superbubble"],
    # ["ss_model/ver03311", "with cool-off"],
    # ["ss_model/ver03311/med", "high-reso"],
    # ["ss_model/ver03311_NoCoolOff", "without cool-off"],
    # ["NoFB", "No feedback"],
    ["ss_model/ver05291/fiducial", "high-reso"],
    ]
# %%
function W3(r, h)
  r = abs(r)/h
  C = 8/h^3/pi
  if r > 1
      return 0
  elseif r > 1/2
      return C*2*(1-r)^3
  else
      return C*(1 - 6*r^2 + 6*r^3)
  end
end
function crosssection(hsml, z)
    function integrand(x)
        return W3(sqrt(z^2 + x^2),hsml)*2*pi*x
    end
    return quadgk(integrand, 0, sqrt(hsml^2 - z^2))[1]
end
function getcenter(sph)
  x = y = z = mtot = 0
  for i in 1:length(sph.mass)
    x += sph.mass[i] * sph.pos[1, i]
    y += sph.mass[i] * sph.pos[2, i]
    z += sph.mass[i] * sph.pos[3, i]
    mtot += sph.mass[i]
  end
  x /= mtot
  y /= mtot
  z /= mtot
  return [x, y, z]
end

function main(H, sph, galaxy)
    SFR = [0.0 for i in 1:4*Range*Range]
    Gas = [0.0 for i in 1:4*Range*Range]
    MassOutflowHot = [0.0 for i in 1:4*Range*Range]
    MassOutflowCold = [0.0 for i in 1:4*Range*Range]
    EnergyOutflowHot = [0.0 for i in 1:4*Range*Range]
    EnergyOutflowCold = [0.0 for i in 1:4*Range*Range]
    Count = [0.0 for i in 1:4*Range*Range]
    Rho = [0.0 for i in 1:4*Range*Range]
    for i in 1:length(sph.mass)
      z = abs(sph.pos[3,i] - galaxy.pos[3,1])
      x = sph.pos[1,i] - galaxy.pos[1,1]
      y = sph.pos[2,i] - galaxy.pos[2,1]
      dz = abs(z - H)
      if abs(x) < rlim && abs(y) < rlim
        index = Int(ceil(x/patchsize)+Range + (ceil(y/patchsize)+Range)*Range)
        area = patchsize*patchsize
        SFR[index] += sph.sfr[i]/area
        Gas[index] += sph.mass[i]/area
        Rho[index] += sph.rho[i]
        Count[index] += 1
        if dz < sph.hsml[i]
          v = abs(sph.vel[3,i] - galaxy.vel[3,1])
          massoutflow = sph.mass[i]*2e43*crosssection(sph.hsml[i], dz)*v*kmsinkpcyr
          Energy = sph.mass[i]*(0.5*v^2 + sph.u[i])
          energyoutflow = Energy*2e53*crosssection(sph.hsml[i], dz)*v*kmsinkpcyr
          if sph.temp[i] > 1e5
            MassOutflowHot[index] += massoutflow/area
            EnergyOutflowHot[index] += energyoutflow/area
          else
            MassOutflowCold[index] += massoutflow/area
            EnergyOutflowCold[index] += energyoutflow/area
          end
        end
      end
    end
    Rho = Rho./Count
  return (SFR, MassOutflowHot, MassOutflowCold, EnergyOutflowHot, EnergyOutflowCold, Gas, Rho)
end
struct Sph
    pos::Array{Float32, 2}
    vel::Array{Float32, 2}
    hsml::Array{Float32, 1}
    mass::Array{Float32, 1}
    rho::Array{Float32, 1}
    temp::Array{Float32, 1}
    u::Array{Float32, 1}
    sfr::Array{Float32, 1}
end
struct Galaxy
    pos::Array{Float32, 1}
    vel::Array{Float32, 1}
    sfr::Float32
    stellarmass::Float32
end
# %%
SigmaSFR = [[] for i in 1:length(models)]
SigmaGas = [[] for i in 1:length(models)]
EnergyOutflowRate = [[[],[]] for i in 1:length(models)]
MassOutflowRate = [[[],[]] for i in 1:length(models)]
EnergyLoadingFactor = [[[],[]] for i in 1:length(models)]
SpecificEnergy = [[[],[]] for i in 1:length(models)]
Density = [[] for i in 1:length(models)]
DiskHeight = [[] for i in 1:length(models)]
for i in 1:length(models)
  j = time
        path2snapshot = @sprintf("%s/%s/snapshot_%03d/snapshot_%03d.hdf5", rootdirectory, models[i][1], j, j)
        path2subfind = @sprintf("%s/%s/snapshot_%03d/groups_%03d/sub_%03d.hdf5", rootdirectory, models[i][1], j, j, j)
        h5open(path2snapshot, "r") do fp
            pos = read(fp,"PartType0/Coordinates")
            vel = read(fp,"PartType0/Velocities")
            hsml = read(fp, "PartType0/SmoothingLength")
            mass = read(fp, "PartType0/Masses")
            rho = read(fp, "PartType0/Density")
            temp = read(fp, "PartType0/Temperature")
            u = read(fp, "PartType0/InternalEnergy")
            global sfr = read(fp, "PartType0/StarFormationRate")
            global massstars = read(fp, "PartType4/Masses")
            global sph = Sph(pos, vel, hsml, mass, rho, temp, u, sfr)
        end
        galpos = getcenter(sph)
        galvel = [0, 0, 0]
        galsfr = sum(sfr)
        stellarmass = sum(massstars)
        global galaxy = Galaxy(galpos, galvel, galsfr, stellarmass)

      # h5open(path2subfind, "r") do gp
      #   galpos = read(gp,"Group/GroupPos")
      #   galvel = read(gp,"Subhalo/SubhaloVel")
      #   galsfr = sum(sfr)
      #   stellarmass = sum(massstars)
      #   global galaxy = Galaxy(galpos, galvel, galsfr, stellarmass)
      # end
        results = main(H, sph, galaxy)
        SigmaSFR[i] = results[1]
        SigmaGas[i] = results[6]
        Density[i] = results[7]
        MassOutflowRate[i][1] = results[2]
        MassOutflowRate[i][2] = results[3]
        EnergyOutflowRate[i][1] = results[4]
        EnergyOutflowRate[i][2] = results[5]
        EnergyLoadingFactor[i][1] = EnergyOutflowRate[i][1]./SigmaSFR[i]/1e49
        EnergyLoadingFactor[i][2] = EnergyOutflowRate[i][2]./SigmaSFR[i]/1e49
        SpecificEnergy[i][1] = EnergyOutflowRate[i][1]./MassOutflowRate[i][1]
        SpecificEnergy[i][2] = EnergyOutflowRate[i][2]./MassOutflowRate[i][2]
        DiskHeight[i] = SigmaGas[i]./Density[i]
end
# %%
eta_hot = readdlm("eta_hot.txt")
eta_cold = readdlm("eta_cold.txt")
# %%

using Plots
gr()
font = Plots.font(16)
lfont = Plots.font(12)
Plots.plot(legend=:bottomleft,xlim=(1e-4, 1), ylim=(1e-7, 10), xscale=:log10, yscale=:log10, xlabel=latexstring("\$\\Sigma_{\\mathrm{SFR}} [M_{\\odot} \\mathrm{kpc}^{-2} \\mathrm{yr}^{-1}]\$"), ylabel=L"\eta_{\rm E}", xtickfont=font, ytickfont=font, xlabelfont=font, ylabelfont=font, legendfont=lfont, guidefont=font, topmargin=3mm, bottommargin=3mm, rightmargin=3mm)
for i in 1:length(models)
Plots.plot!(SigmaSFR[i], EnergyLoadingFactor[i][1], label="hot outflow (This work)", markersize=5, markercolor=:red, st=:scatter)
Plots.plot!(SigmaSFR[i], EnergyLoadingFactor[i][2], label="cold outflow (This work)", markersize=5, markercolor=:blue, st=:scatter)
Plots.scatter!(eta_hot[:,1], eta_hot[:,2], label="hot outflow (Li & Bryan 20)", markersize=5, markercolor=:red, marker=:xcross)
Plots.scatter!(eta_cold[:,1], eta_cold[:,2], label="cold outflow (Li & Bryan 20)", markersize=5, markercolor=:blue, marker=:xcross)
end
# Plots.savefig("Energyloading.png")
Plots.savefig("results/Energyloading.pdf")


text:
using PowerModelsReliability
using PowerModels
using InfrastructureModels
using Ipopt
using Mosek
using Juniper
using Cbc
using CPLEX

mosek = MosekSolver()
ipopt = IpoptSolver()
juniper = JuniperSolver(IpoptSolver(print_level=0); mip_solver=CplexSolver())

data_tf = PowerModels.parse_file("./test/data/case5_tf.m")
data_tf["multinetwork"] = false
resultACtf = run_tfopf(data_tf, ACPPowerModel, ipopt; multinetwork=false, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => true))
data_tf = PowerModels.parse_file("./test/data/case5_tf.m")
data_tf["multinetwork"] = false
resultDCtf = run_tfopf(data_tf, DCPPowerModel, ipopt; multinetwork=false, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => true))

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
include("../src/numint_module.jl")

module ntest
using ..numint
using BenchmarkTools

# fcn(x) = x^3
# fcn(x) = x^4 * log(x + sqrt(x^2 + 1))
fcn(x) = exp(-3*x) - cos(5*pi*x)

interval = numint.Interval(0.0, 8.0)

N = 10000
# I = numint.trapezoidal(fcn, interval, N)
bt = @benchmark I = numint.trapezoidal(fcn, interval, N)
show(stdout, MIME("text/plain"), bt)

integration_options = numint.RombergOptions(abs_tol = 1e-8, rel_tol = 1e-8, N = 250)
br = @benchmark stat_romb = numint.romberg($fcn, $interval, $integration_options)
show(stdout, MIME("text/plain"), br)

bgk = @benchmark stat_gk = numint.gauss_kronrod($fcn, $interval)
show(stdout, MIME("text/plain"), bgk)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
mutable struct TreeSet{LABEL, RT, NL, N, B, TREE <:
                       AbstractTree{OneTree, RT, NL, N, B}} <:
    AbstractTree{ManyTrees, RT, NL, N, B}
    trees::Dict{LABEL, TREE}
    treeinfo::Dict{LABEL, Dict{String, Any}}

    TreeSet(treedict::Dict{LABEL, TREE},
            info::Dict{LABEL, Dict{String, Any}} =
            Dict(broadcast(key -> (key, Dict{String, Any}()), keys(treedict)))) where
    {LABEL, RT, NL, N, B, TREE <: AbstractTree{OneTree, RT, NL, N, B}} =
        new{LABEL, RT, NL, N, B, TREE}(treedict, info)
end

TreeSet(trees::AbstractVector{T}) where T <: AbstractTree{OneTree} =
    TreeSet(Dict(Pair.(Base.OneTo(length(trees)), trees)))

import Base.IteratorSize
function IteratorSize(::Type{TreeSet})
    return HasLength()
end

import Base.IteratorEltype
function IteratorEltype(::Type{TreeSet})
    return HasEltype()
end

import Phylo.API: _getleafinfo
_getleafinfo(ts::TreeSet) = _getleafinfo(first(gettrees(ts)))
_getleafinfo(ts::TreeSet, leafname) = _getleafinfo(first(gettrees(ts)), leafname)

import Phylo.API: _ntrees
_ntrees(ts::TreeSet) = length(ts.trees)

import Phylo.API: _gettrees
_gettrees(ts::TreeSet) = values(ts.trees)

import Phylo.API: _gettreenames
_gettreenames(ts::TreeSet) = keys(ts.trees)

import Phylo.API: _gettreeinfo
_gettreeinfo(ts::TreeSet) = ts.treeinfo
_gettreeinfo(ts::TreeSet{LABEL}, name::LABEL) where LABEL = ts.treeinfo[name]

import Phylo.API: _gettree
_gettree(ts::TreeSet, name) = ts.trees[name]

import Base.eltype
eltype(::TreeSet{LABEL, RT, NL, N, B, TREE}) where
    {LABEL, RT, NL, N, B,
     TREE <: AbstractTree{OneTree, RT, NL, N, B}} = TREE

import Base.length
length(ts::TreeSet) = length(ts.trees)

import Base.getindex
getindex(ts::TreeSet{LABEL}, idx::LABEL) where LABEL = ts.trees[idx]

import Phylo.API: _getleafnames
function _getleafnames(ts::TREESET) where {LABEL, RT, NL, N, B,
                    TREE <: AbstractTree{OneTree, RT, NL, N, B},
                    TREESET <: TreeSet{LABEL, RT, NL, N, B, TREE}}
    lns = unique(_getleafnames(t) for t in values(ts.trees))
    if length(lns) > 1
        error("Inconsistent leaf names in TreeSet")
    elseif isempty(lns)
        return NL[]
    end
    return first(lns)
end

import Phylo.API: _nleaves
function _nleaves(ts::TREESET) where {LABEL, RT, NL, N, B,
                    TREE <: AbstractTree{OneTree, RT, NL, N, B},
                    TREESET <: TreeSet{LABEL, RT, NL, N, B, TREE}}
    nls = unique(_nleaves(t) for t in values(ts.trees))
    if length(nls) > 1
        error("Inconsistent leafinfo in TreeSet")
    elseif isempty(nls)
        return 0
    end
    return first(nls)
end

import Base.keys
keys(ts::TreeSet) = keys(ts.trees)


text:
module ToStruct

function tostruct(T::Type, x::Any)
    try
        x::T
    catch
        try
            convert(T, x)
        catch
            T(x)
        end
    end
end

function tostruct(T::Union, x::Any)
    # NOTE: T.b is more specific in most cases
    #       For instance union of String and Nothing types wile automatically ordered to be Union{Nothing,String}
    try
        tostruct(T.b, x)
    catch
        tostruct(T.a, x)
    end
end

function tostruct(T::Type{U} where U<:AbstractVector, x::AbstractVector)
    ET = eltype(T)
    T(collect(tostruct(ET, e) for e in x))
end

function tostruct(T::Type{U} where U<:AbstractDict, x::AbstractDict)
    KT, VT = eltype(T()).types
    T(tostruct(KT, k) => tostruct(VT, v) for (k, v) in x)
end

function tostruct(T::DataType, x::AbstractDict{AbstractString})
    args = map(fieldnames(T)) do fname
        FT = fieldtype(T, fname)
        v = getdefault(FT, x, String(fname))
        tostruct(FT, v)
    end
    T(args...)
end

tostruct(T::DataType, x::AbstractDict) = tostruct(T, convert(Dict{AbstractString,Any}, x))

function getdefault(T::Type, x::AbstractDict, k::Any)
    if T >: Nothing
        get(x, k, nothing)
    elseif T >: Missing
        get(x, k, missing)
    else
        x[k]
    end
end

end # module


text:
# This file was generated, do not modify it. # hide
fig = figure(figsize=(15, 15))

subplot(221)
scatter(data.Wind_gen, data.Radiation_dir)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Solar radiation - directional", fontsize=14)

subplot(222)
scatter(data.Wind_gen, data.Radiation_dif)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Solar radiation - diffuse", fontsize=14)

subplot(223)
scatter(data.Wind_gen, data.Windspeed)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Wind speed (m/s)", fontsize=14)

subplot(224)
scatter(data.Wind_gen, data.Temperature)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Temperature (C)", fontsize=14)

savefig(joinpath(@OUTPUT, "wind_scatter.png"), bbox_inches="tight") # hide


text:
module JPLEphemeris

using Reexport

@reexport using AstroBase.Bodies

include("daf.jl")
include("spk.jl")

end


text:
#################################################
# Rendering
# This controls how failures are displayed
abstract type RenderMode end
struct Diff <: RenderMode end

abstract type BeforeAfter <: RenderMode end
struct BeforeAfterLimited <: BeforeAfter end
struct BeforeAfterFull <: BeforeAfter end
struct BeforeAfterImage <: BeforeAfter end

render_item(::RenderMode, item) = println(item)
function render_item(::BeforeAfterLimited, item)
    show(IOContext(stdout, :limit=>true, :displaysize=>(20,80)), "text/plain", item)
    println()
end
function render_item(::BeforeAfterImage, item)
    str_item = @withcolor ImageInTerminal.encodeimg(ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), item, 20, 40)[1]
    println("eltype: ", eltype(item))
    println("size: ", map(length, axes(item)))
    println("thumbnail:")
    println.(str_item)
end

## 2 arg form render for comparing
function render(mode::BeforeAfter, reference, actual)
    println("- REFERENCE -------------------")
    render_item(mode, reference)
    println("-------------------------------")
    println("- ACTUAL ----------------------")
    render_item(mode, actual)
    println("-------------------------------")
end
function render(::Diff, reference, actual)
    println("- DIFF ------------------------")
    @withcolor println(deepdiff(reference, actual))
    println("-------------------------------")
end

## 1 arg form render for new content
function render(mode::RenderMode, actual)
    println("- NEW CONTENT -----------------")
    render_item(mode, actual)
    println("-------------------------------")
end

"""
    default_rendermode(::DataFormat, actual)

Infer the most appropriate render mode according to type of reference file and `actual`.
"""
default_rendermode(::Type{<:DataFormat}, ::Any) = BeforeAfterLimited()
default_rendermode(::Type{<:DataFormat}, ::AbstractString) = Diff()
default_rendermode(::Type{<:DataFormat}, ::AbstractArray{<:Colorant}) = BeforeAfterImage()

# plain TXTs
default_rendermode(::Type{DataFormat{:TXT}}, ::Any) = Diff()
default_rendermode(::Type{DataFormat{:TXT}}, ::AbstractString) = Diff()
default_rendermode(::Type{DataFormat{:TXT}}, ::Number) = BeforeAfterFull()
default_rendermode(::Type{DataFormat{:TXT}}, ::AbstractArray{<:Colorant}) = BeforeAfterImage()

# SHA256
default_rendermode(::Type{DataFormat{:SHA256}}, ::Any) = BeforeAfterFull()
default_rendermode(::Type{DataFormat{:SHA256}}, ::AbstractString) = BeforeAfterFull()
default_rendermode(::Type{DataFormat{:SHA256}}, ::AbstractArray{<:Colorant}) = BeforeAfterLimited()


text:
module day7

using ..aoc19: run_program!

using DelimitedFiles, Combinatorics

readInput() = readdlm(joinpath(@__DIR__, "input.txt"), ',', Int) |> vec

part1(input) = maximum(run_amplifiers(input, p) for p in permutations(0:4))

part2(input) = maximum(run_amplifiers(input, p) for p in permutations(5:9))

function run_amplifiers(input, phases)
    channels = [Channel{Int}(2) for i = 1:5]
    put!.(channels, phases) #phase setting
    put!(first(channels), 0) #initial input
    tasks = [@async run_program!(copy(input), channels[i], channels[mod1(i+1, 5)]) for i = 1:5]
    wait.(tasks)
    take!(first(channels))
end

end


text:
# __precompile__(false)
module GraphFusedElasticNet

include("Utils.jl")
include("L1Solver.jl")
include("L2Solver.jl")
include("BinomialGFEN.jl")
include("GaussianGFEN.jl")
include("DensityTree.jl")
include("GraphTrails.jl")
include("BayesianOpt.jl")
include("ARS.jl")
include("BayesianBinomialGFEN.jl")

export BinomialGFEN, GaussianGFEN
export BayesianBinomialGFEN
export sample_chain
export filter1D!, filter1Dl2!, filter1D, filter1Dl2
export fit!, predict
export Trails
export graph_from_edgelist, find_trails, grid_trails
export GaussianProcessSampler, RandomGaussianProcessSampler
export addobs!, gpsample, gpeval
export eval_logprob, make_tree_from_bfs, uniform_binary_splits, eval_logdens, make_tree

end


text:
@testset "FieldElement Operations" begin
    import Secp256k1: 𝐹, P
    @testset "Addition" begin
        @test 𝐹(2)   + 𝐹(15) == 𝐹(17)
        @test 𝐹(P-1) + 𝐹(2)  == 𝐹(1)
    end
    @testset "Substraction" begin
        @test 𝐹(29) - 𝐹(4) == 𝐹(25)
        @test 𝐹(1)  - 𝐹(2) == 𝐹(P-1)
    end
    @testset "Multiplication" begin
        a = 𝐹(2)
        b = 𝐹(P-1)
        @test a * b == 𝐹(P-2)
        @test a * 2 == 𝐹(4)
        @test 3 * b == 𝐹(P-3)
    end
    @testset "Power" begin
        @test 𝐹(2)^3 == 𝐹(8)
        @test 𝐹(P-1)^5 * 𝐹(18) == 𝐹(P-18)
    end
    @testset "Division" begin
        a = 𝐹(P-3)
        b = 𝐹(3)
        @test 𝐹(P-3) / 𝐹(3)    == 𝐹(P-1)
        @test 𝐹(1)^-rand(UInt) == 𝐹(1)
        @test 𝐹(4)^-2 == 𝐹(big"0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc3")
    end
end


text:
struct Color
    name::String
end

Color(::Nothing) = nothing
Color(c::Color) = c


text:
using LightGraphs, Gadfly, GraphPlot, LinearAlgebra, SimpleWeightedGraphs;

# Test network
N = 8
    Pij = zeros(N,N)
    Pij[1,2] = Pij[2,1] = 1
    Pij[3,2] = Pij[2,3] = 1
    Pij[3,4] = Pij[4,3] = 1
    Pij[1,5] = Pij[5,1] = 1
    Pij[6,5] = Pij[5,6] = 1
    Pij[6,7] = Pij[7,6] = 1
    Pij[6,8] = Pij[8,6] = 1

    display(Pij)

# Creating ancestor vector
    A_i = zeros(N)
    for i in 2:N
        for j in 1:i-1
            if Pij[i,j] ≠ 0
                A_i[i] = j
            end
        end
    end
    println("A_i = "*string(A_i))

    locs_x, locs_y = collect(1.0:N), zeros(N);

    for n in 2:N
        if A_i[n] ≠ n-1
            locs_x[n] = locs_x[convert(Int,A_i[n])] + 1
            locs_y[n] = locs_y[convert(Int,A_i[n])] - 0.5
        else
            locs_x[n] = locs_x[n-1]+1
            locs_y[n] = locs_y[n-1]
        end
    end

    n_label = collect(1:N);
    n = 2
    while n ≤ N
    global n, N, n_label, locs_x, locs_y
    println("n=$n, N=$N")
        if A_i[n] ≠ n-1
            println("Branching @ n = $n")
            # Adding phantom node on x+0.5 and same y
            insert!(locs_x,n,0.5+locs_x[convert(Int,A_i[n])])
            insert!(locs_y,n,locs_y[convert(Int,A_i[n])])
            insert!(A_i,    n,A_i[n])
            insert!(n_label,n,0)

            # Adding phantom node at x and y-0.5
            insert!(locs_x,n+1,locs_x[n])
            insert!(locs_y,n+1,locs_y[n]-0.5)
            insert!(A_i,    n+1,n)
            insert!(n_label,n+1,0)

            n = n+2

            # Adjusting the parent of the branching node
            A_i[n] = n-1

            N=N+2
            for j in n+1:N
                A_i[j] = A_i[j]+2
            end
            display([A_i n_label])
        end
        n += 1
        println("New n: $n")
    end

    g = SimpleGraph(N);
        NodeSize = ones(N)
        for i = 2:N
            add_edge!(g, i, convert(Int,A_i[i]))
            if n_label[i]==0
            NodeSize[i] = 0.2
            end
        end

    L = length(A_i);

        linw = collect(1:L);

        edgelinewidth = 40.0;
        nodelabeldist = 2;
        nodelabelangleoffset = π/4;
        NODELABELSIZE = 6;
        nodefillc = "black";
        edgestrokec = "black";
        EDGELINEWIDTH = 5/sqrt(N);
        # ==================================================
        # Function
        gplot(g, Float32.(locs_x), Float32.(locs_y),
        nodelabel = n_label,
        nodelabeldist = nodelabeldist,
        nodelabelangleoffset = nodelabelangleoffset,
        NODELABELSIZE = NODELABELSIZE,
        nodelabelsize = NodeSize,
        nodesize = NodeSize,
        nodefillc = nodefillc,
        edgestrokec = edgestrokec,
        EDGELINEWIDTH = EDGELINEWIDTH,
        );


text:
let
    function f_lm(x)
      [x[1], 2.0 - x[2]]
    end
    function g_lm(x)
      [1.0 0.0; 0.0 -1.0]
    end

    initial_x = [100.0, 100.0]

    results = Optim.levenberg_marquardt(f_lm, g_lm, initial_x)
    @assert norm(Optim.minimizer(results) - [0.0, 2.0]) < 0.01

    function rosenbrock_res(x, r)
        r[1] = 10.0 * (x[2] - x[1]^2 )
        r[2] =  1.0 - x[1]
        return r
    end

    function rosenbrock_jac(x, j)
        j[1, 1] = -20.0 * x[1]
        j[1, 2] =  10.0
        j[2, 1] =  -1.0
        j[2, 2] =   0.0
        return j
    end

    r = zeros(2)
    j = zeros(2,2)

    frb(x) = rosenbrock_res(x, r)
    grb(x) = rosenbrock_jac(x, j)

    initial_xrb = [-1.2, 1.0]

    results = Optim.levenberg_marquardt(frb, grb, initial_xrb)

    @assert norm(Optim.minimizer(results) - [1.0, 1.0]) < 0.01

    # check estimate is within the bound PR #278
     result = Optim.levenberg_marquardt(frb, grb, [150.0, 150.0]; lower = [10.0, 10.0], upper = [200.0, 200.0])
     @test Optim.minimizer(result)[1] >= 10.0
     @test Optim.minimizer(result)[2] >= 10.0

    # tests for #178, taken from LsqFit.jl, but stripped
    let
        srand(12345)

        model(x, p) = p[1]*map(exp,-x.*p[2])

        xdata = linspace(0,10,20)
        ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))

        f_lsq = p -> p[1]*map(exp, -xdata.*p[2])-ydata
        g_lsq = Calculus.jacobian(f_lsq)
        results = Optim.levenberg_marquardt(f_lsq, g_lsq, [0.5, 0.5])

        @assert norm(Optim.minimizer(results) - [1.0, 2.0]) < 0.05
    end

    let
        srand(12345)

        model(x, p) = p[1]*map(exp, -x./p[2])+p[3]

        xdata = 1:100
        ydata = model(xdata, [10.0, 10.0, 10.0]) + 0.1*randn(length(xdata))

        f_lsq = p -> model(xdata,p)-ydata
        g_lsq = Calculus.jacobian(f_lsq)

        # tests for box constraints, PR #196
        @test_throws ArgumentError Optim.levenberg_marquardt(f_lsq, g_lsq, [15.0, 15.0, 15.0], lower=[5.0, 11.0])
        @test_throws ArgumentError Optim.levenberg_marquardt(f_lsq, g_lsq, [5.0, 5.0, 5.0], upper=[15.0, 9.0])
        @test_throws ArgumentError Optim.levenberg_marquardt(f_lsq, g_lsq, [15.0, 10.0, 15.0], lower=[5.0, 11.0, 5.0])
        @test_throws ArgumentError Optim.levenberg_marquardt(f_lsq, g_lsq, [5.0, 10.0, 5.0], upper=[15.0, 9.0, 15.0])

        lower=[5.0, 11.0, 5.0]
        results = Optim.levenberg_marquardt(f_lsq, g_lsq, [15.0, 15.0, 15.0], lower=lower)
        Optim.minimizer(results)
        @test Optim.converged(results)
        @test all(Optim.minimizer(results) .>= lower)

        upper=[15.0, 9.0, 15.0]
        results = Optim.levenberg_marquardt(f_lsq, g_lsq, [5.0, 5.0, 5.0], upper=upper)
        Optim.minimizer(results)
        @test Optim.converged(results)
        @test all(Optim.minimizer(results) .<= upper)

        # tests for PR #267
        Optim.levenberg_marquardt(f_lsq, g_lsq, [15.0, 15.0, 15.0], show_trace=true)
    end
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Extract and organize json files into subdirectories according to decision id and
members from clustering results.
Use RLESCAS json_update_recur (if needed) and recursive_plot to generate the plots
recursively on the subdirectories.
Example:
using GrammarExpts,ExtractMembers
extract_members(...)
using RLESCAS
include_visualize()
json_update1_recur("Clusters") #if needed to update format
recursive_plot("Clusters")
"""
module ExtractMembers

export extract_members

using DataFrames
using RLESCAS #this is a big dependency, for is_nmac() only

const DATADIR = Pkg.dir("Datasets/data")
const METAFILE = "_META.csv.gz"

#TODO: don't access meta file directly
"""
N is number of members to extract for each node
"""
function extract_members(memberfile::AbstractString, jsondir::AbstractString,
    jsonroot::AbstractString, dataset::AbstractString, outdir::AbstractString,
    N::Int64; jsonext::AbstractString="json.gz")

    metafile = joinpath(DATADIR, dataset, METAFILE)
    meta = readtable(metafile)
    D = readtable(memberfile)
    for i = 1:nrow(D)
        id = D[i, 1]
        members = string2array(D[i, 2])
        subdir = joinpath(outdir, string(id))
        mkpath(subdir)
        NN = 0 #track number of members copied
        for m in members
            enc = meta[m, :encounter_id]
            fname = "$jsonroot$enc.$jsonext"
            srcpath = joinpath(jsondir, fname)
            if is_nmac(srcpath) #only copy nmacs
                cp(srcpath, joinpath(subdir, fname))
                NN += 1
                NN >= N && break #only process first N of each cluster
            end
        end
    end
end

function string2array(str_array::AbstractString; T::Type=Int64)
    str = replace(str_array, "[", "")
    str = replace(str, "]", "")
    A = split(str, ',')
    A = map(s->parse(T, s), A)
    A
end

end


text:
# Compare with tabulated values,  classical hypergeometric formulas

@testset "Bessel" begin
    T = Float64

    for α in (3 / 2, 1 / 2, 1, 2)
        P = Bessel{α,T}
        β = 2

        p0, p1, p2, p3, p4 = basis.(P, 0:4)

        # Krall and Frink
        for x in range(0, 1, length=5)
            @test p0(x) ≈ 1
            @test p1(x) ≈ 1 + α * (x / β)
            @test p2(x) ≈ 1 + 2(α + 1) * (x / β) + (α + 1) * (α + 2) * (x / β)^2
            @test p3(x) ≈
                  1 +
                  3(α + 2) * (x / β) +
                  3(α + 2) * (α + 3) * (x / β)^2 +
                  (α + 2) * (α + 3) * (α + 4) * (x / β)^3
            @test p4(x) ≈
                  1 +
                  4(α + 3) * (x / β) +
                  6(α + 3) * (α + 4) * (x / β)^2 +
                  4(α + 3) * (α + 4) * (α + 5) * (x / β)^3 +
                  (α + 3) * (α + 4) * (α + 5) * (α + 6) * (x / β)^4
        end
    end

    x = variable(Polynomial)
    for α in (1 / 4, 1 / 2, 3 / 4, 1, 2)
        P = Bessel{α}
        for n in 0:5
            p = basis(P, n)
            @test p(x) ≈ SP.classical_hypergeometric(P, n, x)
        end
    end
end

@testset "ChebyshevT" begin
    p0 = Chebyshev([1])  # 1
    p1 = Chebyshev([0, 1]) # x
    p2 = Chebyshev([0, 0, 1]) #
    p3 = Chebyshev([0, 0, 0, 1]) #
    p4 = Chebyshev([0, 0, 0, 0, 1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ x
        @test p2(x) ≈ 2x^2 - 1
        @test p3(x) ≈ 4x^3 - 3x
        @test p4(x) ≈ 8x^4 - 8x^2 + 1
    end
end

@testset "ChebyshevU" begin
    p0 = ChebyshevU([1])  # 1
    p1 = ChebyshevU([0, 1]) # x
    p2 = ChebyshevU([0, 0, 1]) #
    p3 = ChebyshevU([0, 0, 0, 1]) #
    p4 = ChebyshevU([0, 0, 0, 0, 1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ 2x
        @test p2(x) ≈ 4x^2 - 1
        @test p3(x) ≈ 8x^3 - 4x
        @test p4(x) ≈ 16x^4 - 12x^2 + 1
    end
end

@testset "Hermite" begin
    p0 = Hermite([1])  # 1
    p1 = Hermite([0, 1]) # x
    p2 = Hermite([0, 0, 1]) #
    p3 = Hermite([0, 0, 0, 1]) #
    p4 = Hermite([0, 0, 0, 0, 1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ 2x
        @test p2(x) ≈ 4x^2 - 2
        @test p3(x) ≈ 8x^3 - 12x
        @test p4(x) ≈ 16x^4 - 48x^2 + 12
    end

    x = variable(LaurentPolynomial)
    P = Hermite
    for n in 2:5
        p = basis(P, n)
        @test p(x) ≈ SP.classical_hypergeometric(P, n, x)
    end
end

@testset "ChebyshevHermite" begin
    p0 = ChebyshevHermite([1])  # 1
    p1 = ChebyshevHermite([0, 1]) # x
    p2 = ChebyshevHermite([0, 0, 1]) #
    p3 = ChebyshevHermite([0, 0, 0, 1]) #
    p4 = ChebyshevHermite([0, 0, 0, 0, 1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ x
        @test p2(x) ≈ x^2 - 1
        @test p3(x) ≈ x^3 - 3x
        @test p4(x) ≈ x^4 - 6x^2 + 3
    end

    x = variable(LaurentPolynomial)
    P = ChebyshevHermite
    for n in 2:5
        p = basis(P, n)
        @test p(x) ≈ SP.classical_hypergeometric(P, n, x) atol = 1e-8
    end
end

@testset "Legendre" begin
    p0 = Legendre([1])  # 1
    p1 = Legendre([0, 1]) # x
    p2 = Legendre([0, 0, 1]) #
    p3 = Legendre([0, 0, 0, 1]) #
    p4 = Legendre([0, 0, 0, 0, 1]) #

    for x in range(-1, stop=1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ x
        @test p2(x) ≈ 1 / 2 * (3x^2 - 1)
        @test p3(x) ≈ 1 / 2 * (5x^3 - 3x)
        @test p4(x) ≈ 1 / 8 * (35x^4 - 30x^2 + 3)
    end
end

@testset "Shifted Legendre" begin
    T = Float64
    P = ShiftedLegendre{T}

    p0 = P([1])  # 1
    p1 = P([0, 1]) # x
    p2 = P([0, 0, 1]) #
    p3 = P([0, 0, 0, 1]) #
    p4 = P([0, 0, 0, 0, 1]) #

    for x in range(0, stop=1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ 2x - 1
        @test p2(x) ≈ 6x^2 - 6x + 1
        @test p3(x) ≈ 20x^3 - 30x^2 + 12x - 1
        @test p4(x) ≈ 70x^4 - 140x^3 + 90x^2 - 20x + 1
    end
end

@testset "Laguerre" begin
    P = Laguerre{0}
    p0 = P([1])  # 1
    p1 = P([0, 1]) # x
    p2 = P([0, 0, 1]) #
    p3 = P([0, 0, 0, 1]) #
    p4 = P([0, 0, 0, 0, 1]) #

    for x in range(0, stop=5, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ -x + 1
        @test p2(x) ≈ 1 / 2 * (x^2 - 4x + 2)
        @test p3(x) ≈ 1 / 6 * (-x^3 + 9x^2 - 18x + 6)
        @test p4(x) ≈ 1 / 24 * (x^4 - 16x^3 + 72x^2 - 96x + 24)
    end

    x = variable(Polynomial)
    for α in (0, 1 / 4, 1 / 2, 3 / 4)
        P = Laguerre{α}
        for n in 2:5
            p = basis(P, n)
            @test p(x) ≈ SP.classical_hypergeometric(P, n, x)
        end
    end
end

@testset "Gegenbauer" begin
    for α in (1 / 4, 1 / 3, 1 / 2, 1, 2)
        P = Gegenbauer{α,Float64}
        p0 = P([1])  # 1
        p1 = P([0, 1]) # x
        p2 = P([0, 0, 1]) #
        p3 = P([0, 0, 0, 1]) #
        p4 = P([0, 0, 0, 0, 1]) #

        for x in range(-1, stop=1, length=5)
            @test p0(x) ≈ 1
            @test p1(x) ≈ 2α * x
            @test p2(x) ≈ -α + 2α * (1 + α) * x^2
            @test p3(x) ≈ -2α * (1 + α) * x + 4 / 3 * α * (1 + α) * (2 + α) * x^3
        end
    end

    x = variable(Polynomial)
    for α in (1 / 4, 1 / 2, 3 / 4, 1)
        P = Gegenbauer{α}
        for n in 2:5
            p = basis(P, n)
            @test p(x) ≈ SP.classical_hypergeometric(P, n, x)
        end
    end
end

@testset "Jacobi" begin
    T = Float64
    x = variable()

    # ChebyshevT = J(-1/2, -1/2)
    for (alpha_beta, Q) in (
        ((-1 / 2, -1 / 2), Chebyshev{T}),
        ((1 / 2, 1 / 2), ChebyshevU{T}),
        ((0, 0), Legendre{T}),
    )
        P = Jacobi{alpha_beta...,T}
        for i in 2:6
            p = basis(P, i)
            q = basis(Q, i)
            @test SP.monic(p)(x) ≈ SP.monic(q)(x)
        end
    end

    x = variable(Polynomial)
    for (α, β) in ((1 / 2, 1 / 2), (-1 / 2, 1 / 2), (1 / 2, -1 / 2), (-1 / 2, -1 / 2))
        P = Jacobi{α,β}
        for n in 2:5
            p = basis(P, n)
            @test p(x) ≈ SP.classical_hypergeometric(P, n, x)
        end
    end

    ## We have `jacobi_eval` for α and β outside [-1, ∞)
    α, β = 1 / 2, -1 / 2
    for n in 2:5
        x0 = 1 / 2
        @test SP.jacobi_eval(α, β, n, x0) ≈ basis(Jacobi{α,β}, n)(x0)
    end
end

## CDOP

@testset "Charlier" begin
    x = variable()
    for μ in (1 / 4, 1 / 2, 1, 2)
        P = Charlier{μ}
        @test basis(P, 0)(x) ≈ one(x)
        @test basis(P, 1)(x) ≈ 1 - x / μ
        @test basis(P, 2)(x) ≈ (x^2 + μ^2 - x * (1 + 2μ)) / μ^2
    end

    for μ in (1 / 4, 1 / 2, 1, 2)
        P = Charlier{μ}
        for i in 0:5
            @test basis(P, i)(x) ≈ SP.classical_hypergeometric(P, i, x)
        end
    end
end

@testset "Hahn" begin
    x = variable()
    for α in (1 / 4, 1 / 2)
        for β in (1 / 4, 1 / 2)
            𝐍 = 5

            P = Hahn{α,β,𝐍}
            for i in 0:𝐍
                @test basis(P, i)(x) ≈ SP.classical_hypergeometric(P, i, x)
            end
        end
    end
end

@testset "Krawchouk" begin
    x = variable()
    for p in (1 / 4, 1 / 2, 1, 2)
        for N in (3, 6)
            P = Krawchouk{p,N}
            @test basis(P, 0)(x) ≈ one(x)
            @test basis(P, 1)(x) ≈ -N * p + x
            @test basis(P, 2)(x) ≈ 1 / 2 * (N^2 * p^2 + x * (2p + x - 1) - N * p * (p + 2x))
        end
    end

    for p in (1 / 4, 1 / 2, 1, 2)
        for N in (3, 6)
            P = Krawchouk{p,N}
            for i in 0:N
                @test basis(P, i)(x) ≈ SP.classical_hypergeometric(P, i, x)
            end
        end
    end
end

@testset "Meixner" begin
    x = variable()
    for γ in (1 / 4, 1 / 2, 3 / 4)
        for μ in (1 / 4, 1 / 2, 3 / 4)
            P = Meixner{γ,μ}
            for i in 0:4
                @test basis(P, i)(x) ≈ SP.classical_hypergeometric(P, i, x)
            end
        end
    end
end


text:
function compareGeneratedCode(groundTruth::String, value::Expr)
	string(Meta.parse(groundTruth)) == join(filter(s -> !startswith(s, "#="), strip.(split(string(value), "\n")[2:end-1])), "\n");
end

function testThreadedForeach()::Bool
	return compareGeneratedCode("tforeach(identity, [])", @macroexpand @threaded foreach(identity, []));
end

function testThreadedMapfoldl()::Bool
	return compareGeneratedCode("tmapfoldl(identity, [])", @macroexpand @threaded mapfoldl(identity, []));
end

function testThreadedMapfoldr()::Bool
	return compareGeneratedCode("tmapfoldr(identity, [])", @macroexpand @threaded mapfoldr(identity, []));
end

function testThreadedMap()::Bool
	return compareGeneratedCode("tmap(identity, [])", @macroexpand @threaded map(identity, []));
end

function testThreadedMap!()::Bool
	return compareGeneratedCode("tmap!(identity, [])", @macroexpand @threaded map!(identity, []));
end

function testThreadedMapreduce()::Bool
	return compareGeneratedCode("tmaptreduce(identity, [])", @macroexpand @threaded mapreduce(identity, []));
end

function testThreadedReduce()::Bool
	return compareGeneratedCode("treduce(identity, [])", @macroexpand @threaded reduce(identity, []));
end

# Test that the macro doesn't mess up other functions
function testThreadedOther()::Bool
	return compareGeneratedCode("eltype(identity, [])", @macroexpand @threaded eltype(identity, []));
end

function testThreadedNested()::Bool
	return compareGeneratedCode("tmap(tmap!(identity), [[]])", @macroexpand @threaded map(map!(identity), [[]]));
end


text:
module AGCD
using Polynomials
using LinearAlgebra
#using SparseArrays # didn't get boost I thought we might with JF

## This provide AGCD.agcd for finding an *approximate* GCD of two polynomials. The most common use in this package
## is to reduce a polynomial `p` to a square free polynomial `q=p/gcd(p, p')`.

## This isn't too accurate for higher order polys, but
## seems to be okay now for degree 20 or less
## Can we do better?
## Follow up on these:
# * [Winkler&Hasan](http://ac.els-cdn.com/S0377042712003135/1-s2.0-S0377042712003135-main.pdf?_tid=126266d0-ea31-11e6-a1a7-00000aab0f6b&acdnat=1486140896_2119601f37cf58ff0d3eefbc05fecbb4)
# * [LiLiuZhi](http://www.sciencedirect.com/science/article/pii/S0377042707000271?np=y&npKey=f5e3b7a1bc9583214348aad066170ea389212ebba443454938cacfdbf96177f8)
# * [Winkler](http://ml.dcs.shef.ac.uk/summer_school/ln.pdf)
# * [Matonardi&...](http://www.sciencedirect.com/science/article/pii/S0024379515005649)
# * [FarmerPhD](http://www.learninglink.bbk.ac.uk/site/assets/files/1025/farmer.pdf)
# * [WinklerLaoHasan](http://www.sciencedirect.com/science/article/pii/S0377042712000751)
# * [Sagraloff&...](https://arxiv.org/pdf/1308.4088v2.pdf)
# * [](https://arxiv.org/pdf/1207.0630v3.pdf)

## ------------------------

# we use vectors, not polynomials. Some of this if copied and modified from Polynomials.jl
# originally by @vtjnash
_degree(p::Vector) = length(p) - 1
function _polyval(p::Vector{T}, x::S) where {T, S}
    R = promote_type(T,S)
    y = convert(R, p[end])
    for i in lastindex(p)-1:-1:1
        y = p[i] + x * y
    end
    y
end
function _monic!(ps::Vector{T}) where T
    lambda = iszero(ps[end]) ? one(T) : 1 / ps[end]
    ps .= ps .* lambda
end
function _monic(p::Vector)
    p1 = copy(p)
    _monic!(p1)
    p1
end

function _polymul(p::Vector{T}, q) where {T}
    R = promote_type(T,eltype(q))
    n = length(p)-1
    m = length(q)-1
    a = zeros(R, m+n+1)

    for i = 1:n+1
        for j = 1:m+1
            a[i+j-1] += p[i] * q[j]
        end
    end
    a
end

function _polyder(p::Vector{T}) where {T}
    S = eltype(float(p[1]))
    n = length(p)
    a2 = Vector{S}(undef, n-1)
    for i = 1:n-1
        a2[i] = p[i+1] * i
    end
    a2
end

## ------------------------

## Special matrices
## p = [p_0, p_1, ..., p_n]
function cauchy_matrix_size(p::Vector, k::Integer)
    n = length(p)
    (n + k - 1, k)
end
function cauchy_matrix!(M, p::Vector{T}, k::Integer) where {T}
    n = length(p)
    for i in 1:k
        M[(1:n) .+ (i-1), i] = p
    end
end
function cauchy_matrix(p::Vector{T}, k::Integer) where {T}
    n,m = cauchy_matrix_size(p, k)
    out = zeros(T, n, m)
    cauchy_matrix!(out, p, k)
    out
end

## Jacobian F(u,v,w) = [p,p'] is J(u,v,w)
function JF_size(u, v, w)

    m, k, j = _degree(u), _degree(v), _degree(w)
    n, l = m + k, m + j

    ai, aj = cauchy_matrix_size(v, m + 1)
    bi, bj = cauchy_matrix_size(u, k + 1)
    di, dj = cauchy_matrix_size(w, m + 1)
    fi, fj = cauchy_matrix_size(u, j + 1)
    ci, cj = ai, fj
    ei, ej = di, bj

    (1 + ai + di, aj + bj + cj)

end

function JF!(M, u::Vector{T}, v, w) where {T}
    m, k, j = _degree(u), _degree(v), _degree(w)
    n, l = m + k, m + j

    # JF_size should return these
    ai, aj = cauchy_matrix_size(v, m + 1)
    bi, bj = cauchy_matrix_size(u, k + 1)
    di, dj = cauchy_matrix_size(w, m + 1)
    fi, fj = cauchy_matrix_size(u, j + 1)
    ci, cj = ai, fj
    ei, ej = di, bj

    M[1,1] = one(T)
    cauchy_matrix!(view(M, 1 .+ (1:ai), 1:aj),        v, m + 1)
    cauchy_matrix!(view(M, 1 .+ (1:bi), aj .+ (1:bj)), u, k + 1)
#    M[1 + (1:ci), aj + bj + (1:cj)] = zeros(T, n+1, m+1-j)
    cauchy_matrix!(view(M, 1 + ai .+ (1:di), 1:dj),   w, m + 1)
#    M[1 + bi + (1:ei), dj + (1:ej)] = zeros(T, m + 1, n + 1 - j)
    cauchy_matrix!(view(M, 1 + ci .+ (1:fi), dj + ej .+ (1:fj)), u, j + 1)

end
function JF(u::Vector{U}, v::Vector{V}, w::Vector{W}) where {U,V,W}
    R = promote_type(U,V, W)
    M = zeros(R, JF_size(u, v, w)...)
    JF!(M, u, v, w)
    M
end

# eqn (12) to get weights
function _weights(p::Vector{T},q) where {T}
    n, m = length(p), length(q)
    wts = ones(T, 1 + n + m)
    for (i,pj) in enumerate(p)
        wts[1+i] = min(one(real(T)), inv(abs(pj)))
    end
    for (i,qj) in enumerate(q)
        wts[1 + n + i] = min(one(real(T)), inv(abs(qj)))
    end
    wts
end

# reduce residual by Gauss-Newton algorithm
# updates u,v,w; returns residual error and flag for convergence
function reduce_residual!(u,v,w, p::Vector{T}, q, wts, ρ) where {T}
    m, n, l = map(_degree, (u, v, w))
    # preallocate
    #R = big(T)
    R = T
    A = zeros(R, JF_size(u, v, w)...)
    b = zeros(R, 1 + length(p) + length(q))
    inc = zeros(R, m + n + l + 3) #weighted_least_square(A, b, wts)
    up,vp,wp = copy(u),copy(v),copy(w)

    ρm, ρm_1 = one(T), residual_error(p,q,u,v,w, wts)
    MAXSTEPS = 10
    ctr = 0
    flag = :not_converged

    rho_0 = ρm_1

    while ctr <= MAXSTEPS

        ρm =  agcd_update!(p, q, u,v,w, m, n, A, b, inc, up, vp, wp, wts, ρm_1)

        if ρm <  ρm_1
            ctr += 1
            ρm_1 = ρm
        else
            break
        end
    end

    # threshold here is really important, but not clear what it should be.
    # We empirically get better results---it seems---with sqrt(norm(p,2))
    # though recomendation if norm(u,2); Might also make sense to include
    # q in the computation
    if ρm <= ρ * norm(u,2)
        flag = :converged
    end
    if  ρm < rho_0
    end
    return (ρm, flag)
end

## Lemma 4.1, (25)
## solve weighted least squares problem W*(Ax - b) = 0
function weighted_least_square(A, b, w)
    W = diagm(0 => w)
    (W * A) \ (W * b)
end

function weighted_least_square!(M, A, b, w)
    W = diagm(0 => w)
    M[:] .= (W * A) \ (W * b)
end

## compute F(u,v,w) - [p, p'] = [u*v, u*w] - [p, p']
Fmp(p,q,u,v,w) =  [u[end]; _polymul(u,v); _polymul(u,w)] .- [1; p; q]
function Fmp!(b, p,q,u,v,w)
    b[1] = u[end] - 1
    offset = 1
    for (i,v) in enumerate(_polymul(u,v) .- p)
        b[offset + i] = v
    end
    offset = 1 + length(p)
    for (i,v1) in enumerate(_polymul(u,w) .- q)
        b[offset + i] = v1
    end
end

residual_error(p,q,u,v,w, wts=ones(1 + length(p) + length(q))) = norm( ([_polymul(u,v); _polymul(u,w)] .- [p; q]) .* wts[2:end], 2)

### refine u,v,w estimate using
## Newton's method and weighted least squares
# updates in place u,v,w
function agcd_update!(p, q, u,v,w,m,n, A, b, inc, up, vp, wp, wts, err0)

    JF!(A, u,v,w)
    Fmp!(b, p,q,u,v,w)

    weighted_least_square!(inc, A, b, wts)

    Δu = inc[1:(1+m)]
    Δv = inc[(m+2):(m+n+2)]
    Δw = inc[(m+n+3):end]

    up .= u; vp .= v; wp .= w

    up .-=  Δu; _monic!(up)
    vp .-=  Δv; _monic!(vp)
    wp .-=  Δw; _monic!(wp)

    err = residual_error(p,q,up,vp,wp, wts)
    # update if we improve
    if err  < err0
        u .= up; v .= vp; w .= wp

    else
        err = err0
    end

    # return error estimate

    err
end

# used below
function _mult(Gs, A)
    for G in Gs
        A = G * A
    end
    A
end

function qr_sylvester!(Gs, A0, k)

        A = copy(A0)
        m,n = size(A)
        for i in m:-1:2
            Gi, r = givens(A[1,1],A[i,1], 1, i)
            A = Gi * A
            push!(Gs, Gi)  # [G1, G2, ..., Gn]
        end

        for i in m:-1:3
            Gi, r = givens(A[2,2],A[i,2], 2, i)
            A = Gi * A
            push!(Gs, Gi)
        end
        return A[1:2,1:2]
end

function qr_sylvester!(Gs, A0, k, R)
    T = eltype(A0)
    N = (k-1)
    Zs = zeros(T, N)
    B = vcat(hcat(Zs, Zs), A0)
    B = _mult(Gs, B)
    m,n = size(B)
    for i in m:-1:(2k)
        Gi, r = givens(B[2(k-1)+1, 1],   B[i, 1], 2(k-1)+1, i)
        B = Gi * B
        push!(Gs, Gi)
    end
    for i in m:-1:(2k+1)
            Gi, r = givens(B[2k, 2], B[i, 2], 2k, i)
        B = Gi * B
        push!(Gs, Gi)
    end
    return [vcat(R, zeros(T, 2, 2*(k-1))) B[1:(2k),:]]
end

# return sigma, x
# converge on smallest eigenvalue of (A'*A) using power rule
# A=QR; (A'*A)^-1) = (R'*Q'*Q*R)^(-1) = (R'*R)^(-1)
# instead of computing x_{i+1} = (R'*R)^{-1} xi; we solve R'y=xi; Rz=y;x=z/|z|
function smallest_eigval(R::LinearAlgebra.UpperTriangular{T},
                         thresh=sqrt(eps(real(T)))) where {T}
    k = size(R)[1]/2

    if iszero(det(R))
        return (:iszero, zero(T), T[])
    end

    m,n = size(R)

    x = ones(T, n)
    y = zeros(T, m)
    z = zeros(T, n)
    σ, σ1 = Inf*one(real(T)), Inf*one(real(T))
    flag = :ispositive

    for cnt in 1:10
        y .= R' \ x
        z .= R  \ y
        nz = 1/norm(z,2)
        x .= z .* nz
        sigma = norm(R * x, 2)
        σ1 = abs(sigma)

        if σ1 < 0.99 *  σ
            σ = σ1
            continue
        end

        if σ1 < thresh
            flag = :ispossible
            break
        end
        if  abs(σ - σ1)  < 1.1 * σ1
            flag = :ispossible
            break
        end

        σ = σ1
    end

    return (flag, σ1, x)
end

## return k, sigma for possible k values
function sylvester_matrix_singular_values(p::Poly, q::Poly=polyder(p);kwargs...)
    sylvester_matrix_singular_values(Polynomials.coeffs(p), Polynomials.coeffs(q); kwargs...)
end

function sylvester_matrix_singular_values(ps::Vector{T},
                                          qs::Vector{S}=_polyder(ps),
                                          θ=sqrt(eps(real(T)))) where {T,S}
    _monic!(ps); _monic!(qs)
    n, m = length(ps), length(qs)
    wts = _weights(ps,qs)

    R = promote_type(T,S)
    A0::Matrix{R} = R[ps vcat(zeros(T, n-m),qs)]
    Gs = LinearAlgebra.Givens{T}[]
    k = 1
    sigmas = R[]
    R = qr_sylvester!(Gs, A0, k)

    local x::Vector{T}

    while k <= m-1
        V = UpperTriangular(R)
        flag, sigma, x = smallest_eigval(V, θ * norm(ps,2))
        push!(sigmas, sigma)
        k += 1
        R =  qr_sylvester!(Gs, A0, k, R)
    end
    return [1:(k-1) sigmas]
end

# return R for sylvester matrix of size k
function sylvester_matrix_k(ps::Vector{T},
                            K) where {T}
    qs = _polyder(ps)
    _monic!(ps); _monic!(qs)
    n, m = length(ps), length(qs)
    wts = _weights(ps,qs)

    A0::Matrix{T} = [ps vcat(zeros(T, n-m),qs)]
    Gs = LinearAlgebra.Givens{T}[]
    R = qr_sylvester!(Gs, A0, 1)
    k = 2
    while k <= K
        R =  qr_sylvester!(Gs, A0, k, R)
        k += 1
    end
    V = UpperTriangular(R)
end

"""
    agcd(p, q, θ=sqrt(eps(real(T))), ρ=cbrt(eps(real(T))) * θ)

Find an approximate GCD for polynomials `p` and `q` using an algorithm of [Zeng](https://doi.org/10.1090/S0025-5718-04-01692-8).

Returns u,v,w, err where:

* `u*v ≈ monic(p)`
* `u*w ≈ monic(q)`
* The total residual error in these approximations is bounded by `err`.

Further,

* `v` and `w` should share no common roots (`u` is a gcd of `u*v` and `u*w`)
* the roots of `v` should exhaust the unique values of roots of `p`.

(If `p` and `q` are specified as vectors, the returned values will be
vectors. If specified as objects of `Polynomial` type then the
returned values will be as well.)

The tolerances are:

* θ: singular value threshold. Used to identify if smallest singular value of Sylvester matrix is ≈ 0. We use `sqrt(eps(T))`, close to that of Zeng's `1e-8`, but more flexible should `BigFloat` values be used.

* ρ: initial residual tolerance. If we can get (u,v,w) error less than this, we stop. Zeng uses `1e-10`, we use as a default the smaller `eps(T)^(5/6)`.

The algorithm looks for the first `k` for which the corresponding
Sylvester matrix is rank deficient. This follows Lemma 2.4 of the paper, which
finds the smallest singular value. In the process, an estimate for
`u`,`v`, and `w` is produced. If this singular value is approximately
0 (as determined with the parameter θ) *and* the estimates of the
decomposition can be refined via a Gauss-Jordan process to a close
enough approximation (as determined with the parameter ρ) then the
desired values are found.

The act of identifying the value of `k` depends on determining if the
smallest singular value is 0 which in practice is quite hard. The
simple threshold of `norm(p,2) ⋅ θ` will fail for higher-degree
polynomials. For manually identifying this value, the complete set of
estimated singular values for different values of `k` are returned by
`sylvester_matrix_singular_values`. For a given `k` the refined values
of `u`, `v`, and `w` are returned by `rank_k_agcd`.

"""
function agcd(ps::Vector{T}, qs::Vector{S}=_polyder(ps);
              θ=sqrt(eps(real(T))),
              ρ=cbrt(eps(real(T))) * θ,
              maxk=length(qs)) where {T,S}
    _monic!(ps); _monic!(qs)
    n, m = length(ps), length(qs)
    wts = _weights(ps,qs)

    U = promote_type(T,S)
    A0::Matrix{U} = U[ps vcat(zeros(U, n-m),qs)]

    nm = norm(ps, 2)

    thresh = nm * θ ## this is sensitive

    Gs = LinearAlgebra.Givens{T}[]
    k = 1
    R = qr_sylvester!(Gs, A0, k)

    local x::Vector{U}
    while k <=  maxk
        V = UpperTriangular(R)
        flag, sigma, x = smallest_eigval(V, thresh)

        if flag == :iszero
            # exact zero means we need to identify initial guess for u,v,w
            u = qs
            w = ones(U, 1)
            v = cauchy_matrix(u, length(ps) - m + 1) \ ps
            _monic!(v)

            return (u,v,w, zero(U))

        end

        if flag != :iszero &&  sigma < thresh
            v = x[2:2:end]; _monic!(v)
            w = x[(1+2(n-m)):2:end]; _monic!(w)
            A =  cauchy_matrix(v, n - length(v) + 1)
            u = A \ ps
            _monic!(u)

            ρm, flag = reduce_residual!(u, v, w, ps, qs, wts, ρ)

            if flag == :converged || flag == :not_updated
                return (u, v, w, ρm)
            end

        end

        # else bumpup k
        k += 1
        if k > m
            u = ones(U,1)
            w = qs
            v = ps
            return (u, v, w, zero(U))
        end
        R =  qr_sylvester!(Gs, A0, k, R)
    end

    # Okay, we gave it a good shot, but didn't find a perfect candidate
    # let's move on...
    v = x[2:2:end]; _monic!(v)
    w = x[3:2:end]; _monic!(w)
    A =  cauchy_matrix(v, n - length(v) + 1)
    u = A \ ps
    _monic!(u)

    wts = _weights(ps,qs)
    ρm, flag = reduce_residual!(u,v,w, ps, qs, wts, ρ)
    return (u, v, w, ρm)

end

function agcd(p::Poly, q::Poly=polyder(p); kwargs...)
    u,v,w,rho = agcd(Polynomials.coeffs(p), Polynomials.coeffs(q); kwargs...)
    (Poly(u), Poly(v), Poly(w), rho)
end

## --------------
## preconditioning code
## taken from https://who.rocq.inria.fr/Jan.Elias/pdf/je_masterthesis.pdf

## compute q(x) = p(phi*x)
function Tphi!(p, ps::Vector{T}, phi) where {T}
    p .= ps .* (phi^(i-1) for i in eachindex(ps))
end
function Tphi(ps::Vector{T}, phi) where {T}
    p = copy(ps)
    p .= ps .* (phi^(i-1) for i in eachindex(ps))
    p
end
## c
## compute q(x) = p(x-alpha)
function T_alpha(ps::Vector{T}, alpha) where {T}
end

function geometric_mean(a::Vector{T})::T where {T}
    b = filter(!iszero, a)
    n = length(b)
    prod(abs(u)^(one(T)/n) for u in b)
end

# compute max(max(ps), max(qs)) / min(min(ps), min(qs))
function ratio(p::Vector{T}, q::Vector{T}) where {T}
    m,M = Inf*one(real(T)), zero(real(T))
    for ps in (p, q)
        for a in ps
            aa = abs(a)
            if 0 < aa < m
                m = aa
            end
            if aa > M
                M = aa
            end
        end
    end

    iszero(m) ? Inf*one(T) : M/m

end

## find alpha, gamma that minimize ratio of max coefficients to min coefficients, for getting zeros
## 1.12 writes this as a linear programming problem, we just ...
function precondition(p::Vector{T}, q::Vector) where {T}

    S = real(T)
    m = ratio(p,q)

    alphas = ((2*one(S))^i for i in -5:5)
    phis = ((2*one(S))^i for i in -5:5)

    out = (one(S), one(S))

    for α in alphas
        for ϕ in phis
            r = ratio(Tphi(p, ϕ), α * Tphi(q, ϕ))
            if r < m
                out = (α, ϕ)
            end
        end
    end

    α, ϕ = out

    p = Tphi(p, ϕ)
    q = α * Tphi(q, ϕ)

    p = p * (1/geometric_mean(p))
    q = q * (1/geometric_mean(q))

    p, q, ϕ, α

end

end


text:
module day8

using UnicodePlots

readInput() = [parse(Int, i) for i in readline(joinpath(@__DIR__, "input.txt"))]

function part1(input, w=25, h=6)
    layers = reshape(input, w, h, :)
    layermin = last(findmin(mapslices(x -> count(==(0), x), layers, dims=(1,2))))[3]
    count(==(1), layers[:,:,layermin]) * count(==(2), layers[:,:,layermin])
end

function part2(input, w=25, h=6)
    layers = reshape(input, w, h, :)
    img = fill(-1, w, h)
    for i = 1:w, j = 1:h
        k = 1
        while layers[i,j,k] == 2
            k += 1
        end
        img[i, j] = layers[i,j,k]
    end
    coords = findall(==(1), img)
    scatterplot([c[1] for c in coords], [h-c[2] for c in coords] ; canvas = BlockCanvas, height = h)
end

end


text:
export ScoreBrier

"""
    ScoreBrier <: ScoringRule

 Calculate Brier Scores for a set of outcomes and predictions

## References
* https://en.wikipedia.org/wiki/Scoring_rule
* https://en.wikipedia.org/wiki/Brier_score
* https://www2.cs.duke.edu/courses/spring17/compsci590.2/proper_scoring.pdf

## Parameters
* `normalise::false`: if true, normalise scores so that for binary choice
     + S( [0.5,0.5], 1) = 0.0
     + S( [0.0,1.0], 1) = 1.0

```
"""
struct ScoreBrier <: ScoringRule
    normalise::Bool
end
ScoreBrier(; normalise::Bool=false) = ScoreBrier(normalise)

function score_direction(srule::ScoreBrier)
    if srule.normalise
        return 1 # normalising puts it in the positive direction
    else
        return -1
    end
end

function scoring_function(srule::ScoreBrier, predicted_probabilities::Array{Float64,1}, outcome::Int)
    # general model allows C classes, so number classes 1,2,...,C
    #    predicted_probabilities is a C-element vector, which must sum to 1
    #    outcome is the index of the event that actually happens
    if abs( sum(predicted_probabilities) - 1.0 ) > 1.0e-2
        error("probabilities must sum to 1")
    end
    C = length(predicted_probabilities)
    e_i = zeros(Float64, C)
    e_i[outcome] = 1.0
    result = sum( (predicted_probabilities - e_i).^2 )
    if srule.normalise
        return 1-2*result
    else
        return result
    end
end


text:
export Shamanskii

function Shamanskii(nlp_stop,
					PData :: PDataFact,
					dₙ :: Vector,
					gt :: Vector;
					λfact2 :: Bool = false)
	# printstyled("on est dans Shamanskii \n", color = :bold)
	nlp_at_x = nlp_stop.current_state
	x = nlp_at_x.x
	n = length(x)
	xdemi = x + dₙ
	gtemp = grad(nlp_stop.pb, xdemi)

	ϵ2 =  sqrt(eps(eltype(gt)))
   	Γ = max.(abs.(PData.Δ), ϵ2)

	if (PData.λ == 0) && !(λfact2)
		d̃ = PData.L \ gtemp[PData.pp]
    	d̂ = PData.L' \ (PData.Q * (PData.Q' * d̃ ./ Γ))
		d = -d̂[invperm(PData.pp)]
	else
		# println("on est ici 🐵↓")
		# @show (nlp_at_x.Hx + PData.λ * Matrix(1.0I, n, n))
		L, D, pp, ρ, ncomp = ldlt_symm((nlp_at_x.Hx + PData.λ * Matrix(1.0I, n, n)))
		# println("On a L")
		D = Hermitian(D)
	    X = eigen(D)
	    Δ = X.values
	    Q =  X.vectors
		d̃ = L \ gtemp[pp]
    	d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
		d = -d̂[invperm(pp)]
	end

	dₕₒ = dₙ .+ d

	return dₕₒ
end


text:
# transforms.jl
# =============
# Transformations on DAG:s
#

load("utils/req.jl")
req("utils/cached.jl")
req("dag/dag.jl")

# -- evaluation ---------------------------------------------------------------

type Context{C}
    c::C
    cache::Cache

    Context(c::C) = new(c, Cache())
    Context() = new(C(), Cache())
end
Context{C}(c::C) = Context{C}(c)

cacheentry(f::Function, c::Context, args...) = cacheentry(f, c.cache, args...)

evaluate(c::Context, ns::Nodes) = { (@cached evaluate(c, node))|node in ns }

default_evaluate(c::Context, node::Node) = Node(node, evaluate(c, node.args))
evaluate(c::Context, node::Node) = default_evaluate(c, node)

# -- topsort iteration --------------------------------------------------------

type Topsorter; end
typealias TopsortContext Context{Topsorter}

forward(sink::Node) = @task evaluate(TopsortContext(), sink)

function evaluate(c::TopsortContext, node::Node)
    evaluate(c, node.args)
    produce(node)
    nothing
end

# -- nodewise rewrite ---------------------------------------------------------

type Rewriter; rewrite::Function; end
typealias RewriteContext Context{Rewriter}

rewrite_dag(sink::Node, fun::Function) = evaluate(Context(Rewriter(fun)), sink)
function evaluate(c::RewriteContext, node::Node)
    args = evaluate(c, node.args)
    c.c.rewrite(node, args)
end

# -- Name nodes with fanout > 1 -----------------------------------------------

name_fanout_nodes(sink::Node) = rewrite_dag(sink, namefanout_rewrite)

function namefanout_rewrite(oldnode::Node, args::Vector)
    for arg in args
        nu = (arg.num_uses += 1)
        if !has_name(arg) && (nu >= 2) && isa(arg, OpNode)
            arg.name = gensym()
        end
    end
    Node(oldnode, args)
end

# -- Make sure named nodes have unique names ----------------------------------

function name_nodes_uniquely(sink::Node)
    names = Set{Symbol}()
    rewrite_dag(sink, (node, args)->nameuniquely_rewrite(names, node, args))
end

function nameuniquely_rewrite(names::Set{Symbol}, oldnode::Node, args::Vector)
    node = Node(oldnode, args)
    if has_name(node)
        if has(names, get_name(node))
            node.name = gensym()
        end
        add(names, node.name)
    end
    node
end

# -- Collect SymNode names ----------------------------------------------------

function collect_symnode_names(sink::Node)
    allnames = Dict{Symbol,Vector{Symbol}}()
    for node in forward(sink)
        if isa(node, SymNode)
            names = @setdefault allnames[node.val.kind] Symbol[]
            push(names, node.val.name)
        end
    end
    allnames
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@doc raw"""
    scatter_ray(
        ::Type{DiffuseBRDF},
        pcg::PCG,
        incoming_dir::Vec,
        interaction_point::Point,
        normal::Normal,
        depth::Int64,
        ) :: Ray

Return a `Ray` scattered by a material with a `DiffuseBRDF`.

A `DiffuseBRDF` has a uniform BRDF, i.e.
``f_r(\mathbf{x}, \mathbf{\Psi}\rightarrow\mathbf{\Theta}) = \rho_d / \pi``;
the importance sampling for the `PathTracer` algorithm use consequently the
following PDF:

```math
p(\omega) \propto
    f_r(\mathbf{x}, \mathbf{\Psi}\rightarrow\mathbf{\Theta}) \, \cos(\vartheta)
    = \frac{\rho_d}{\pi} \, \cos(\vartheta)
    \propto \cos(\vartheta)
```
```math
    \Rightarrow \quad
p(\omega) = \frac{\cos(\vartheta)}{\pi}
    \quad \Rightarrow \quad
p(\vartheta ,\varphi) = \frac{\cos(\vartheta) \, \sin(\vartheta) }{2\pi}
```
```math
\Rightarrow \quad
\begin{aligned}
    &p(\vartheta) = 2 \, \cos(\vartheta) \, \sin(\vartheta) \\
    &p(\varphi | \vartheta) = \frac{1}{2 \pi}
\end{aligned}
```

See also: [`DiffuseBRDF`](@ref), [`Ray`](@ref), [`Vec`](@ref),
[`Point`](@ref), [`Normal`](@ref), [`PCG`](@ref)
"""
function scatter_ray(
            ::Type{DiffuseBRDF},
            pcg::PCG,
            incoming_dir::Vec,
            interaction_point::Point,
            normal::Normal,
            depth::Int64,
            )

    e1, e2, e3 = create_onb_from_z(normal)
    cos_θ_sq = random(pcg)
    cos_θ, sin_θ = √(cos_θ_sq) , √(1.0 - cos_θ_sq)
    ϕ = 2.0 * pi * random(pcg)

    return Ray(
                interaction_point,
                e1*cos(ϕ)*cos_θ + e2*sin(ϕ)*cos_θ + e3*sin_θ,
                1.0e-3,   # tmin, be generous here
                Inf,
                depth
            )
end

@doc raw"""
    scatter_ray(
        ::Type{SpecularBRDF},
        pcg::PCG,
        incoming_dir::Vec,
        interaction_point::Point,
        normal::Normal,
        depth::Int64,
        ) :: Ray

Return a `Ray` scattered by a material with a `SpecularBRDF`.

A `SpecularBRDF` has a Dirac delta BRDF, i.e.:

```math
f_r(\mathbf{x}, \mathbf{\Psi} \rightarrow \mathbf{\Theta})
    \propto
\frac{\delta(\sin^2\theta_r - \sin^2\theta) \,
    \delta(\psi_r \pm \pi - \psi)}{\cos\theta},
```

The importance sampling for the `PathTracer` algorithm use consequently the
following PDF:

```math
p(\omega) \propto
    f_r(\mathbf{x}, \mathbf{\Psi}\rightarrow\mathbf{\Theta}) \, \cos(\vartheta)
    \propto \frac{1}{\cos(\vartheta)} \, \cos(\vartheta)
    \propto cost
```

See also: [`SpecularBRDF`](@ref), [`Ray`](@ref), [`Vec`](@ref),
[`Point`](@ref), [`Normal`](@ref), [`PCG`](@ref)
"""
function scatter_ray(
            ::Type{SpecularBRDF},
            pcg::PCG,
            incoming_dir::Vec,
            interaction_point::Point,
            normal::Normal,
            depth::Int64,
            )
    #ray_dir = normalize(Vec(incoming_dir.x, incoming_dir.y, incoming_dir.z))
    ray_dir = normalize(incoming_dir)

    return Ray(
                interaction_point,
                ray_dir - normal*2.0*(Vec(normal) ⋅ ray_dir),
                1e-3,
                Inf,
                depth
            )
end

#=
scatter_ray(
            b::BRDF,
            pcg::PCG,
            incoming_dir::Vec,
            interaction_point::Point,
            normal::Normal,
            depth::Int64,
        ) =
        scatter_ray(
            typeof(b),
            pcg,
            incoming_dir,
            interaction_point,
            normal,
            depth,
        )
=#


text:
const text_bb_map = Dict(
    'A' => '𝔸',
    'B' => '𝔹',
    'C' => 'ℂ',
    'D' => '𝔻',
    'E' => '𝔼',
    'F' => '𝔽',
    'G' => '𝔾',
    'H' => 'ℍ',
    'I' => '𝕀',
    'J' => '𝕁',
    'K' => '𝕂',
    'L' => '𝕃',
    'M' => '𝕄',
    'N' => 'ℕ',
    'O' => '𝕆',
    'P' => 'ℙ',
    'Q' => 'ℚ',
    'R' => 'ℝ',
    'S' => '𝕊',
    'T' => '𝕋',
    'U' => '𝕌',
    'V' => '𝕍',
    'W' => '𝕎',
    'X' => '𝕏',
    'Y' => '𝕐',
    'Z' => 'ℤ',
    'a' => '𝕒',
    'b' => '𝕓',
    'c' => '𝕔',
    'd' => '𝕕',
    'e' => '𝕖',
    'f' => '𝕗',
    'g' => '𝕘',
    'h' => '𝕙',
    'i' => '𝕚',
    'j' => '𝕛',
    'k' => '𝕜',
    'l' => '𝕝',
    'm' => '𝕞',
    'n' => '𝕟',
    'o' => '𝕠',
    'p' => '𝕡',
    'q' => '𝕢',
    'r' => '𝕣',
    's' => '𝕤',
    't' => '𝕥',
    'u' => '𝕦',
    'v' => '𝕧',
    'w' => '𝕨',
    'x' => '𝕩',
    'y' => '𝕪',
    'z' => '𝕫',
    '0' => '𝟘',
    '1' => '𝟙',
    '2' => '𝟚',
    '3' => '𝟛',
    '4' => '𝟜',
    '5' => '𝟝',
    '6' => '𝟞',
    '7' => '𝟟',
    '8' => '𝟠',
    '9' => '𝟡'
)

const text_it_map = Dict(
    'A' => '𝐴',
    'B' => '𝐵',
    'C' => '𝐶',
    'D' => '𝐷',
    'E' => '𝐸',
    'F' => '𝐹',
    'G' => '𝐺',
    'H' => '𝐻',
    'I' => '𝐼',
    'J' => '𝐽',
    'K' => '𝐾',
    'L' => '𝐿',
    'M' => '𝑀',
    'N' => '𝑁',
    'O' => '𝑂',
    'P' => '𝑃',
    'Q' => '𝑄',
    'R' => '𝑅',
    'S' => '𝑆',
    'T' => '𝑇',
    'U' => '𝑈',
    'V' => '𝑉',
    'W' => '𝑊',
    'X' => '𝑋',
    'Y' => '𝑌',
    'Z' => '𝑍',
    'a' => '𝑎',
    'b' => '𝑏',
    'c' => '𝑐',
    'd' => '𝑑',
    'e' => '𝑒',
    'f' => '𝑓',
    'g' => '𝑔',
    'h' => 'ℎ',
    'i' => '𝑖',
    'j' => '𝑗',
    'k' => '𝑘',
    'l' => '𝑙',
    'm' => '𝑚',
    'n' => '𝑛',
    'o' => '𝑜',
    'p' => '𝑝',
    'q' => '𝑞',
    'r' => '𝑟',
    's' => '𝑠',
    't' => '𝑡',
    'u' => '𝑢',
    'v' => '𝑣',
    'w' => '𝑤',
    'x' => '𝑥',
    'y' => '𝑦',
    'z' => '𝑧',
    'Α' => '𝛢',
    'Β' => '𝛣',
    'Γ' => '𝛤',
    'Δ' => '𝛥',
    'Ε' => '𝛦',
    'Ζ' => '𝛧',
    'Η' => '𝛨',
    'Θ' => '𝛩',
    'Ι' => '𝛪',
    'Κ' => '𝛫',
    'Λ' => '𝛬',
    'Μ' => '𝛭',
    'Ν' => '𝛮',
    'Ξ' => '𝛯',
    'Ο' => '𝛰',
    'Π' => '𝛱',
    'Ρ' => '𝛲',
    'ϴ' => '𝛳',
    'Σ' => '𝛴',
    'Τ' => '𝛵',
    'Υ' => '𝛶',
    'Φ' => '𝛷',
    'Χ' => '𝛸',
    'Ψ' => '𝛹',
    'Ω' => '𝛺',
    '∇' => '𝛻',
    'α' => '𝛼',
    'β' => '𝛽',
    'γ' => '𝛾',
    'δ' => '𝛿',
    'ε' => '𝜀',
    'ζ' => '𝜁',
    'η' => '𝜂',
    'θ' => '𝜃',
    'ι' => '𝜄',
    'κ' => '𝜅',
    'λ' => '𝜆',
    'μ' => '𝜇',
    'ν' => '𝜈',
    'ξ' => '𝜉',
    'ο' => '𝜊',
    'π' => '𝜋',
    'ρ' => '𝜌',
    'ς' => '𝜍',
    'σ' => '𝜎',
    'τ' => '𝜏',
    'υ' => '𝜐',
    'φ' => '𝜑',
    'χ' => '𝜒',
    'ψ' => '𝜓',
    'ω' => '𝜔',
    '∂' => '𝜕',
    'ϵ' => '𝜖',
    'ϑ' => '𝜗',
    'ϰ' => '𝜘',
    'ϕ' => '𝜙',
    'ϱ' => '𝜚',
    'ϖ' => '𝜛'
)

const text_bf_map = Dict(
    'A' => '𝐀',
    'B' => '𝐁',
    'C' => '𝐂',
    'D' => '𝐃',
    'E' => '𝐄',
    'F' => '𝐅',
    'G' => '𝐆',
    'H' => '𝐇',
    'I' => '𝐈',
    'J' => '𝐉',
    'K' => '𝐊',
    'L' => '𝐋',
    'M' => '𝐌',
    'N' => '𝐍',
    'O' => '𝐎',
    'P' => '𝐏',
    'Q' => '𝐐',
    'R' => '𝐑',
    'S' => '𝐒',
    'T' => '𝐓',
    'U' => '𝐔',
    'V' => '𝐕',
    'W' => '𝐖',
    'X' => '𝐗',
    'Y' => '𝐘',
    'Z' => '𝐙',
    'a' => '𝐚',
    'b' => '𝐛',
    'c' => '𝐜',
    'd' => '𝐝',
    'e' => '𝐞',
    'f' => '𝐟',
    'g' => '𝐠',
    'h' => '𝐡',
    'i' => '𝐢',
    'j' => '𝐣',
    'k' => '𝐤',
    'l' => '𝐥',
    'm' => '𝐦',
    'n' => '𝐧',
    'o' => '𝐨',
    'p' => '𝐩',
    'q' => '𝐪',
    'r' => '𝐫',
    's' => '𝐬',
    't' => '𝐭',
    'u' => '𝐮',
    'v' => '𝐯',
    'w' => '𝐰',
    'x' => '𝐱',
    'y' => '𝐲',
    'z' => '𝐳',
    'Α' => '𝚨',
    'Β' => '𝚩',
    'Γ' => '𝚪',
    'Δ' => '𝚫',
    'Ε' => '𝚬',
    'Ζ' => '𝚭',
    'Η' => '𝚮',
    'Θ' => '𝚯',
    'Ι' => '𝚰',
    'Κ' => '𝚱',
    'Λ' => '𝚲',
    'Μ' => '𝚳',
    'Ν' => '𝚴',
    'Ξ' => '𝚵',
    'Ο' => '𝚶',
    'Π' => '𝚷',
    'Ρ' => '𝚸',
    'ϴ' => '𝚹',
    'Σ' => '𝚺',
    'Τ' => '𝚻',
    'Υ' => '𝚼',
    'Φ' => '𝚽',
    'Χ' => '𝚾',
    'Ψ' => '𝚿',
    'Ω' => '𝛀',
    '∇' => '𝛁',
    'α' => '𝛂',
    'β' => '𝛃',
    'γ' => '𝛄',
    'δ' => '𝛅',
    'ε' => '𝛆',
    'ζ' => '𝛇',
    'η' => '𝛈',
    'θ' => '𝛉',
    'ι' => '𝛊',
    'κ' => '𝛋',
    'λ' => '𝛌',
    'μ' => '𝛍',
    'ν' => '𝛎',
    'ξ' => '𝛏',
    'ο' => '𝛐',
    'π' => '𝛑',
    'ρ' => '𝛒',
    'ς' => '𝛓',
    'σ' => '𝛔',
    'τ' => '𝛕',
    'υ' => '𝛖',
    'φ' => '𝛗',
    'χ' => '𝛘',
    'ψ' => '𝛙',
    'ω' => '𝛚',
    '∂' => '𝛛',
    'ϵ' => '𝛜',
    'ϑ' => '𝛝',
    'ϰ' => '𝛞',
    'ϕ' => '𝛟',
    'ϱ' => '𝛠',
    'ϖ' => '𝛡',
    '0' => '𝟎',
    '1' => '𝟏',
    '2' => '𝟐',
    '3' => '𝟑',
    '4' => '𝟒',
    '5' => '𝟓',
    '6' => '𝟔',
    '7' => '𝟕',
    '8' => '𝟖',
    '9' => '𝟗'
)

const text_frak_map = Dict(
    'A' => '𝔄',
    'B' => '𝔅',
    'C' => 'ℭ',
    'D' => '𝔇',
    'E' => '𝔈',
    'F' => '𝔉',
    'G' => '𝔊',
    'H' => 'ℌ',
    'I' => 'ℑ',
    'J' => '𝔍',
    'K' => '𝔎',
    'L' => '𝔏',
    'M' => '𝔐',
    'N' => '𝔑',
    'O' => '𝔒',
    'P' => '𝔓',
    'Q' => '𝔔',
    'R' => 'ℜ',
    'S' => '𝔖',
    'T' => '𝔗',
    'U' => '𝔘',
    'V' => '𝔙',
    'W' => '𝔚',
    'X' => '𝔛',
    'Y' => '𝔜',
    'Z' => 'ℨ',
    'a' => '𝔞',
    'b' => '𝔟',
    'c' => '𝔠',
    'd' => '𝔡',
    'e' => '𝔢',
    'f' => '𝔣',
    'g' => '𝔤',
    'h' => '𝔥',
    'i' => '𝔦',
    'j' => '𝔧',
    'k' => '𝔨',
    'l' => '𝔩',
    'm' => '𝔪',
    'n' => '𝔫',
    'o' => '𝔬',
    'p' => '𝔭',
    'q' => '𝔮',
    'r' => '𝔯',
    's' => '𝔰',
    't' => '𝔱',
    'u' => '𝔲',
    'v' => '𝔳',
    'w' => '𝔴',
    'x' => '𝔵',
    'y' => '𝔶',
    'z' => '𝔷'
)

const text_cal_map = Dict(
    'A' => '𝓐',
    'B' => '𝓑',
    'C' => '𝓒',
    'D' => '𝓓',
    'E' => '𝓔',
    'F' => '𝓕',
    'G' => '𝓖',
    'H' => '𝓗',
    'I' => '𝓘',
    'J' => '𝓙',
    'K' => '𝓚',
    'L' => '𝓛',
    'M' => '𝓜',
    'N' => '𝓝',
    'O' => '𝓞',
    'P' => '𝓟',
    'Q' => '𝓠',
    'R' => '𝓡',
    'S' => '𝓢',
    'T' => '𝓣',
    'U' => '𝓤',
    'V' => '𝓥',
    'W' => '𝓦',
    'X' => '𝓧',
    'Y' => '𝓨',
    'Z' => '𝓩',
    'a' => '𝓪',
    'b' => '𝓫',
    'c' => '𝓬',
    'd' => '𝓭',
    'e' => '𝓮',
    'f' => '𝓯',
    'g' => '𝓰',
    'h' => '𝓱',
    'i' => '𝓲',
    'j' => '𝓳',
    'k' => '𝓴',
    'l' => '𝓵',
    'm' => '𝓶',
    'n' => '𝓷',
    'o' => '𝓸',
    'p' => '𝓹',
    'q' => '𝓺',
    'r' => '𝓻',
    's' => '𝓼',
    't' => '𝓽',
    'u' => '𝓾',
    'v' => '𝓿',
    'w' => '𝔀',
    'x' => '𝔁',
    'y' => '𝔂',
    'z' => '𝔃'
)

const text_mono_map = Dict(
    'A' => '𝙰',
    'B' => '𝙱',
    'C' => '𝙲',
    'D' => '𝙳',
    'E' => '𝙴',
    'F' => '𝙵',
    'G' => '𝙶',
    'H' => '𝙷',
    'I' => '𝙸',
    'J' => '𝙹',
    'K' => '𝙺',
    'L' => '𝙻',
    'M' => '𝙼',
    'N' => '𝙽',
    'O' => '𝙾',
    'P' => '𝙿',
    'Q' => '𝚀',
    'R' => '𝚁',
    'S' => '𝚂',
    'T' => '𝚃',
    'U' => '𝚄',
    'V' => '𝚅',
    'W' => '𝚆',
    'X' => '𝚇',
    'Y' => '𝚈',
    'Z' => '𝚉',
    'a' => '𝚊',
    'b' => '𝚋',
    'c' => '𝚌',
    'd' => '𝚍',
    'e' => '𝚎',
    'f' => '𝚏',
    'g' => '𝚐',
    'h' => '𝚑',
    'i' => '𝚒',
    'j' => '𝚓',
    'k' => '𝚔',
    'l' => '𝚕',
    'm' => '𝚖',
    'n' => '𝚗',
    'o' => '𝚘',
    'p' => '𝚙',
    'q' => '𝚚',
    'r' => '𝚛',
    's' => '𝚜',
    't' => '𝚝',
    'u' => '𝚞',
    'v' => '𝚟',
    'w' => '𝚠',
    'x' => '𝚡',
    'y' => '𝚢',
    'z' => '𝚣',
    '0' => '𝟶',
    '1' => '𝟷',
    '2' => '𝟸',
    '3' => '𝟹',
    '4' => '𝟺',
    '5' => '𝟻',
    '6' => '𝟼',
    '7' => '𝟽',
    '8' => '𝟾',
    '9' => '𝟿'
)

for (f, map) in (
        :to_blackboardbold => :text_bb_map,
        :to_boldface => :text_bf_map,
        :to_italic => :text_it_map,
        :to_caligraphic => :text_cal_map,
        :to_frakture => :text_frak_map,
        :to_mono => :text_mono_map,
    )
    @eval begin
        $f(x::Char) = get($map, x, x) # default to same char to keep errors down
        $f(io::IO, x::Char) = print(io, $f(x))
        $f(io::IO, x::AbstractString) = for char in x
            print(io, $f(char))
        end
        $f(x::AbstractString) = sprint() do io
            $f(io, x)
        end
    end
end
"""
Underlines a text
"""
function to_underline(x::AbstractString)
    sprint() do io
        to_underline(io, x)
    end
end
function to_underline(io::IO, x::Char)
    print(io, x, Char(0x0332))
end
function to_underline(io::IO, x::AbstractString)
    for c in x
        to_underline(io, c)
    end
end

"""
Overlines a text
"""
function to_overline(x::AbstractString)
    sprint() do io
        to_overline(io, x)
    end
end
function to_overline(io::IO, x::Char)
    print(io, x, Char(0x0305))
end
function to_overline(io::IO, x::AbstractString)
    for c in x
        to_overline(io, c)
    end
end


text:
module SimpleDataStructures

export
    SimpleVector,
    SimpleSparseVector,
    SimpleDict,
    SimpleDefaultDict,
    SimpleSet,
    SimpleLinkedList,
    SimpleStack,
    SimpleQueue

include("vector.jl")
include("sparse_vector.jl")
include("dict.jl")
include("default_dict.jl")
include("set.jl")
include("doubly_linked_list.jl")
include("stack.jl")
include("queue.jl")

end