text:
@testset "Inverse2Reification: $(fct_type), dimension $(dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], dim in [2, 3], T in [Int, Float64]
    base_model = if T == Int
        IntReificationEqualToModel{T}()
    elseif T == Float64
        FloatReificationEqualToModel{T}()
    else
        @assert false
    end
    mock = MOIU.MockOptimizer(base_model)
    model = COIB.Inverse2Reification{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Reification{MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Inverse,
    )

    if T == Int
        x, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
        y, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
    elseif T == Float64
        x = MOI.add_variables(model, dim)
        y = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x..., y...])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x..., y...])
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.Inverse(dim))

    for i in 1:dim
        @test MOI.is_valid(model, x[i])
        @test MOI.is_valid(model, y[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.Inverse}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Inverse) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}),
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 2 * dim ^ 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == 2 * dim ^ 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}()) == 2 * dim ^ 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == dim ^ 2

        @test Set(MOI.get(bridge, MOI.ListOfVariableIndices())) == Set([bridge.vars_first..., bridge.vars_second...])
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}())) == Set([bridge.vars_first_bin..., bridge.vars_second_bin...])
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}())) == Set([bridge.cons_first_reif..., bridge.cons_second_reif...])
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == bridge.cons_equivalence
    end

    @testset "Set of variables" begin
        @test length(bridge.vars_first) == dim ^ 2
        @test length(bridge.vars_second) == dim ^ 2

        for i in 1:dim
            for j in 1:dim
                @test MOI.is_valid(model, bridge.vars_first[i, j])
                @test MOI.is_valid(model, bridge.vars_second[i, j])

                @test MOI.is_valid(model, bridge.vars_first_bin[i, j])
                @test MOI.is_valid(model, bridge.vars_second_bin[i, j])

                @test MOI.get(model, MOI.ConstraintFunction(), bridge.vars_first_bin[i, j]) == bridge.vars_first[i, j]
                @test MOI.get(model, MOI.ConstraintSet(), bridge.vars_first_bin[i, j]) == MOI.ZeroOne()

                @test MOI.get(model, MOI.ConstraintFunction(), bridge.vars_second_bin[i, j]) == bridge.vars_second[i, j]
                @test MOI.get(model, MOI.ConstraintSet(), bridge.vars_second_bin[i, j]) == MOI.ZeroOne()
            end
        end
    end

    @testset "Reification" begin
        @test length(bridge.cons_first_reif) == dim ^ 2
        @test length(bridge.cons_second_reif) == dim ^ 2

        for i in 1:dim
            for j in 1:dim
                @test MOI.is_valid(model, bridge.cons_first_reif[i, j])
                @test MOI.is_valid(model, bridge.cons_second_reif[i, j])

                f1 = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_first_reif[i, j])
                @test length(f1.terms) == 2
                @test f1.constants == T[zero(T), -T(j)]
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_first_reif[i, j]) == CP.Reification(MOI.EqualTo(zero(T)))

                t1 = f1.terms[1]
                @test t1.output_index == 1
                @test t1.scalar_term.coefficient === one(T)
                @test t1.scalar_term.variable == bridge.vars_first[i, j]

                t2 = f1.terms[2]
                @test t2.output_index == 2
                @test t2.scalar_term.coefficient === one(T)
                @test t2.scalar_term.variable == x[i]

                f2 = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_second_reif[i, j])
                @test length(f2.terms) == 2
                @test f2.constants == T[zero(T), -T(i)]
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_second_reif[i, j]) == CP.Reification(MOI.EqualTo(zero(T)))

                t1 = f2.terms[1]
                @test t1.output_index == 1
                @test t1.scalar_term.coefficient === one(T)
                @test t1.scalar_term.variable == bridge.vars_second[i, j]

                t2 = f2.terms[2]
                @test t2.output_index == 2
                @test t2.scalar_term.coefficient === one(T)
                @test t2.scalar_term.variable == y[j]
            end
        end
    end

    @testset "Equivalence" begin
        @test length(bridge.cons_equivalence) == dim ^ 2

        for i in 1:dim
            for j in 1:dim
                @test MOI.is_valid(model, bridge.cons_equivalence[i, j])
                f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_equivalence[i, j])
                @test length(f.terms) == 2
                @test f.constant == zero(T)
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_equivalence[i, j]) == MOI.EqualTo(zero(T))

                t1 = f.terms[1]
                @test t1.coefficient === one(T)
                @test t1.variable == bridge.vars_first[i, j]

                t2 = f.terms[2]
                @test t2.coefficient === -one(T)
                @test t2.variable == bridge.vars_second[i, j]
            end
        end
    end
end


text:
using Weave
using DrWatson

# for f in readdir(joinpath(@__DIR__, "../src/weave/RCall/"), join=true)
#     weave(f, outpath="../notebooks/RCall")
# end

# for f in readdir(joinpath(@__DIR__, "../src/weave/RCall/"), join=true)
#     notebook(f, outpath="../notebooks/RCall")
# end

nbname(f) = splitext(basename(f))[1] * ".ipynb"

for pkg in ["Turing", "Soss", "RCall"]
    mkpath(projectdir("notebooks", pkg))
    cd(projectdir("notebooks", pkg)) do
        for f in readdir(srcdir("weave", pkg), join=true)
            convert_doc(f, nbname(f))
        end
    end
end

# for pkg in ["RCall", "Soss"]
#     for f in readdir(joinpath(@__DIR__, "../src/weave", pkg), join=true)
#         weave(f, out_path=joinpath(@__DIR__, "../html", pkg))
#     end
# end

ext: jl


text:
using NNN
using Base.Test
using FactCheck

include("sample.jl")
include("rnn.jl")


text:
@system FloralInitiation(Stage, Germination) begin
    critPPD: critical_photoperiod => 12.5 ~ preserve(u"hr", parameter)

    floral_initiateable(germinated) ~ flag

    #FIXME: implement Sun
    floral_initiated(critPPD, day_length=sun.day_length, day=sun.day) => begin
        #FIXME solstice consideration is broken (flag turns false after solstice) and maybe unnecessary
        # w = self.pheno.weather
        # solstice = w.time.tz.localize(datetime.datetime(w.time.year, 6, 21))
        # # no MAX_LEAF_NO implied unlike original model
        # w.time <= solstice and w.day_length >= self.critical_photoperiod
        day_length >= critPPD && day <= 171u"d"
    end ~ flag(once)

    floral_initiating(floral_initiateable & !floral_initiated) ~ flag

    # #FIXME postprocess similar to @produce?
    # def finish(self):
    #     GDD_sum = self.pheno.gdd_recorder.rate
    #     print(f"* Floral initiation: time = {self.time}, GDDsum = {GDD_sum}")
end


text:
"""
    ParseType

- PARSE_NONE: Only collect information
- PARSE_MINIMAL: Parse only useful information.
- PARSE_OPPONENT: Parse all information of the opponent
- PARSE_BOTH_PLAYERS: Parse all information for both players.
"""
@enum ParseType begin
    PARSE_NONE = 0
    PARSE_MINIMAL = 1
    PARSE_OPPONENT = 2
    PARSE_BOTH_PLAYERS = 3
end
VsRecorderBase.enum_prefix(::Type{ParseType}) = "PARSE_"

"""
    PokemonBattle

Only works for double ranked battles right now.
"""
Base.@kwdef struct PokemonBattle <: AbstractVsSource
    language::GameLanguage
    double::Bool
    # (height, width)
    image_size::Tuple{Int, Int}
    parse_type::ParseType
end

const GlobalVsConfig = Ref{Nullable{
    VsConfig{DefaultStrategy, PokemonBattle}
}}(nothing)

function VsRecorderBase.vs_setup(
    ::Type{PokemonBattle};
    language::GameLanguage = get_game_language(),
    double = true,
    num_skip_frames = 59,
    use_gray_image = false,
    feature_size = (360, 640),
    image_size = (960, 540),
    parse_type::ParseType = PARSE_MINIMAL,
    match_threshold = 0.1,
    set_global = true,
    init_descriptors = true
)

if use_gray_image
    @warn "Using gray image is strongly unrecommended for its low accuracy."
end

set_language!(GlobalI18nContext[], enum_to_string(language), String[])
ocr_language = VsI18n.OCR_LANGUAGES[language]
strategy = DefaultStrategy(
    match_threshold = match_threshold,
    height = feature_size[1],
    width = feature_size[2],
    init_descriptors = init_descriptors
)
source = PokemonBattle(
    language = language,
    double = double,
    image_size = image_size,
    parse_type = parse_type
)

config = VsConfig(
    num_skip_frames = num_skip_frames,
    use_gray_image = use_gray_image,
    ocr_language = ocr_language,
    strategy = strategy,
    source = source
)

if set_global
    GlobalVsConfig[] = config
end

config

end

default_config() = isnothing(GlobalVsConfig[]) ? vs_setup(
    PokemonBattle, parse_type = PARSE_BOTH_PLAYERS,
    set_global = false,
    init_descriptors = false
) : GlobalVsConfig[]

default_context() = VsRecorderBase.initialize(default_config())


text:
# This file was autogenerated in src/version_git.jl
immutable GitVersionInfo
    commit::AbstractString
    commit_short::AbstractString
    branch::AbstractString
    build_number::Int
    date_string::AbstractString
    tagged_commit::Bool
    fork_master_distance::Int
    fork_master_timestamp::Float64
end

const GIT_VERSION_INFO = GitVersionInfo(
    "8e92541a5f5c7c3b2eae117a0b9daa471009d3a6",
    "8e92541*",
    "master",
    1,
    "2015-02-13 11:47 UTC",
    false,
    0,
    1423828038.
)

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using Base.Test, MathOptInterface
using LinQuadOptInterface

const MOI= MathOptInterface
const MOIT = MathOptInterface.Test
const LQOI = LinQuadOptInterface

@testset "LinQuadOptInterface" begin
    @testset "Unit Tests" begin
        config = MOIT.TestConfig(solve = false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.basic_constraint_tests(solver, config;
            exclude = [
                (MOI.SingleVariable, MOI.EqualTo{Float64}),
                (MOI.SingleVariable, MOI.Integer),
                (MOI.SingleVariable, MOI.LessThan{Float64}),
                (MOI.SingleVariable, MOI.Interval{Float64}),
                (MOI.SingleVariable, MOI.GreaterThan{Float64})
            ]
        )
        MOIT.unittest(solver, config, [
            "solve_affine_interval",
            "solve_qp_edge_cases",
            "solve_qcp_edge_cases",
            "solve_affine_deletion_edge_cases"
        ])
    end

    @testset "Linear tests" begin
        linconfig = MOIT.TestConfig(solve = false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.contlineartest(solver , linconfig)
    end

    @testset "Quadratic tests" begin
        quadconfig = MOIT.TestConfig(solve=false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.contquadratictest(solver, quadconfig)
    end

    @testset "Linear Conic tests" begin
        linconfig = MOIT.TestConfig(solve=false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.lintest(solver, linconfig)
    end

    @testset "Integer Linear tests" begin
        intconfig = MOIT.TestConfig(solve=false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.intlineartest(solver, intconfig, ["int2"])
    end

    @testset "ModelLike tests" begin
        intconfig = MOIT.TestConfig()
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.nametest(solver)
        @testset "validtest" begin
            MOIT.validtest(solver)
        end
        @testset "emptytest" begin
            MOIT.emptytest(solver)
        end
        @testset "orderedindicestest" begin
            # MOIT.orderedindicestest(solver)
        end
        @testset "canaddconstrainttest" begin
            MOIT.canaddconstrainttest(solver, Float64, Complex{Float64})
        end
        @testset "copytest" begin
            solver2 = LQOI.MockLinQuadOptimizer()
            MOIT.copytest(solver,solver2)
        end
    end
end


text:
module CLRBridge

import ...DotNET:HRESULT,isfailed,CLRHostError,BStr,CLRHost,create_host,create_delegate
export CLRObject,null,isnull,gethandle,CLRException,BindingFlags

const Handle = UInt

mutable struct CLRObject
    handle::Handle
end

gethandle(obj::CLRObject) = getfield(obj, :handle)

const null = CLRObject(0)

isnull(obj::CLRObject) = gethandle(obj) == 0

function track(handle)
    obj = CLRObject(handle)
    if handle != 0
        finalizer(obj) do _
            Release(handle)
        end
    end
    obj
end

Base.string(x::CLRObject) = GetString(gethandle(x))

struct CLRException <: Exception
    message::String
    object::CLRObject
end

function Base.showerror(io::IO, ex::CLRException)
    print(io, "CLRException: $(ex.message)")
end

function track_and_throw(exhandle)
    exobj = track(exhandle)
    msg = try
        GetString(exhandle)
    catch ex
        "Failed to get exception description, caused by $ex"
    end
    throw(CLRException(msg, exobj))
end

empty_fp() = Ref{Ptr{Cvoid}}(0)

const fp_Release = empty_fp()
const fp_Duplicate = empty_fp()
const fp_CreateArray = empty_fp()
const fp_PutString = empty_fp()
const fp_GetString = empty_fp()
const fp_FreeString = empty_fp()
const fp_PutChar = empty_fp()
const fp_GetChar = empty_fp()
const fp_PutBool = empty_fp()
const fp_GetBool = empty_fp()
const fp_PutInt8 = empty_fp()
const fp_GetInt8 = empty_fp()
const fp_PutUInt8 = empty_fp()
const fp_GetUInt8 = empty_fp()
const fp_PutInt16 = empty_fp()
const fp_GetInt16 = empty_fp()
const fp_PutUInt16 = empty_fp()
const fp_GetUInt16 = empty_fp()
const fp_PutInt32 = empty_fp()
const fp_GetInt32 = empty_fp()
const fp_PutUInt32 = empty_fp()
const fp_GetUInt32 = empty_fp()
const fp_PutInt64 = empty_fp()
const fp_GetInt64 = empty_fp()
const fp_PutUInt64 = empty_fp()
const fp_GetUInt64 = empty_fp()
const fp_PutFloat32 = empty_fp()
const fp_GetFloat32 = empty_fp()
const fp_PutFloat64 = empty_fp()
const fp_GetFloat64 = empty_fp()

const fp_GetType = empty_fp()
const fp_GetObjectType = empty_fp()
const fp_InvokeMember = empty_fp()

const fp_SetCallbackHandler = empty_fp()
const fp_CreateDelegate = empty_fp()

const registered_callbacks = Dict{UInt,Function}()

function handle_callback(context, argc)
    if !haskey(registered_callbacks, context)
        error("Callback function not registered")
    end
    fn = registered_callbacks[context]
    return Handle(fn(argc))
end

function init(host::CLRHost)
    fp_primitive(x) = create_delegate(host, "CLRBridge", "CLRBridge.Primitive", x)
    fp_Release[] = fp_primitive("Release")
    fp_Duplicate[] = fp_primitive("Duplicate")
    fp_CreateArray[] = fp_primitive("CreateArray")
    fp_PutString[] = fp_primitive("PutString")
    fp_GetString[] = fp_primitive("GetString")
    fp_FreeString[] = fp_primitive("FreeString")
    fp_PutChar[] = fp_primitive("PutChar")
    fp_GetChar[] = fp_primitive("GetChar")
    fp_PutBool[] = fp_primitive("PutBool")
    fp_GetBool[] = fp_primitive("GetBool")
    fp_PutInt8[] = fp_primitive("PutInt8")
    fp_GetInt8[] = fp_primitive("GetInt8")
    fp_PutUInt8[] = fp_primitive("PutUInt8")
    fp_GetUInt8[] = fp_primitive("GetUInt8")
    fp_PutInt16[] = fp_primitive("PutInt16")
    fp_GetInt16[] = fp_primitive("GetInt16")
    fp_PutUInt16[] = fp_primitive("PutUInt16")
    fp_GetUInt16[] = fp_primitive("GetUInt16")
    fp_PutInt32[] = fp_primitive("PutInt32")
    fp_GetInt32[] = fp_primitive("GetInt32")
    fp_PutUInt32[] = fp_primitive("PutUInt32")
    fp_GetUInt32[] = fp_primitive("GetUInt32")
    fp_PutInt64[] = fp_primitive("PutInt64")
    fp_GetInt64[] = fp_primitive("GetInt64")
    fp_PutUInt64[] = fp_primitive("PutUInt64")
    fp_GetUInt64[] = fp_primitive("GetUInt64")
    fp_PutFloat32[] = fp_primitive("PutFloat32")
    fp_GetFloat32[] = fp_primitive("GetFloat32")
    fp_PutFloat64[] = fp_primitive("PutFloat64")
    fp_GetFloat64[] = fp_primitive("GetFloat64")
    fp_meta(x) = create_delegate(host, "CLRBridge", "CLRBridge.Meta", x)
    fp_GetType[] = fp_meta("GetType")
    fp_GetObjectType[] = fp_meta("GetObjectType")
    fp_InvokeMember[] = fp_meta("InvokeMember")
    fp_callback(x) = create_delegate(host, "CLRBridge", "CLRBridge.Callback", x)
    fp_SetCallbackHandler[] = fp_callback("SetCallbackHandler")
    fp_CreateDelegate[] = fp_callback("CreateDelegate")
    handler = @cfunction(handle_callback,Handle,(UInt, UInt32))
    SetCallbackHandler(handler)
    nothing
end

function Release(handle)
    ccall(fp_Release[], Bool, (Handle,), handle)
end

function Duplicate(handle)
    track(ccall(fp_Duplicate[], Handle, (Handle,), handle))
end

function CreateArray(argc)
    track(ccall(fp_CreateArray[], Handle, (UInt32,), argc))
end

function PutString(handle, value)
    ccall(fp_PutString[], Handle, (Handle, BStr), handle, value)
end

function GetString(handle)
    bstr = ccall(fp_GetString[], BStr, (Handle,), handle)
    s = unsafe_string(bstr)
    FreeString(bstr)
    return s
end

function FreeString(bstr)
    ccall(fp_FreeString[], Cvoid, (BStr,), bstr)
end

function PutChar(handle, value)
    ccall(fp_PutChar[], Handle, (Handle, UInt16), handle, value)
end

function GetChar(handle)
    Char(ccall(fp_GetChar[], UInt16, (Handle,), handle))
end

function PutBool(handle, value)
    ccall(fp_PutBool[], Handle, (Handle, Bool), handle, value)
end

function GetBool(handle)
    ccall(fp_GetBool[], Bool, (Handle,), handle)
end

function PutInt8(handle, value)
    ccall(fp_PutInt8[], Handle, (Handle, Int8), handle, value)
end

function GetInt8(handle)
    ccall(fp_GetInt8[], Int8, (Handle,), handle)
end

function PutUInt8(handle, value)
    ccall(fp_PutUInt8[], Handle, (Handle, UInt8), handle, value)
end

function GetUInt8(handle)
    ccall(fp_GetUInt8[], UInt8, (Handle,), handle)
end

function PutInt16(handle, value)
    ccall(fp_PutInt16[], Handle, (Handle, Int16), handle, value)
end

function GetInt16(handle)
    ccall(fp_GetInt16[], Int16, (Handle,), handle)
end

function PutUInt16(handle, value)
    ccall(fp_PutUInt16[], Handle, (Handle, UInt16), handle, value)
end

function GetUInt16(handle)
    ccall(fp_GetUInt16[], UInt16, (Handle,), handle)
end

function PutInt32(handle, value)
    ccall(fp_PutInt32[], Handle, (Handle, Int32), handle, value)
end

function GetInt32(handle)
    ccall(fp_GetInt32[], Int32, (Handle,), handle)
end

function PutUInt32(handle, value)
    ccall(fp_PutUInt32[], Handle, (Handle, UInt32), handle, value)
end

function GetUInt32(handle)
    ccall(fp_GetUInt32[], UInt32, (Handle,), handle)
end

function PutInt64(handle, value)
    ccall(fp_PutInt64[], Handle, (Handle, Int64), handle, value)
end

function GetInt64(handle)
    ccall(fp_GetInt64[], Int64, (Handle,), handle)
end

function PutUInt64(handle, value)
    ccall(fp_PutUInt64[], Handle, (Handle, UInt64), handle, value)
end

function GetUInt64(handle)
    ccall(fp_GetUInt64[], UInt64, (Handle,), handle)
end

function PutFloat32(handle, value)
    ccall(fp_PutFloat32[], Handle, (Handle, Float32), handle, value)
end

function GetFloat32(handle)
    ccall(fp_GetFloat32[], Float32, (Handle,), handle)
end

function PutFloat64(handle, value)
    ccall(fp_PutFloat64[], Handle, (Handle, Float64), handle, value)
end

function GetFloat64(handle)
    ccall(fp_GetFloat64[], Float64, (Handle,), handle)
end

function GetType(typename)
    exception = Ref{Handle}()
    ret = ccall(fp_GetType[], Handle, (BStr, Ptr{Handle}), typename, exception)
    if exception[] != 0
        track_and_throw(exception[])
    end
    return track(ret)
end

function GetObjectType(handle)
    ret = ccall(fp_GetObjectType[], Handle, (Handle,), handle)
    return track(ret)
end

baremodule BindingFlags

using Base:@enum
@enum BindingFlag::UInt32 begin
    Default = 0
    IgnoreCase = 1
    DeclaredOnly = 2
    Instance = 4
    Static = 8
    Public = 16
    NonPublic = 32
    FlattenHierarchy = 64
    InvokeMethod = 256
    CreateInstance = 512
    GetField = 1024
    SetField = 2048
    GetProperty = 4096
    SetProperty = 8192
    PutDispProperty = 16384
    PutRefDispProperty = 32768
    ExactBinding = 65536
    SuppressChangeType = 131072
    OptionalParamBinding = 262144
    IgnoreReturn = 16777216
end
end

import .BindingFlags:BindingFlag

Base.:(|)(a::BindingFlag, b::BindingFlag) = UInt32(a) | UInt32(b)
Base.:(|)(a::BindingFlag, b) = UInt32(a) | b
Base.:(|)(a, b::BindingFlag) = b | a

function InvokeMember(type, name, bindingFlags, binder, target, providedArgs)
    exception = Ref{Handle}()
    ret = ccall(fp_InvokeMember[], Handle,
    (Handle, BStr, UInt32, Handle, Handle, Ptr{Handle}, UInt32, Ptr{Handle}),
    type, name, bindingFlags, binder, target, providedArgs, length(providedArgs), exception)
    if exception[] != 0
        track_and_throw(exception[])
    end
    return track(ret)
end

function SetCallbackHandler(fp)
    ccall(fp_SetCallbackHandler[], Cvoid, (Ptr{Cvoid},), fp)
end

function CreateDelegate(htype, context)
    exception = Ref{Handle}()
    ret = ccall(fp_CreateDelegate[], Handle, (Handle, UInt, Ptr{Handle}), htype, context, exception)
    if exception[] != 0
        track_and_throw(exception[])
    end
    return track(ret)
end

end


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
np = pyimport("numpy")

β = 1/4; γ = 1/2
a = b = 0.1
m = 40
n = 20
h = 0.01
NT = 200
Δt = 5/NT
bdedge = []
for j = 1:n
  push!(bdedge, [(j-1)*(m+1)+m+1 j*(m+1)+m+1])
end
bdedge = vcat(bdedge...)

bdnode = Int64[]
for j = 1:n+1
  push!(bdnode, (j-1)*(m+1)+1)
end

M = compute_fem_mass_matrix1(m, n, h)
S = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = [M S;S M]

H = [1.0 0.0 0.0
    0.0 1.0 0.0
    0.0 0.0 0.5]
K = compute_fem_stiffness_matrix(H, m, n, h)
C = a*M + b*K # damping matrix

L = M + γ*Δt*C + β*Δt^2*K
L, Lbd = fem_impose_Dirichlet_boundary_condition(L, bdnode, m, n, h)

a = zeros(2(m+1)*(n+1))
v = zeros(2(m+1)*(n+1))
d = zeros(2(m+1)*(n+1))
U = zeros(2(m+1)*(n+1),NT+1)

Sigma = zeros(NT+1, 4m*n, 3)
Varepsilon = zeros(NT+1, 4m*n, 3)
for i = 1:NT
    global a, v, d
    T = eval_f_on_boundary_edge((x,y)->0.1, bdedge, m, n, h)
    T = [T zeros(length(T))]
    rhs = compute_fem_traction_term(T, bdedge, m, n, h)
    if i*Δt>3.0
      rhs = zero(rhs)
    end

    td = d + Δt*v + Δt^2/2*(1-2β)*a
    tv = v + (1-γ)*Δt*a
    rhs = rhs - C*tv - K*td
    rhs[[bdnode; bdnode.+(m+1)*(n+1)]] .= 0.0

    a = L\rhs
    d = td + β*Δt^2*a
    v = tv + γ*Δt*a
    U[:,i+1] = d

    Varepsilon[i+1,:,:] = eval_strain_on_gauss_pts(U[:,i+1], m, n, h)
    Sigma[i+1,:,:] = Varepsilon[i+1,:,:] * H
end

visualize_displacement(U, m, n, h; name = "_linear", xlim_=[-0.01,0.5], ylim_=[-0.05,0.22])
visualize_displacement(U, m, n, h;  name = "_linear")
visualize_stress(H, U, m, n, h;  name = "_linear")

close("all")
figure(figsize=(15,5))
subplot(1,3,1)
idx = div(n,2)*(m+1) + m+1
plot((0:NT)*Δt, U[idx,:])
xlabel("time")
ylabel("x displacement")

subplot(1,3,2)
idx = 4*(div(n,2)*m + m)
plot((0:NT)*Δt, Sigma[:,idx,1])
xlabel("time")
ylabel("x stress")

subplot(1,3,3)
idx = 4*(div(n,2)*m + m)
plot((0:NT)*Δt, Varepsilon[:,idx,1])
xlabel("time")
ylabel("x strain")
savefig("linear.png")


text:
"""
It represents a FE function at a set of time steps. It is a wrapper of a ODE
solution for free values combined with data for Dirichlet values. Thus, it is a
lazy iterator that computes the solution at each time step when accessing the
solution.
"""
struct TransientFESolution
  odesol::ODESolution
  trial
end

function TransientFESolution(
  solver::TransientFESolver, op::TransientFEOperator, uh0, t0::Real, tF::Real)

  ode_solver = solver.odes
  ode_op = get_algebraic_operator(op)
  u0 = get_free_values(uh0)
  ode_sol = solve(ode_solver,ode_op,u0,t0,tF)
  trial = get_trial(op)

  TransientFESolution(ode_sol, trial)
end

#@fverdugo this is a general implementation of iterate for TransientFESolution
# We could also implement another one for the very common case that the
# underlying ode_op is a ODEOpFromFEOp object

function Base.iterate(sol::TransientFESolution)

  odesolnext = Base.iterate(sol.odesol)

  if odesolnext === nothing
    return nothing
  end

  (uf, tf), odesolstate = odesolnext

  Uh = allocate_trial_space(sol.trial)
  Uh = evaluate!(Uh,sol.trial,tf)
  uh = FEFunction(Uh,uf)

  state = (Uh, odesolstate)

  (uh, tf), state
end

function Base.iterate(sol::TransientFESolution, state)

  Uh, odesolstate = state

  odesolnext = Base.iterate(sol.odesol,odesolstate)

  if odesolnext === nothing
    return nothing
  end

  (uf, tf), odesolstate = odesolnext

  Uh = evaluate!(Uh,sol.trial,tf)
  uh = FEFunction(Uh,uf)

  state = (Uh, odesolstate)

  (uh, tf), state

end

function test_transient_fe_solution(fesol::TransientFESolution)
  for (uhn,tn) in fesol
    @test isa(uhn,FEFunction)
    @test isa(tn,Real)
  end
  true
end


text:
matrixdepot()
groups = ["symmetric", "inverse", "ill-cond", "pos-def", "eigen","sparse", "random", "regprob", "all", "data"]

for group in groups
    matrixdepot(group)
end

try
    matrixdepot("something")
catch ArgumentError
    println("matrixdepot: no information for `something`")
end

n = rand(1:55)
name = matrixdepot(n)
matrixdepot(name)

list = matrixdepot(1, 3, 4:20)
m = length(matrixdepot("all"))
try
    matrixdepot(m+1)
catch ArgumentError
    println("error: access the $(m+1) matrix")
end

@addgroup newlist = matrixdepot(3:6, 20)

workspace()
using MatrixDepot
using Base.Test

println(matrixdepot("newlist"))

@rmgroup newlist


text:
export plot_boundarymap, plot_boundarymap_portion

"""
    plot_boundarymap(bmap, intervals; kwargs...)

Plots the boundary map.
The input arguments are the return values of [`boundarymap`](@ref) (also applies
for the parallelized version, where `bmap` is a vector of vectors of 2-vectors).

## Keyword Arguments
* `ax = (figure(); gca())` : The axis to plot on.
* `color = "C0"` : The color to use for the plotted points. Can be either a
  color for `PyPlot.plot` or a vector of colors of length `length(bmap)`, in
  order to give each initial condition a different color (for parallelized version).
* `ms = 1.0` : Marker size of the points.
* `bordercolor = "C3"` : The color of the vertical lines that denote the obstacle
  borders.
* `obstacleindices = true`: show obstacle indices above plot
* Any other keyword argument is passed to `PyPlot.plot` which plots the points of
  the section.
"""
function plot_boundarymap(bmap, intervals;
    ax = (PyPlot.figure(); PyPlot.gca()),
    color = "C0", bordercolor = "C3", ms = 1.0, obstacleindices = true, kwargs...)

    # Plot the points
    if typeof(bmap) <: Vector{<:SVector}
        c = typeof(color) <: AbstractVector ? color[1] : color
        ξs = [b[1] for b in bmap]; sφs = [b[2] for b in bmap]

        ax.plot(ξs, sφs; marker="o", color = c,
        linestyle="None", ms = ms, kwargs...)
    else
        for (i, bmapp) in enumerate(bmap)
            ξs = [b[1] for b in bmapp]; sφs = [b[2] for b in bmapp]

            c = typeof(color) <: AbstractVector ? color[i] : color

            ax.plot(ξs, sφs; marker="o", color = c,
            linestyle="None", ms = ms, kwargs...)
        end
    end

    # Plot obstacle limits
    xmax = 0.0
    for intv ∈ intervals
        for xval ∈ intv
            ax.plot([xval,xval], [-1, 1], linewidth = 1.5, color = bordercolor,
            alpha = 0.5)
            xmax = (xval > xmax) ? xval : xmax
        end
    end
    ax.set_xlim(0,xmax)
    ax.set_ylim(-1,1)
    ax.set_xlabel("\$\\xi\$")
    ax.set_ylabel("\$\\sin(\\phi_n)\$")

    #number obstacles by index
    if obstacleindices
        #introduce twin axis
        ax2 = ax.twiny()

        #non-labelled major tics at every obstacle border
        ax2.xaxis.set_major_formatter(PyPlot.matplotlib.ticker.NullFormatter())
        ax2.set_xticks(intervals)

        #zero-length minor tics in between borders, labelled with the obstacle index
        ax2.xaxis.set_minor_locator(PyPlot.matplotlib.ticker.FixedLocator(
            [(intervals[i] + intervals[i+1])/2 for i in 1:length(intervals)-1]
        ))

        ax2.xaxis.set_minor_formatter(
            PyPlot.matplotlib.ticker.FixedFormatter(map(x->string(x),1:length(intervals)))
        )

        ax2.tick_params(axis="x", which="minor", length=0)
        ax2.set_xlabel("obstacle index")

        return (ax, ax2)
    end
    return ax
end

"""
    plot_bmportion(d, δξ, δφ = δξ; kwargs...)

Plots histograms in boundary coordinates.
The input arguments are the dictionary `d` generated by [`boundarymap_portion`](@ref)
and the histogram box size `δ`

## Keyword Arguments
* `ax = (figure(); gca())` : The axis to plot on.
* `cb = true` : show a colorbar next to the plot
* `transp = true`: don't plot non-visited boxes
* Any other keyword argument is passed to `PyPlot.pcolormesh` which plots the histogram.
"""
function plot_boundarymap_portion(d, δξ, δφ = δξ;
    ax = (PyPlot.figure(); PyPlot.gca()), cb::Bool = true,
                        transp::Bool = true, kwargs...)
    ξmax = maximum(x->x[1], keys(d)) + 2
    φmax = maximum(x->x[2], keys(d)) + 2

    data = ones((ξmax, φmax)).* (transp ? NaN : 1)

    for (key, val) ∈ d
        data[key[1] + 1, key[2] + 1] = val
    end

    ξs = repeat(δξ.*(0:ξmax-1), 1, φmax)
    φs = repeat((δφ.*(0:1:φmax-1) .- 1)',ξmax, 1)

    plot = ax.pcolormesh(ξs, φs, data; kwargs...)

    cb && PyPlot.colorbar(plot)
    return ax
end


text:
using Plots

long_name = Dict("T" =>"Temperature [°C]", "wT"=>"Temperature flux [°C⋅m/s]")
x_lims = Dict("T" =>(18,20), "wT"=>(-1e-5,4e-5))

"""
plot_profile(gp::GP, data::ProfileData, V_name, time_index, gpr_prediction)
# Description
Used in interact.jl
Plots the simulation profile at a single given index in the data timeseries along with the
corresponding GP prediction (prediciting on the training data). i.e. produces a snapshot of the profile evolution.
The gpr_prediction is computed outside this function for efficiency.
# Arguments
- `ℳ` (GP or NN). The model.
- `𝒟` (ProfileData). The ProfileData object used for training and testing.
- `V_name` (ProfileData). The ProfileData object whose starting profile will be evolved forward using ℳ.
- `time_index` (Int). The time index
- `gpr_prediction` (Array). Output of get_gpr_pred (which should only be computed once) on ℳ and 𝒟.
"""
function plot_profile(ℳ, 𝒟::ProfileData, time_index, gpr_prediction)
    day_string = string(floor(Int, 𝒟.t[time_index]/86400))
    variable = 𝒟.problem.variable # "T" or "wT"
    xlims = x_lims[variable]

    exact = 𝒟.v[:,time_index]

    p = scatter(gpr_prediction[time_index], 𝒟.zavg, label = "GP", xlims=xlims)
    plot!(exact, 𝒟.z, legend = :topleft, label = "LES", xlabel = "$(long_name[variable])", ylabel = "Depth [m]", title = "day " * day_string, xlims=xlims)
    return p
end

function animate_profile(ℳ, 𝒟)

    variable = 𝒟.problem.variable # "T" or "wT"
    xlims = x_lims[variable]

    predi = predict(ℳ, 𝒟; postprocessed=true)

    animation_set = 1:30:(𝒟.Nt-2)
    anim = @animate for i in animation_set
        exact = 𝒟.v[:,i]
        day_string = string(floor(Int, 𝒟.t[i]/86400))
        scatter(predi[i], 𝒟.zavg, label = "GP")
        plot!(exact, 𝒟.z, legend = :topleft, label = "LES", xlabel = "$(long_name[variable])", ylabel = "Depth [m]", title = "day " * day_string, xlims=xlims)
    end

    return anim
end

# Instead of plotting the full profile computed from the model output, returns the model output directly instead.
# For example, if the problem is a residual problem, plots the residual predicted by the model, not the profile computed from the residual.
function plot_model_output(ℳ, 𝒟, time_index, model_outputs)

    day_string = string(floor(Int, 𝒟.t[time_index]/86400))
    variable = 𝒟.problem.variable # "T" or "wT"
    xlims2 = (minimum(minimum(𝒟.y)), maximum(maximum(𝒟.y)))

    exact = 𝒟.y[time_index]

    p = scatter(model_outputs[time_index], 𝒟.zavg, label = "GP", xlims=xlims)
    plot!(exact, 𝒟.zavg, legend = false, label = "LES", xlabel = "$(long_name[variable])", ylabel = "Depth [m]", title = "Direct model output", xlims=xlims)
    return p
end

# v = 𝒟.v
# v = [v[:,i] for i in size(v)[2]]
# _ , exact_y = get_predictors_targets(v, 𝒟.problem) # what the model output should be exactly.
function animate_profile_and_model_output(ℳ, 𝒟)

    variable = 𝒟.problem.variable # "T" or "wT"
    xlims1 = x_lims[variable]
    # xlims2 = (minimum(minimum(𝒟.y))-0.005, maximum(maximum(𝒟.y))+0.02)
    xlims2 = (minimum(minimum.(𝒟.y)), maximum(maximum.(𝒟.y)))

    model_output, predi = predict(ℳ, 𝒟; postprocessed="both")

    ###
    # f = closure_free_convection_kpp(length(𝒟.z), 𝒟.t[2]-𝒟.t[1], get_les_data("general_strat_32_profiles.jld2")) # full resolution KPP prediction
    # kpp_data = f(OceanTurb.KPP.Parameters())
    # scaling = 𝒟.problem.scaling
    # v = scale(𝒟.v, scaling)
    # kpp_data = scale(kpp_data, scaling) # kpp(T[i])
    # target2 = v .- kpp_data
    ###

    animation_set = 1:30:(𝒟.Nt-2)
    anim = @animate for i in animation_set
        day_string = string(floor(Int, 𝒟.t[i]/86400))

        exact = 𝒟.v[:,i]
        p1 = scatter(predi[i], 𝒟.zavg, label = "GP mean prediction");
        plot!(exact, 𝒟.z, legend = :topleft, label = "LES", xlabel = "$(long_name[variable])", ylabel = "Depth [m]", title = "day " * day_string, xlims=xlims1);

        # exact = target2[:,i]
        exact = 𝒟.y[i]
        p2 = scatter(model_output[i], 𝒟.zavg, label = "Model output");
        plot!(exact, 𝒟.zavg, label = "Target", ylabel = "Depth [m]", title = "Direct model output", xlims=xlims2);

        plot(p1, p2, layout=(@layout [a b]), size=(1000,500));
    end

    anim
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
# Autogenerated wrapper script for CompilerSupportLibraries_jll for x86_64-apple-darwin-libgfortran5
export libatomic, libgcc_s, libgfortran, libgomp, libstdcxx

JLLWrappers.@generate_wrapper_header("CompilerSupportLibraries")
JLLWrappers.@declare_library_product(libatomic, "@rpath/libatomic.1.dylib")
JLLWrappers.@declare_library_product(libgcc_s, "libgcc_s.1.dylib")
JLLWrappers.@declare_library_product(libgfortran, "@rpath/libgfortran.5.dylib")
JLLWrappers.@declare_library_product(libgomp, "@rpath/libgomp.1.dylib")
JLLWrappers.@declare_library_product(libstdcxx, "@rpath/libstdc++.6.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libatomic,
        "lib/libatomic.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgcc_s,
        "lib/libgcc_s.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgfortran,
        "lib/libgfortran.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgomp,
        "lib/libgomp.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libstdcxx,
        "lib/libstdc++.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#TODO: more detailed documentation?
"""
BIPARTITE MATCHING
------------------
    return a maximum weight bipartite matching of a graph

Functions
---------
- Matching_Setup = bipartite_matching_setup{T}(A::SparseMatrixCSC{T,Int64})
- Matching_Setup = bipartite_matching_setup{T}(x::Vector{T},ei::Vector{Int64},ej::Vector{Int64},m::Int64,n::Int64)
- Matching_Output = bipartite_matching_primal_dual{T}(rp::Vector{Int64}, ci::Vector{Int64},ai::Vector{T}, m::Int64, n::Int64)
- Matching_Output = bipartite_matching{T}(A::SparseMatrixCSC{T,Int64})
- Matching_Output = bipartite_matching{T}(w::Vector{T},ei::Vector{Int64},ej::Vector{Int64},m::Int64,n::Int64)
- Matching_Output = bipartite_matching{T}(w::Vector{T},ei::Vector{Int64},ej::Vector{Int64})
- ind = bipartite_matching_indicator{T}(w::Vector{T},ei::Vector{Int64},ej::Vector{Int64})
- (m1,m2) = edge_list(M_output::Matching_output)
- S = create_sparse(M_output::Matching_output)\n
You can check the documentation of each of the output modifiers functions separately.

Example
-------
~~~
W = sprand(10,8,0.5)
bipartite_matching(W)
ei = [1;2;3]
ej = [3;2;4]
Matching_Output = bipartite_matching([10;12;13],ei,ej)
Matching_Output.weight
Matching_Output.cardinality
Matching_Output.match
S = create_sparse(bipartite_matching(W)) # get the sparse matrix
(m1,m2) = edge_list(bipartite_matching(W)) # get the edgelist
~~~
"""

:bipartite_matching

###########################
##    Type Definitions    #
###########################
type Matching_setup
    rp::Array{Int64,1}
    ci::Array{Int64,1}
    ai::Array{Float64,1}
    tripi::Array{Int64,1}
    m::Int64
    n::Int64
end

type Matching_output
    m::Int64
    n::Int64
    weight::Float64
    cardinality::Int64
    match::Array{Int64,1}
end

######################
#   setup  funtions  #
######################

function bipartite_matching_setup_phase1{T}(A::SparseMatrixCSC{T,Int64})
    (nzi,nzj,nzv) = findnz(A)
    return (nzi,nzj,nzv)
end

function bipartite_matching_setup{T}(A::SparseMatrixCSC{T,Int64})
    (nzi,nzj,nzv) = bipartite_matching_setup_phase1(A)
    nedges = length(nzi)
    (m,n)=size(A)

    rp = ones(Int64,m+1) # csr matrix with extra edges
    ci = zeros(Int64,nedges+m)
    ai = zeros(Float64,nedges+m)

    rp[1]=0
    for i=1:nedges
        rp[nzi[i]+1]=rp[nzi[i]+1]+1
    end
    rp=cumsum(rp)

    for i=1:nedges
        ai[rp[nzi[i]]+1]=nzv[i]
        ci[rp[nzi[i]]+1]=nzj[i]
        rp[nzi[i]]=rp[nzi[i]]+1
    end

    for i=1:m # add the extra edges
        ai[rp[i]+1]=0
        ci[rp[i]+1]=n+i
        rp[i]=rp[i]+1
    end

    # restore the row pointer array
    for i=m:-1:1
        rp[i+1]=rp[i]
    end
    rp[1]=0
    rp=rp+1

    #check for duplicates in the data
    colind = zeros(Int64,m+n)
    for i=1:m
        for rpi=rp[i]:rp[i+1]-1
            if colind[ci[rpi]] == 1
                error("bipartite_matching:duplicateEdge")
            end
        colind[ci[rpi]]=1
        end

        for rpi=rp[i]:rp[i+1]-1
            colind[ci[rpi]]=0
        end # reset indicator
    end

    M_setup = Matching_setup(rp,ci,ai,[],m,n)
    return M_setup
end

function bipartite_matching_setup{T}(x::Vector{T},ei::Vector{Int64},
                                        ej::Vector{Int64},m::Int64,n::Int64)
    (nzi,nzj,nzv) = (ei,ej,x)
    nedges = length(nzi)

    rp = ones(Int64,m+1) # csr matrix with extra edges
    ci = zeros(Int64,nedges+m)
    ai = zeros(Float64,nedges+m)
    tripi = zeros(Int64,nedges+m)
    # 1. build csr representation with a set of extra edges from vertex i to
    # vertex m+i

    rp[1]=0
    for i=1:nedges
        rp[nzi[i]+1]=rp[nzi[i]+1]+1
    end

    rp=cumsum(rp)
    for i=1:nedges
        tripi[rp[nzi[i]]+1]=i
        ai[rp[nzi[i]]+1]=nzv[i]
        ci[rp[nzi[i]]+1]=nzj[i]
        rp[nzi[i]]=rp[nzi[i]]+1
    end
    for i=1:m # add the extra edges
        tripi[rp[i]+1]=-1
        ai[rp[i]+1]=0
        ci[rp[i]+1]=n+i
        rp[i]=rp[i]+1
    end

    # restore the row pointer array
    for i=m:-1:1
        rp[i+1]=rp[i]
    end

    rp[1]=0
    rp=rp+1

    #check for duplicates in the data
    colind = zeros(Int64,m+n)

    for i=1:m
        for rpi=rp[i]:rp[i+1]-1
            if colind[ci[rpi]] == 1
                error("bipartite_matching:duplicateEdge")
            end
        colind[ci[rpi]]=1
        end

        for rpi=rp[i]:rp[i+1]-1
            colind[ci[rpi]]=0
        end # reset indicator
    end
    M_setup = Matching_setup(rp,ci,ai,tripi,m,n)
    return M_setup
end

##################
#   primal-dual  #
##################

function bipartite_matching_primal_dual{T}(rp::Vector{Int64}, ci::Vector{Int64},
                    ai::Vector{T}, m::Int64, n::Int64)

    # variables used for the primal-dual algorithm
    alpha=zeros(Float64,m)
    bt=zeros(Float64,m+n)#beta
    queue=zeros(Int64,m)
    t=zeros(Int64,m+n)
    match1=zeros(Int64,m)
    match2=zeros(Int64,m+n)
    tmod = zeros(Int64,m+n)
    ntmod=0

    # initialize the primal and dual variables
    for i=1:m
        for rpi=rp[i]:rp[i+1]-1
            if ai[rpi] > alpha[i]
               alpha[i]=ai[rpi]
            end
        end
    end

    # dual variables (bt) are initialized to 0 already
    # match1 and match2 are both 0, which indicates no matches

    i=1
    while i<=m
        for j=1:ntmod
            t[tmod[j]]=0
        end
        ntmod=0
        # add i to the stack
        head=1
        tail=1
        queue[head]=i
        while head <= tail && match1[i]==0
            k=queue[head]
            for rpi=rp[k]:rp[k+1]-1
                j = ci[rpi]
                if ai[rpi] < alpha[k] + bt[j] - 1e-8
                    continue
                end # skip if tight
                if t[j]==0
                    tail=tail+1
                    if tail <= m
                        queue[tail]=match2[j]
                    end
                    t[j]=k
                    ntmod=ntmod+1
                    tmod[ntmod]=j
                    if match2[j]<1
                        while j>0
                            match2[j]=t[j]
                            k=t[j]
                            temp=match1[k]
                            match1[k]=j
                            j=temp
                        end
                        break
                    end
                end
            end
            head=head+1
        end
        if match1[i] < 1
            theta=Inf
            for j=1:head-1
                t1=queue[j]
                for rpi=rp[t1]:rp[t1+1]-1
                    t2=ci[rpi]
                    if t[t2] == 0 && alpha[t1] + bt[t2] - ai[rpi] < theta
                        theta = alpha[t1] + bt[t2] - ai[rpi]
                    end
                end
            end
            for j=1:head-1
                alpha[queue[j]] = alpha[queue[j]] - theta
            end
            for j=1:ntmod
                bt[tmod[j]] = bt[tmod[j]] + theta
            end
            continue
        end
        i=i+1
    end
    val=0
    for i=1:m
        for rpi=rp[i]:rp[i+1]-1
            if ci[rpi]==match1[i]
                val=val+ai[rpi]
            end
        end
    end
    noute = 0
    for i=1:m
        if match1[i]<=n
            noute=noute+1
        end
    end

    M_output = Matching_output(m,n,val,noute,match1)
    return M_output
end

####################
##    Functions    #
####################

function bipartite_matching{T}(A::SparseMatrixCSC{T,Int64})
    M_setup = bipartite_matching_setup(A)
    return bipartite_matching_primal_dual(M_setup.rp, M_setup.ci, M_setup.ai,
                                          M_setup.m, M_setup.n)
end

function bipartite_matching{T}(w::Vector{T},ei::Vector{Int64},
                                        ej::Vector{Int64},m::Int64,n::Int64)
    M_setup = bipartite_matching_setup(w,ei,ej,m,n)
    return bipartite_matching_primal_dual(M_setup.rp, M_setup.ci, M_setup.ai,
                                          M_setup.m, M_setup.n)
end

function bipartite_matching{T}(w::Vector{T},ei::Vector{Int64},
                                        ej::Vector{Int64})
    M_setup = bipartite_matching_setup(w,ei,ej,maximum(ei),maximum(ej))
    return bipartite_matching_primal_dual(M_setup.rp, M_setup.ci, M_setup.ai,
                                          M_setup.m, M_setup.n)
end

####################
##    Indicator    #
####################

"""
Returns the matching indicator of a matrix stored in triplet format
Example:
bipartite_matching_indicator([10;12;13],[1;2;3],[3;2;4])
"""
function bipartite_matching_indicator{T}(w::Vector{T},ei::Vector{Int64},
                                        ej::Vector{Int64})
    M_setup = bipartite_matching_setup(w, ei, ei, maximum(ei), maximum(ej))
    M_out = bipartite_matching_primal_dual(M_setup.rp, M_setup.ci, M_setup.ai,
                                           M_setup.m, M_setup.n)
    return edge_indicator(M_out,ei, ej)
end

#########################
#   Output Modifiers    #
#########################
"""
Returns the edge list of a matching output
Example:
M_out = bipartite_matching([10;12;13],[1;2;3],[3;2;4])
edge_list(M_out)
"""
function edge_list(M_output::Matching_output)
    m1=zeros(Int64,M_output.cardinality)
    m2=zeros(Int64,M_output.cardinality)
    noute=1
    for i=1:M_output.m
        if M_output.match[i]<=M_output.n
            m1[noute]=i
            m2[noute]=M_output.match[i]
            noute=noute+1
        end
    end
    return (m1,m2)
end

##########
"""
Creates and returns a sparse matrix that represents the outputed matching
Example:
M_out = bipartite_matching([10;12;13],[1;2;3],[3;2;4])
create_sparse(M_out)
"""
function create_sparse(M_output::Matching_output)
    (ei,ej) = edge_list(M_output)
    A = sparse(ei,ej,1,M_output.m,M_output.n)
    return A
end

##########
function edge_indicator(M_output::Matching_output, ei::Vector, ej::Vector)
    assert(length(ei) == length(ej))
    ind = zeros(Int64,length(ei))
    for i=1:length(ei)
        if M_output.match[ei[i]] == ej[i]
            ind[i] = 1
        end
    end
    return ind
end


text:
struct OpString
	idx::String

	function OpString(a::String)
		return new(a)
	end

	function OpString(a::Vector{String})

		if length(a)>1
			res = OpString(a[1]) + OpString(a[2:end])
		else
			res = OpString(a[1])
		end

		return new(res.idx)

	end
end

#
# overwritting base operators for OpString matrix multiplication
#

# in case either a or b is empty "" creat new string by just combining otherwise add a "→" between them
import Base.+
(+)(a::OpString, b::OpString) = a.idx=="" || b.idx=="" ? OpString(string(a.idx, b.idx)) : OpString(string(a.idx,"→ ",b.idx))

# in case either a or b is empty "" the new string is also empty otherwise split operator string and add a ⋅ inbetween
# note that a an b OpStrings have to have same length
import Base.*
(*)(a::OpString, b::OpString) = a.idx == "" || b.idx == "" ? OpString("") : OpString(string.(split(a.idx,"→ "),"⋅", split(b.idx,"→ ")))

import Base.zero
zero(a::OpString) = OpString("")

function StringArray_to_OpStringArray(A::AbstractArray{<:String})

	B = Array{OpString}
	B = fill(OpString(""),size(A)...)
	B[:,:] = OpString.(A[:,:])

	return B
end

function MPO_to_StringArray(mpo_A::MPO)
	bDim = getBondDim(mpo_A)
	bin = 0
	bout = 0
	try
		bin = bDim[1]
		bout = bDim[2]
	catch
		bin = bDim
		bout = bDim
	end

	B = Array{String}
	B = fill("", bin, bout)
	Op_index = mpo_A.Op_index

	Idx = LinearIndices(fill(1,bin,bout))
	mpo_idx = [findfirst(x-> x == idx, Idx) for idx in Op_index]

	B[mpo_idx] .= [string(i) for i in 1:size(mpo_idx)[1]]

	return B
end

MPO_to_OpString(mpoA::MPO) = StringArray_to_OpStringArray(MPO_to_StringArray(mpoA))

# matrix multiplication

function stringM_mul(A::AbstractArray{OpString}, B::AbstractArray{OpString})
	return A*B
end

stringM_mul(A::AbstractArray{String}, B::AbstractArray{String}) = stringM_mul(StringArray_to_OpStringArray(A), StringArray_to_OpStringArray(B))

"""
   convert_OpStringArrayToTuple(A::AbstractArray{<:String})

converts and Array respresenting the MPO product as strings to an Array containing
tuples where each int corresponds to the operator index of the local MPO

#

"""
function convert_OpStringArrayToTuple(A::AbstractArray{OpString})
	bin, bout = size(A)
	B = Array{Any,2}
	B = fill([],bin,bout)
	for i in 1:bin
		row = []
		for j in 1:bout
			try
				#decompose string
				terms = split(A[i,j].idx,"→ ")
				terms = [Tuple(parse.(Int, split(x,"⋅"))) for x in terms]

				B[i,j] = terms
			catch
				continue
			end

		end
	end

	return B
end

"""
    operatorString(MPOvec; op_String = nothing)

Returns all operator strings contained in an mpo chain in form of a list of integer tuples.
Example:

# Arguments:

"""
function stringMPO(MPOvec::Vector{MPO{localOp{DT,D}}}; op_String::Union{Array{OpString,2},Nothing} = nothing) where {DT<:Number, D}
	mpo = MPOvec[1]
	B = MPO_to_OpString(mpo)

	op_String = op_String == nothing ? B : stringM_mul(op_String,B)

	if size(MPOvec[2:end])[1] > 0
		stringMPO(MPOvec[2:end], op_String = op_String)
	else
		return convert_OpStringArrayToTuple(op_String)

	end
end

max_stars_repo_name: Gawatz/MPOmodule

max_issues_repo_name: Gawatz/MPOmodule

max_forks_repo_name: Gawatz/MPOmodule


text:
# MaxPlus Semiring Implementation
# MaxPlus is over real numbers with max as addition and + as multiplication
# The zero element is -inf and the one element is 0

struct MaxPlusSemiringElement{T<:Real} <: AbstractSemiringElement{T}
    val::T
end

# interface methods
val(x::MaxPlusSemiringElement) = x.val
+(l::MaxPlusSemiringElement, r::MaxPlusSemiringElement) = MaxPlusSemiringElement(max(l.val, r.val))
*(l::MaxPlusSemiringElement, r::MaxPlusSemiringElement) = MaxPlusSemiringElement(l.val+r.val)
zero(::Type{MaxPlusSemiringElement{T}}) where T = MaxPlusSemiringElement(-Inf)
one(::Type{MaxPlusSemiringElement{T}}) where T = MaxPlusSemiringElement{T}(0)

iscommutative(::Type{<:MaxPlusSemiringElement}) = true
isidempotent(::Type{<:MaxPlusSemiringElement}) = true
isstar(::Type{<:MaxPlusSemiringElement}) = true

# MinPlus Semiring Implementation
# MinPlus is over real numbers with min as addition and + as multiplication
# The zero element is +inf and the one element is 0

struct MinPlusSemiringElement{T<:Real} <: AbstractSemiringElement{T}
    val::T
end

# interface methods
val(x::MinPlusSemiringElement) = x.val
+(l::MinPlusSemiringElement, r::MinPlusSemiringElement) = MinPlusSemiringElement(min(l.val, r.val))
*(l::MinPlusSemiringElement, r::MinPlusSemiringElement) = MinPlusSemiringElement(l.val+r.val)
zero(::Type{MinPlusSemiringElement{T}}) where T = MinPlusSemiringElement(Inf)
one(::Type{MinPlusSemiringElement{T}}) where T = MinPlusSemiringElement{T}(0)

iscommutative(::Type{<:MinPlusSemiringElement}) = true
isidempotent(::Type{<:MinPlusSemiringElement}) = true
isstar(::Type{<:MinPlusSemiringElement}) = true

# Real (or Probabilistic) Semiring Implementation
# The Real Semiring is over real numbers with + as addition and * as multiplication
# The zero element is 0 and the one element is 1

struct RealSemiringElement{T<:Real} <: AbstractSemiringElement{T}
    val::T
end

# interface methods
val(x::RealSemiringElement) = x.val
+(l::RealSemiringElement, r::RealSemiringElement) = RealSemiringElement(l.val+r.val)
*(l::RealSemiringElement, r::RealSemiringElement) = RealSemiringElement(l.val*r.val)
zero(::Type{RealSemiringElement{T}}) where T = RealSemiringElement(0)
one(::Type{RealSemiringElement{T}}) where T = RealSemiringElement(1)

iscommutative(::Type{<:RealSemiringElement}) = true
isstar(::Type{<:RealSemiringElement}) = true

# optional methods
-(x::RealSemiringElement) = RealSemiringElement(-x.val)
/(l::RealSemiringElement, r::RealSemiringElement) = RealSemiringElement(l.val/r.val)
function star(x::RealSemiringElement)
    0 <= x.val < 1 ? RealSemiringElement(1/(1-x.val)) : throw(DomainError("$x must be within [0, 1)"))
end

for op in (:conj,) # required for matrix inversion?
    @eval begin
        $op(x::RealSemiringElement) = RealSemiringElement($op(x.val))
    end
end

# Boolean Semiring Implementation
# The Boolean Semiring is over {True, False} with OR as addition and AND as multiplication
# The zero element is False and the one element is True

struct BooleanSemiringElement{Bool} <: AbstractSemiringElement{Bool}
    val::Bool
end

# interface methods
val(x::BooleanSemiringElement) = x.val
+(l::BooleanSemiringElement, r::BooleanSemiringElement) = BooleanSemiringElement(l.val | r.val)
*(l::BooleanSemiringElement, r::BooleanSemiringElement) = BooleanSemiringElement(l.val & r.val)
zero(::Type{BooleanSemiringElement{T}}) where T = BooleanSemiringElement{T}(false)
one(::Type{BooleanSemiringElement{T}}) where T = BooleanSemiringElement{T}(true)

iscommutative(::Type{<:BooleanSemiringElement}) = true
isidempotent(::Type{<:BooleanSemiringElement}) = true
isstar(::Type{<:BooleanSemiringElement}) = true

#optional methods
star(x::R) where R<:BooleanSemiringElement = one(R)


text:
using LinearAlgebra
using StatsBase

function best_approx(A, k)
    U, _, _ = svd(A);
    return U[:, 1:k]
end

function ell(A, i; λ=0)
    ai = A[:, i]
    return ai' * pinv(A * A' + λ^2 * I) * ai
end

function ridge_weights(A, λ)
    _, n = size(A)
    weights = zeros(n)
    for i = 1:n
        weights[i] = ell(A, i, λ=λ)
    end
    return weights
end

function uniform_weights(A)
    _, n = size(A)
    return ones(n) / n
end

function squared_norms_weights(A)
    _, n = size(A)
    return dropdims(sum(A .* A, dims=1), dims=1)
end

function samplemat(A, weights, k)
    columns = sample(1:n, Weights(weights), k)
    return A[:, columns]
end

function ridge_sampling(A, k, λ)
    return samplemat(A, ridge_weights(A, λ), k)
end

function uniform_sampling(A, k)
    return samplemat(A, uniform_weights(A), k)
end

function squared_norms_sampling(A, k)
    return samplemat(A, squared_norms_weights(A), k)
end

function hilbert(n)
    H = zeros(n, n)
    for i = 1:n
        for j = 1:n
            H[i, j] = 1 / (i + j - 1)
        end
    end
    return H
end

function robustproj(A, tol=1e-12)
    U, s, _ = svd(A)
    k = sum(s .> tol)
    return U[:, 1:k] * U[:, 1:k]'
end

function approxerror(A, C)
    return norm(A - robustproj(C) * A)
end


text:
function geometry(Ny, Ly, method)

    if method=="cheb"
        Dy,y   = cheb(Ny);
        y      = y*(0.5*Ly);
        Dy     = (2/Ly)*Dy;
        Dy2    = Dy*Dy;
        return Dy, Dy2, y
    end

    #=
    if method=="FD2" # BROKEN
        y = range(-Ly/2, Ly/2, length=Ny+1);
        hy=y[2]-y[1];
        e = ones(Ny+1,1);
        Dy = Tridiagonal(-1*e[1:Ny]/(2*hy), 0*e[:], e[1:Ny]/(2*hy));
        Dy[1,1:2] = [-1 1]/hy;
        Dy[Ny+1,Ny:Ny+1] = [-1 1]/hy;
        Dy2 = Tridiagonal(e[1:Ny]/(hy^2), -2*e[:]/(hy^2), e[1:Ny]/(hy^2));
        #Dy2[1,1:3] = [1 -2 1]/hy^2;
        #Dy2[Ny+1,Ny-1:Ny+1] = [1 -2 1]/hy^2;
        return Dy,Dy2,y
    end
    =#

end


text:
using ..Utils

"Map from literal to LogicalΔNode"
const LitCache = Dict{Lit, LogicalΔNode}

"Use literal to represent constraint (1 to X, -1 to not X), 0 to represent true"
const ⊤ = convert(Lit, 0)

"""
Learning from data a structured-decomposable circuit with several structure learning algorithms
"""
function learn_struct_prob_circuit(data::Union{XData, WXData};
        pseudocount = 1.0, algo = "chow-liu", algo_kwargs=(α=1.0, clt_root="graph_center"), vtree = "chow-liu", vtree_kwargs=(vtree_mode="balanced",))
    if algo == "chow-liu"
        clt = learn_chow_liu_tree(data; algo_kwargs...)
        vtree = learn_vtree_from_clt(clt; vtree_kwargs...);
        pc = compile_psdd_from_clt(clt, vtree);
        estimate_parameters(pc, convert(XBatches,data); pseudocount = pseudocount)
        pc, vtree
    else
        error("Cannot learn a structured-decomposable circuit with algorithm $algo")
    end
end

#############
# Learn PlainVtree from CLT
#############

"
Learn a vtree from clt,
with strategy (close to) `linear` or `balanced`
"
function learn_vtree_from_clt(clt::CLT; vtree_mode::String)::PlainVtree
    roots = [i for (i, x) in enumerate(parent_vector(clt)) if x == 0]
    rootnode = construct_children(Var.(roots), clt, vtree_mode)

    return node2dag(rootnode)
end

function construct_node(v::Var, clt::CLT, strategy::String)::PlainVtreeNode
    children = Var.(outneighbors(clt, v))
    if isempty(children) # leaf node
        return PlainVtreeLeafNode(v)
    else
        right = construct_children(children, clt, strategy)
        return add_parent(v, right)
    end
end

function construct_children(children::Vector{Var}, clt::CLT, strategy::String)::PlainVtreeNode
    sorted_vars = sort(collect(children))
    children_nodes = Vector{PlainVtreeNode}()
    foreach(x -> push!(children_nodes, construct_node(x, clt, strategy)), sorted_vars)

    if strategy == "linear"
        construct_children_linear(children_nodes, clt)
    elseif strategy == "balanced"
        construct_children_balanced(children_nodes, clt)
    else
        throw("Unknown type of strategy")
    end
end

function construct_children_linear(children_nodes::Vector{PlainVtreeNode}, clt::CLT)::PlainVtreeNode
    children_nodes = Iterators.Stateful(reverse(children_nodes))

    right = popfirst!(children_nodes)
    for left in children_nodes
        right = PlainVtreeInnerNode(left, right)
    end
    return right
end

function construct_children_balanced(children_nodes::Vector{PlainVtreeNode}, clt::CLT)::PlainVtreeNode
    if length(children_nodes) == 1
        return children_nodes[1]
    elseif length(children_nodes) == 2
        return PlainVtreeInnerNode(children_nodes[1], children_nodes[2])
    else
        len = trunc(Int64, length(children_nodes) / 2)
        left = construct_children_balanced(children_nodes[1 : len], clt)
        right = construct_children_balanced(children_nodes[len + 1 : end], clt)
        return PlainVtreeInnerNode(left, right)
    end
end

function add_parent(parent::Var, children::PlainVtreeNode)
    return PlainVtreeInnerNode(PlainVtreeLeafNode(parent), children)
end

#####################
# Compile PSDD from CLT and vtree
#####################

"Compile a psdd circuit from clt and vtree"
function compile_psdd_from_clt(clt::MetaDiGraph, vtree::PlainVtree)
    order = node2dag(vtree[end])
    parent_clt = Var.(parent_vector(clt))

    lin = Vector{ProbΔNode}()
    prob_cache = ProbCache()
    lit_cache = LitCache()
    v2p = Dict{PlainVtreeNode, ProbΔ}()

    get_params(cpt::Dict) = length(cpt) == 2 ? [cpt[1], cpt[0]] : [cpt[(1,1)], cpt[(0,1)], cpt[(1,0)], cpt[(0,0)]]
    function add_mapping!(v::PlainVtreeNode, circuits::ProbΔ)
        if !haskey(v2p, v); v2p[v] = Vector{ProbΔNode}(); end
        foreach(c -> if !(c in v2p[v]) push!(v2p[v], c);end, circuits)
    end

    # compile vtree leaf node to terminal/true node
    function compile_from_vtree_node(v::PlainVtreeLeafNode)
        var = v.var
        children = Var.(outneighbors(clt, var))
        cpt = get_prop(clt, var, :cpt)
        parent = parent_clt[var]
        if isequal(children, [])
            circuit = compile_true_nodes(var, v, get_params(cpt), lit_cache, prob_cache, lin)
        else
            circuit = compile_literal_nodes(var, v, get_params(cpt), lit_cache, prob_cache, lin)
        end
        add_mapping!(v, circuit)
    end

    # compile to decision node
    function compile_from_vtree_node(v::PlainVtreeInnerNode)
        left_var = left_most_child(v.left).var
        right_var = left_most_child(v.right).var
        left_circuit = v2p[v.left]
        right_circuit = v2p[v.right]

        if parent_clt[left_var] == parent_clt[right_var] # two nodes are independent, compile to seperate decision nodes
            circuit = [compile_decision_node([l], [r], v, [1.0], prob_cache, lin) for (l, r) in zip(left_circuit, right_circuit)]
        elseif left_var == parent_clt[right_var] # conditioned on left
            cpt = get_prop(clt, left_var, :cpt)
            circuit = compile_decision_nodes(left_circuit, right_circuit, v, get_params(cpt), prob_cache, lin)
        else
            throw("PlainVtree are not learned from the same CLT")
        end
        add_mapping!(v, circuit)
    end

    foreach(compile_from_vtree_node, vtree)
    return lin
end

#####################
# Construct probabilistic circuit node
#####################

prob_children(n, prob_cache) =
    copy_with_eltype(map(c -> prob_cache[c], n.children), ProbΔNode{<:StructLogicalΔNode})

"Add leaf nodes to circuit `lin`"
function add_prob_leaf_node(var::Var, vtree::PlainVtreeLeafNode, lit_cache::LitCache, prob_cache::ProbCache, lin)
    pos = StructLiteralNode{PlainVtreeNode}( var2lit(var), vtree)
    neg = StructLiteralNode{PlainVtreeNode}(-var2lit(var), vtree)
    lit_cache[var2lit(var)] = pos
    lit_cache[-var2lit(var)] = neg
    pos2 = ProbLiteral(pos)
    neg2 = ProbLiteral(neg)
    prob_cache[pos] = pos2
    prob_cache[neg] = neg2
    push!(lin, pos2)
    push!(lin, neg2)
    return (pos2, neg2)
end

"Add prob⋀ node to circuit `lin`"
function add_prob⋀_node(children::ProbΔ, vtree::PlainVtreeInnerNode, prob_cache::ProbCache, lin)::Prob⋀
    logic = Struct⋀Node{PlainVtreeNode}([c.origin for c in children], vtree)
    prob = Prob⋀(logic, prob_children(logic, prob_cache))
    prob_cache[logic] = prob
    push!(lin, prob)
    return prob
end

"Add prob⋁ node to circuit `lin`"
function add_prob⋁_node(children::ProbΔ, vtree::PlainVtreeNode, thetas::Vector{Float64}, prob_cache::ProbCache, lin)::Prob⋁
    logic = Struct⋁Node{PlainVtreeNode}([c.origin for c in children], vtree)
    prob = Prob⋁(logic, prob_children(logic, prob_cache))
    prob.log_thetas = log.(thetas)
    prob_cache[logic] = prob
    push!(lin, prob)
    return prob
end

"Construct decision nodes given `primes` and `subs`"
function compile_decision_node(primes::ProbΔ, subs::ProbΔ, vtree::PlainVtreeInnerNode, params::Vector{Float64}, prob_cache::ProbCache, lin)
    elements = [add_prob⋀_node([prime, sub], vtree, prob_cache, lin) for (prime, sub) in zip(primes, subs)]
    return add_prob⋁_node(elements, vtree, params, prob_cache, lin)
end

"Construct literal nodes given variable `var`"
function compile_literal_nodes(var::Var, vtree::PlainVtreeLeafNode, probs::Vector{Float64}, lit_cache::LitCache, prob_cache::ProbCache, lin)
    (pos, neg) = add_prob_leaf_node(var, vtree, lit_cache, prob_cache, lin)
    return [pos, neg]
end

"Construct true nodes given variable `var`"
function compile_true_nodes(var::Var, vtree::PlainVtreeLeafNode, probs::Vector{Float64}, lit_cache::LitCache, prob_cache::ProbCache, lin)
    (pos, neg) = add_prob_leaf_node(var, vtree, lit_cache, prob_cache, lin)
    return [add_prob⋁_node([pos, neg], vtree, probs[i:i+1], prob_cache, lin) for i in 1:2:length(probs)]
end

"Construct decision nodes conditiond on different distribution"
function compile_decision_nodes(primes::ProbΔ, subs::ProbΔ, vtree::PlainVtreeInnerNode, params::Vector{Float64}, prob_cache::ProbCache, lin)
    return [compile_decision_node(primes, subs, vtree, params[i:i+1], prob_cache, lin) for i in 1:2:length(params)]
end

#####################
# Map and cache constraints
#####################

function set_base(index, n::StructLiteralNode, bases)
    if positive(n)
        bases[n][variable(n)] = 1
    else
        bases[n][variable(n)] = -1
    end
end

function set_base(index, n::Struct⋁Node, bases)
    len = num_children(n)
    temp = sum([bases[c] for c in n.children])
    bases[n] = map(x-> if x == len 1; elseif -x == len; -1; else 0; end, temp)
end

function set_base(index, n::Struct⋀Node, bases)
    bases[n] = sum([bases[c] for c in n.children])
end

function calculate_all_bases(circuit::ProbΔ)::BaseCache
    num_var = num_variables(circuit[end].origin.vtree)
    bases = BaseCache()
    foreach(n -> bases[n.origin] = fill(⊤, num_var), circuit)
    foreach(n -> set_base(n[1], n[2].origin, bases), enumerate(circuit))
    @assert all(bases[circuit[end].origin] .== ⊤) "Base of root node should be true"
    return bases
end

#####################
# Compile fully factorized PSDD from vtree, all variables are independent initially
#####################

function compile_fully_factorized_psdd_from_vtree(vtree::PlainVtree)::ProbΔ

    function ful_factor_node(v::PlainVtreeLeafNode, lit_cache::LitCache, prob_cache::ProbCache, v2n, lin)
        var = variables(v)[1]
        pos, neg = add_prob_leaf_node(var, v, lit_cache, prob_cache, lin)
        prob_or = add_prob⋁_node([pos, neg], v, [0.5, 0.5], prob_cache, lin)
        v2n[v] = prob_or
        nothing
    end

    function ful_factor_node(v::PlainVtreeInnerNode, lit_cache::LitCache, prob_cache::ProbCache, v2n, lin)
        left = v2n[v.left]
        right = v2n[v.right]
        prob_and = add_prob⋀_node([left, right], v, prob_cache, lin)
        prob_or = add_prob⋁_node([prob_and], v, [1.0], prob_cache, lin)
        v2n[v] = prob_or
        nothing
    end

    lin = Vector{ProbΔNode}()
    prob_cache = ProbCache()
    lit_cache = LitCache()
    v2n = Dict{PlainVtreeNode, ProbΔNode}()

    for v in vtree
        ful_factor_node(v, lit_cache, prob_cache, v2n, lin)
    end

    lin
end
"""
Learning from data a structured-decomposable credal circuit with several structure learning algorithms
"""
function learn_struct_credal_circuit(data::Union{XData, WXData}, s_idm::Float64;
        pseudocount = 1.0, algo = "chow-liu", algo_kwargs=(α=1.0, clt_root="graph_center"), vtree = "chow-liu", vtree_kwargs=(vtree_mode="balanced",))
    if algo == "chow-liu"
        clt = learn_chow_liu_tree(data; algo_kwargs...)
        vtree = learn_vtree_from_clt(clt; vtree_kwargs...);
        pc = compile_csdd_from_clt(clt, vtree);
        estimate_credal_parameters(pc, convert(XBatches,data),s_idm ; pseudocount = pseudocount)
        pc, vtree
    else
        error("Cannot learn a structured-decomposable circuit with algorithm $algo")
    end
end

#####################
# Compile CSDD from CLT and vtree
#####################

"Compile a psdd circuit from clt and vtree"
function compile_csdd_from_clt(clt::MetaDiGraph, vtree::PlainVtree)
    order = node2dag(vtree[end])
    parent_clt = Var.(parent_vector(clt))

    lin = Vector{CredalΔNode}()
    cred_cache = CredalCache()
    lit_cache = LitCache()
    v2p = Dict{PlainVtreeNode, CredalΔ}()

    get_params_cred(cpt::Dict) = length(cpt) == 2 ? [cpt[1], cpt[0]] : [cpt[(1,1)], cpt[(0,1)], cpt[(1,0)], cpt[(0,0)]]
    function add_mapping_cred!(v::PlainVtreeNode, circuits::CredalΔ)
        if !haskey(v2p, v); v2p[v] = Vector{CredalΔNode}(); end
        foreach(c -> if !(c in v2p[v]) push!(v2p[v], c);end, circuits)
    end

    # compile vtree leaf node to terminal/true node
    function compile_from_vtree_node_cred(v::PlainVtreeLeafNode)
        var = v.var
        children = Var.(outneighbors(clt, var))
        cpt = get_prop(clt, var, :cpt)
        parent = parent_clt[var]
        if isequal(children, [])
            circuit = compile_true_nodes_cred(var, v, get_params_cred(cpt), lit_cache, cred_cache, lin)
        else
            circuit = compile_literal_nodes_cred(var, v, get_params_cred(cpt), lit_cache, cred_cache, lin)
        end
        add_mapping_cred!(v, circuit)
    end

    # compile to decision node
    function compile_from_vtree_node_cred(v::PlainVtreeInnerNode)
        left_var = left_most_child(v.left).var
        right_var = left_most_child(v.right).var
        left_circuit = v2p[v.left]
        right_circuit = v2p[v.right]

        if parent_clt[left_var] == parent_clt[right_var] # two nodes are independent, compile to seperate decision nodes
            circuit = [compile_decision_node_cred([l], [r], v, [1.0], cred_cache, lin) for (l, r) in zip(left_circuit, right_circuit)]
        elseif left_var == parent_clt[right_var] # conditioned on left
            cpt = get_prop(clt, left_var, :cpt)
            circuit = compile_decision_nodes_cred(left_circuit, right_circuit, v, get_params_cred(cpt), cred_cache, lin)
        else
            throw("PlainVtree are not learned from the same CLT")
        end
        add_mapping_cred!(v, circuit)
    end

    foreach(compile_from_vtree_node_cred, vtree)
    return lin
end

#####################
# Construct credal circuit node
#####################

cred_children(n, cred_cache) =
    copy_with_eltype(map(c -> cred_cache[c], n.children), CredalΔNode{<:StructLogicalΔNode})

"Add leaf nodes to circuit `lin`"
function add_cred_leaf_node(var::Var, vtree::PlainVtreeLeafNode, lit_cache::LitCache, cred_cache::CredalCache, lin)
    pos = StructLiteralNode{PlainVtreeNode}( var2lit(var), vtree)
    neg = StructLiteralNode{PlainVtreeNode}(-var2lit(var), vtree)
    lit_cache[var2lit(var)] = pos
    lit_cache[-var2lit(var)] = neg
    pos2 = CredalLiteral(pos)
    neg2 = CredalLiteral(neg)
    cred_cache[pos] = pos2
    cred_cache[neg] = neg2
    push!(lin, pos2)
    push!(lin, neg2)
    return (pos2, neg2)
end

"Add cred⋀ node to circuit `lin`"
function add_cred⋀_node(children::CredalΔ, vtree::PlainVtreeInnerNode, cred_cache::CredalCache, lin)::Credal⋀
    logic = Struct⋀Node{PlainVtreeNode}([c.origin for c in children], vtree)
    cred = Credal⋀(logic, cred_children(logic, cred_cache))
    cred_cache[logic] = cred
    push!(lin, cred)
    return cred
end

"Add cred⋁ node to circuit `lin`"
function add_cred⋁_node(children::CredalΔ, vtree::PlainVtreeNode, thetas::Vector{Float64}, cred_cache::CredalCache, lin)::Credal⋁
    logic = Struct⋁Node{PlainVtreeNode}([c.origin for c in children], vtree)
    cred = Credal⋁(logic, cred_children(logic, cred_cache))
    cred.log_thetas = log.(thetas)
    cred_cache[logic] = cred
    push!(lin, cred)
    return cred
end

"Construct decision nodes given `primes` and `subs`"
function compile_decision_node_cred(primes::CredalΔ, subs::CredalΔ, vtree::PlainVtreeInnerNode, params::Vector{Float64}, cred_cache::CredalCache, lin)
    elements = [add_cred⋀_node([prime, sub], vtree, cred_cache, lin) for (prime, sub) in zip(primes, subs)]
    return add_cred⋁_node(elements, vtree, params, cred_cache, lin)
end

"Construct literal nodes given variable `var`"
function compile_literal_nodes_cred(var::Var, vtree::PlainVtreeLeafNode, probs::Vector{Float64}, lit_cache::LitCache, cred_cache::CredalCache, lin)
    (pos, neg) = add_cred_leaf_node(var, vtree, lit_cache, cred_cache, lin)
    return [pos, neg]
end

"Construct true nodes given variable `var`"
function compile_true_nodes_cred(var::Var, vtree::PlainVtreeLeafNode, probs::Vector{Float64}, lit_cache::LitCache, cred_cache::CredalCache, lin)
    (pos, neg) = add_cred_leaf_node(var, vtree, lit_cache, cred_cache, lin)
    return [add_cred⋁_node([pos, neg], vtree, probs[i:i+1], cred_cache, lin) for i in 1:2:length(probs)]
end

"Construct decision nodes conditiond on different distribution"
function compile_decision_nodes_cred(primes::CredalΔ, subs::CredalΔ, vtree::PlainVtreeInnerNode, params::Vector{Float64}, cred_cache::CredalCache, lin)
    return [compile_decision_node_cred(primes, subs, vtree, params[i:i+1], cred_cache, lin) for i in 1:2:length(params)]
end

max_stars_repo_name: IDSIA/cjuice

max_issues_repo_name: IDSIA/cjuice

max_forks_repo_name: IDSIA/cjuice


text:
include("common.jl")

include("test.jl")

say("bufan")

max_stars_repo_name: innerlee/testRand

max_issues_repo_name: innerlee/testRand

max_forks_repo_name: innerlee/testRand


text:
using Iris
using Test

@testset "Curlcurls" begin
    @testset "1D Symmetric" begin

    end
end


text:
using FinancialMonteCarlo, Test
@show "Test Structs Building"
S0 = 100.0;
K = 100.0;
Kneg = -100.0;
r = 0.02;
T = 1.0;
Tneg = -1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
theta1 = 0.01;
k1 = 0.03;
sigma1 = 0.02;
underlying_name = "ENJ"
mc = MonteCarloConfiguration(Nsim, Nstep, 1);
mc2 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.SerialMode(), 1);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());

toll = 0.8;
Model_enj = BlackScholesProcess(sigma, Underlying(S0, d));
Model_abpl = BlackScholesProcess(sigma, Underlying(S0, d));
FinancialMonteCarlo.Curve([0.0, 0.1], [0.0, 0.1]);
Model = GaussianCopulaNVariateProcess(Model_enj, Model_abpl, 0.0)
param_ = FinancialMonteCarlo.get_parameters(Model_abpl);

@test_throws(ErrorException, FinancialMonteCarlo.set_parameters!(Model_abpl, [1, 2, 3, 4]))

port_ = "eni_aap" → Model;
rfCurve = ZeroRate(r);
rfCurve2 = ZeroRate([0.00, 0.02], T);
@test_throws(ErrorException, port_ + port_)
@test_throws(ErrorException, underlying_name → Model)
@test_throws(ErrorException, Underlying(-S0, rfCurve2.r))
@test_throws(ErrorException, "c_i" → Forward(1.0))
@test_throws(ErrorException, Underlying(-S0))
@test_throws(ErrorException, EuropeanOptionND(-1.0, 1.0))
@test_throws(ErrorException, EuropeanOptionND(1.0, -1.0))
@test_throws(ErrorException, MonteCarloConfiguration(-Nsim, Nstep))
@test_throws(ErrorException, MonteCarloConfiguration(Nsim + 1, Nstep, FinancialMonteCarlo.AntitheticMC()));
@test_throws(ErrorException, MonteCarloConfiguration(Nsim, -Nstep))
#####################################################################
#Payoff Structs Test: Negative TTM
@test_throws(ErrorException, Forward(Tneg));
@test_throws(ErrorException, EuropeanOption(Tneg, K));
@test_throws(ErrorException, AmericanOption(Tneg, K));
@test_throws(ErrorException, BermudanOption([Tneg], K));
@test_throws(ErrorException, BinaryAmericanOption(Tneg, K));
@test_throws(ErrorException, BinaryEuropeanOption(Tneg, K));
@test_throws(ErrorException, BarrierOptionDownOut(Tneg, K, D));
@test_throws(ErrorException, BarrierOptionUpOut(Tneg, K, D));
@test_throws(ErrorException, BarrierOptionUpIn(Tneg, K, D));
@test_throws(ErrorException, BarrierOptionDownIn(Tneg, K, D));
@test_throws(ErrorException, AsianFloatingStrikeOption(Tneg));
@test_throws(ErrorException, AsianFixedStrikeOption(Tneg, K));
@test_throws(ErrorException, DoubleBarrierOption(Tneg, K, K * 10, D));

# Negative Strike and Barriers
@test_throws(ErrorException, EuropeanOption(T, Kneg));
@test_throws(ErrorException, AmericanOption(T, Kneg));
@test_throws(ErrorException, BermudanOption([T], Kneg));
@test_throws(ErrorException, BinaryEuropeanOption(T, Kneg));
@test_throws(ErrorException, BinaryAmericanOption(T, Kneg));
@test_throws(ErrorException, AsianFixedStrikeOption(T, Kneg));
@test_throws(ErrorException, BarrierOptionDownOut(T, Kneg, D));
@test_throws(ErrorException, BarrierOptionDownOut(T, K, Kneg));
@test_throws(ErrorException, BarrierOptionDownIn(T, Kneg, D));
@test_throws(ErrorException, BarrierOptionDownIn(T, K, Kneg));
@test_throws(ErrorException, BarrierOptionUpIn(T, Kneg, D));
@test_throws(ErrorException, BarrierOptionUpIn(T, K, Kneg));
@test_throws(ErrorException, BarrierOptionUpOut(T, Kneg, D));
@test_throws(ErrorException, BarrierOptionUpOut(T, K, Kneg));
@test_throws(ErrorException, DoubleBarrierOption(T, K, K, Kneg));
@test_throws(ErrorException, DoubleBarrierOption(T, Kneg, K, K));
@test_throws(ErrorException, DoubleBarrierOption(T, K, Kneg, K));

@test FinancialMonteCarlo.maturity(Spot()) == 0.0
@test 2.0 * EuropeanOption(T, K) == EuropeanOption(T, K) + 1.0 * EuropeanOption(T, K)
@test 2.0 * EuropeanOption(T, K) == EuropeanOption(T, K) + EuropeanOption(T, K)

########

import FinancialMonteCarlo.AbstractPayoff

struct tmptype2 <: FinancialMonteCarlo.AbstractPayoff{Float64}
    T::Float64
end

@test_throws(ErrorException, payoff([1.0 1.0; 1.0 1.0], tmptype2(1.0), rfCurve));


text:
using JuMP
using AmplNLWriter, Bonmin_jll, Couenne_jll
using Ipopt
using Plots
using MathOptInterface

const MOI = MathOptInterface

include("types.jl")
include("utils.jl")
include("input.jl")
include("model.jl")
include("enumeration.jl")
include("output.jl")

const density = 1000
const g = 9.81
const deepest = 200.0

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
# Copyright 2018 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md

using ModuleInterfaceTools

@api test

@api extend InternedStrings

@api list InternedStrings

@api def testcase begin
    myname = "Scott Paul Jones"
end

@testset "@api def <name> <expr>" begin
    @testcase
    @test myname == "Scott Paul Jones"
end

intern(x::Integer) = 1
intern(x::Float64) = 2

@testset "Function extension" begin
    @test typeof(intern("foo")) == String
    @test intern(1) == 1
    @test intern(2.0) == 2
    @test intern("th") == "th"
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module GslibIO

using FileIO
using Printf
using DelimitedFiles

using GeoStatsBase

"""
    load(file)

Load grid properties from `file`.
"""
function load(file::File{format"GSLIB"})
  open(file) do f
    fs = stream(f)

    # skip header
    skipchars(_ -> false, fs, linecomment='#')

    # read dimensions
    nx, ny, nz = parse.(Int,     split(readline(fs)))
    ox, oy, oz = parse.(Float64, split(readline(fs)))
    sx, sy, sz = parse.(Float64, split(readline(fs)))

    # read property names
    vars = Symbol.(split(readline(fs)))

    # read property values
    X = readdlm(fs)

    # create data dictionary
    data = Dict(vars[j] => reshape(X[:,j], nx, ny, nz) for j in 1:size(X,2))

    RegularGridData(data, (ox,oy,oz), (sx,sy,sz))
  end
end

"""
    save(file, properties, propsize; [optional parameters])

Save 1D `properties`, which originally had 3D size `propsize`.
"""
function save(file::File{format"GSLIB"},
              properties::Vector{V}, propsize::Tuple;
              origin=(0.,0.,0.), spacing=(1.,1.,1.),
              header="", propnames="") where {T<:Real,V<:AbstractArray{T,1}}
  # default property names
  isempty(propnames) && (propnames = ["prop$i" for i=1:length(properties)])
  @assert length(propnames) == length(properties) "number of property names must match number of properties"

  # convert vector of names to a long string
  propnames = join(propnames, " ")

  # collect all properties in a big matrix
  P = hcat(properties...)

  open(file, "w") do f
    # write header
    write(f, "# This file was generated with GslibIO.jl\n")
    !isempty(header) && write(f, "#\n# "*header*"\n")

    # write dimensions
    write(f, @sprintf("%i %i %i\n", propsize...))
    write(f, @sprintf("%f %f %f\n", origin...))
    write(f, @sprintf("%f %f %f\n", spacing...))

    # write property name and values
    write(f, "$propnames\n")
    writedlm(stream(f), P, ' ')
  end
end

"""
    save(file, properties, [optional parameters])

Save 3D `properties` by first flattening them into 1D properties.
"""
function save(file::File{format"GSLIB"}, properties::Vector{A};
              kwargs...) where {T<:Real,A<:AbstractArray{T,3}}
  # sanity checks
  @assert length(Set(size.(properties))) == 1 "properties must have the same size"

  # retrieve grid size
  propsize = size(properties[1])

  # flatten and proceed with pipeline
  flatprops = [prop[:] for prop in properties]

  save(file, flatprops, propsize, kwargs...)
end

"""
    save(file, property)

Save single 3D `property` by wrapping it into a singleton collection.
"""
function save(file::File{format"GSLIB"},
              property::A; kwargs...) where {T<:Real,A<:AbstractArray{T,3}}
  save(file, [property]; kwargs...)
end

"""
    save(file, grid)

Save `grid` of type `RegularGridData` to file.
"""
function save(file::File{format"GSLIB"}, grid::RegularGridData{<:Any,3})
  propnames = sort(collect(keys(variables(grid))))
  properties = [vec(grid[propname]) for propname in propnames]
  save(file, properties, size(grid), origin=origin(grid),
       spacing=spacing(grid), propnames=propnames)
end

end


text:
"""
    ORBIT{N, VT, AM} <: AbstractContinuousPost

Implementation of discrete-time integration for deterministic linear ODEs
for singleton initial conditions and proper input sets.

## Fields

- `δ`            -- step-size of the discretization
- `approx_model` -- (optional, default: `NoBloating`) approximation model

The type fields are:

- `N`   -- number type of the step-size
- `AM`  -- type of the approximation model

## Notes

If the input set is a singleton, the result is the discrete-time sequence of
states that correspond to the analytic solution.
If the input set is a dense set, we sample a new input signal in every step.

Reachability solutions computed with the `ORBIT` algorithm have its own plot
recipe, producing a scatter plot. Use the `markershape` option to change the shape
of the markers used for the scatter plot. To "connect" the points, use the option
`seriestype=:path`. For additional options see the [Plots.jl documentation](https://docs.juliaplots.org/latest/generated/attributes_series/).
"""
struct ORBIT{N, VT, AM} <: AbstractContinuousPost
    δ::N
    approx_model::AM
end

step_size(alg::ORBIT) = alg.δ
numtype(::ORBIT{N}) where {N} = N

# convenience constructor using symbols
function ORBIT(; δ::N, approx_model::AM=NoBloating(exp=:base, setops=:concrete)) where {N, AM}
    VT = Vector{N}
    return ORBIT{N, VT, AM}(δ, approx_model)
end

function setrep(::ORBIT{N, VT, AM}) where {N, VT, AM}
    Singleton{N, VT}
end

function rsetrep(alg::ORBIT{N}) where {N}
    ST = sertrep(alg)
    ReachSet{N, ST}
end

include("post.jl")


text:
Base.:*(times::Integer, w::Wire) = Wire(repeat(w.values, outer=times))

function Base.:+{R,S}(lhs::Wire{R}, rhs::Wire{S})
  (length(R) == length(S)) || throw(SizeMismatchError())

  #for now, disallow bitvectors of length more than 64
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end

  #create a new uninitialized wire of the same length
  result = Wire{range(length(R))}(BitVector(length(R)), trues(length(R)))
  result.values.chunks[1] = lhs.values.chunks[1] + rhs.values.chunks[1]
  result
end

function Base.:-{R,S}(lhs::Wire{R}, rhs::Wire{S})
  (length(R) == length(S)) || throw(SizeMismatchError())

  #for now, disallow bitvectors of length more than 64
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end

  #create a new uninitialized wire of the same length
  result = Wire{range(length(R))}(BitVector(length(R)), trues(length(R)))
  result.values.chunks[1] = lhs.values.chunks[1] - rhs.values.chunks[1]
  result
end

function Base.:*{R,S}(lhs::Wire{R}, rhs::Wire{S})
  new_length = length(R) + length(S)

  if (length(R) + length(S) > 64)
    warn("current > 64 bit results not supported.")
    throw(SizeMismatchError())
  end

  result = Wire{range(new_length)}(BitVector(new_length), trues(new_length))
  result.values.chunks[1] = lhs.values.chunks[1] * rhs.values.chunks[1]
  result
end

function Base.:-{R}(tgt::Wire{R})
  #for now, disallow bitvectors of length more than 64
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end

  #create a new uninitialized wire of the same length
  result = Wire{range(length(R))}(BitVector(length(R)), trues(length(R)))
  result.values.chunks[1] = -tgt.values.chunks[1]
  result
end

function trimInt64(n::Integer)
  ~(-(one(UInt64) << n))
end

function Base.:<{R,S}(lhs::Wire{R}, rhs::Wire{S})
  (length(R) == length(S)) || throw(SizeMismatchError())
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end
  assigned(lhs) && assigned(rhs) || throw(UnassignedError())

  mask = trimInt64(length(R))

  result = Wire((lhs.values.chunks[1] & mask) < (rhs.values.chunks[1] & mask))
end

# >, <=, >= come for free.

function Base.:(==){R,S}(lhs::Wire{R}, rhs::Wire{S})
  (length(R) == length(S)) || throw(SizeMismatchError())
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end
  assigned(lhs) && assigned(rhs) || throw(UnassignedError())

  mask = trimInt64(length(R))

  result = Wire((lhs.values.chunks[1] & mask) == (rhs.values.chunks[1] & mask))
end

################################################################################
## shifters

function Base.:(<<){R,S}(lhs::Wire{R}, rhs::Wire{S})
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end
  if (length(S) > 64)
    warn("> 64 bit wires for shifter RHS not supported.")
    throw(SizeMismatchError())
  end
  assigned(lhs) && assigned(rhs) || throw(UnassignedError())
  result = Wire(lhs.values.chunks[1] << rhs.values.chunks[1], length(R))
end

function Base.:(>>){R,S}(lhs::Wire{R}, rhs::Wire{S})
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end
  if (length(S) > 64)
    warn("> 64 bit wires for shifter RHS not supported.")
    throw(SizeMismatchError())
  end
  assigned(lhs) && assigned(rhs) || throw(UnassignedError())
  result = Wire(lhs.values.chunks[1] >> rhs.values.chunks[1], length(R))
end

function Base.:(>>>){R,S}(lhs::Wire{R}, rhs::Wire{S})
  if (length(R) > 64)
    warn("currently > 64 bit wires not supported.")
    throw(SizeMismatchError())
  end
  if (length(S) > 64)
    warn("> 64 bit wires for shifter RHS not supported.")
    throw(SizeMismatchError())
  end
  assigned(lhs) && assigned(rhs) || throw(UnassignedError())

  #if the top value is true, then we'll have to spoof negativity.
  #the trimInt64() function can help, by masking those invisible values
  #as ones.
  if (lhs.values[end] == true)
    lhs.values.chunks[1] |= ~trimInt64(length(R))
  end

  result = Wire(lhs.values.chunks[1] >> rhs.values.chunks[1], length(R))
end


text:
export
    AbstractTransform,
    transform,
    truncate_modes,
    inverse

"""
    AbstractTransform

## Interface

* `Base.ndims(<:AbstractTransform)`: N dims of modes
* `transform(<:AbstractTransform, 𝐱::AbstractArray)`: Apply the transform to 𝐱
* `truncate_modes(<:AbstractTransform, 𝐱_transformed::AbstractArray)`: Truncate modes that contribute to the noise
* `inverse(<:AbstractTransform, 𝐱_transformed::AbstractArray)`: Apply the inverse transform to 𝐱_transformed
"""
abstract type AbstractTransform end

include("fourier_transform.jl")


text:
# first pass: read whole file, save in binary format, count IDs for
# collation, store AM spell strings in the order of appearance

using AMDB:
    ColSpec, make_firstpass, data_path, firstpass_process_file, all_data_files,
    get_positions, AutoIndex
# using ArgCheck
using ByteParsers:
    Line, DateYYYYMMDD, PosInteger, ViewBytes, parsedtype, ByteVector, FixEmpty
# using DiscreteRanges
using DocStringExtensions: SIGNATURES
using JLD2
using FileIO
using LargeColumns: meta_path

# AMDB.preview_column("ALTER"; N = 10)

colspecs = [
    # 1, person id
    ColSpec("PENR", PosInteger(Int32)),
    # 2:3, spell start, Date
    ColSpec("STARTEND", AMDB.DatePair()),
    # 4, firm id, Int64 for some reason?
    ColSpec("BENR", FixEmpty(-1, PosInteger(Int64)); index_type = Int32), # missing as -1
    # 5, STICHTAG
    ColSpec("STICHTAG", AMDB.ParseDate()),
    # 6, labor market status
    ColSpec("AM", ViewBytes(), index_type = Int8),
    # 11, gender
    ColSpec("GESCHLECHT", ViewBytes(), index_type = Int8),
    # 12, alter
    ColSpec("ALTER", PosInteger(Int8)),
    # 17, number of employees,  seems to be capped at 1000?
    ColSpec("SUM_MA", FixEmpty(-1, PosInteger(Int16))), # missing as -1
    # 19, industry code, 4 digits, has strings like XXXX?
    ColSpec("NACE", ViewBytes(); index_type = Int16),
    # 21, geographical location; 3 digits
    ColSpec("RGS", ViewBytes(); index_type = UInt8), # enough to contain it, has around 130
    # 35, wage data
    ColSpec("AVG_BMG", FixEmpty(-1, PosInteger(Int32))), # missing wage as -1
]

dir = data_path("first_pass")
mkpath(dir)
fp = make_firstpass(dir, colspecs)

# quick check
@assert fp.colnames ==
    [:PENR, :STARTEND, :BENR, :STICHTAG, :AM, :GESCHLECHT, :ALTER, :SUM_MA,
     :NACE, :RGS, :AVG_BMG]

error_io = open(data_path("first_pass_errors.txt"), "w")
println(error_io,
        "first pass started at $(Dates.now()) on machine $(gethostname())")
for filename in all_data_files()
    println("processing $filename")
    err = firstpass_process_file(filename, fp)
    show(error_io, err)
end
close(error_io)
close(fp)

"""
    $SIGNATURES

Convert keys to string.
"""
collect_keys(accumulator::AutoIndex{<: Integer}) = collect(keys(accumulator))

function collect_keys(accumulator::AutoIndex{<:SubArray{UInt8}})
    map(String ∘ copy, keys(accumulator))
end

# save the keys
save(meta_path(dir, "meta.jld2"),
     # "id_counter", fp.orderedcounter,
     AMDB.META_INDEXED_KEYS, map(collect_keys, fp.accumulators),
     AMDB.META_INDEXED_POSITIONS, collect(get_positions(fp.multisubs)),
     AMDB.META_COLUMN_NAMES, fp.colnames)


text:
## test the main entry point in Celeste: the `infer` function
import JLD

"""
test infer with a single (run, camcol, field).
This is basically just to make sure it runs at all.
"""
function test_infer_single()
    # very small patch of sky that turns out to have 4 sources.
    # We checked that this patch is in the given field.
    box = ParallelRun.BoundingBox(164.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]
    result = ParallelRun.one_node_infer(field_triplets, datadir; box=box)
end

function load_surrounding_rcfs()
    wd = pwd()
    cd(datadir)

    # these rcfs are all the rcfs that overlap with (3900,6,269)
    run(`make RUN=3900 CAMCOL=6 FIELD=268`)
    run(`make RUN=3900 CAMCOL=6 FIELD=269`)
    run(`make RUN=3900 CAMCOL=6 FIELD=270`)
    run(`make RUN=4469 CAMCOL=5 FIELD=342`)
    run(`make RUN=4469 CAMCOL=5 FIELD=343`)
    run(`make RUN=4469 CAMCOL=6 FIELD=342`)
    run(`make RUN=4469 CAMCOL=6 FIELD=343`)
    run(`make RUN=4469 CAMCOL=6 FIELD=344`)

    cd(wd)
end

function test_infer_rcf()
    load_surrounding_rcfs()

    resfile = joinpath(datadir, "celeste-003900-6-0269.jld")
    rm(resfile, force=true)

    rcf = RunCamcolField(3900, 6, 269)
    objid = "1237662226208063492"
    ParallelRun.infer_rcf(rcf, datadir, datadir; objid=objid)

    @test isfile(resfile)
    println(filesize(resfile))
    @test filesize(resfile) > 1000  # should be about 15 KB
    rm(resfile)
end

function test_source_division_parallelism()
    box = ParallelRun.BoundingBox(164.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]
    results = ParallelRun.divide_sources_and_infer(box, datadir)
    @test length(results) == 4
end

function test_load_active_pixels()
    images, ea, one_body = gen_sample_star_dataset()

    # these images have 20 * 23 * 5 = 2300 pixels in total.
    # the star is bright but it doesn't cover the whole image.
    # it's hard to say exactly how many pixels should be active,
    # but not all of them, and not none of them.
    Infer.load_active_pixels!(ea.images, ea.patches; min_radius_pix=0.0)

    # most star light (>90%) should be recorded by the active pixels
    num_active_photons = 0.0
    num_active_pixels = 0
    total_pixels = 0
    for n in 1:ea.N
        img = ea.images[n]
        p = ea.patches[1, n]
        H2, W2 = size(p.active_pixel_bitmap)
        total_pixels += length(img.pixels)
        for w2 in 1:W2, h2 in 1:H2
            # (h2, w2) index the local patch, while (h, w) index the image
            h = p.bitmap_offset[1] + h2
            w = p.bitmap_offset[2] + w2
            num_active_photons += img.pixels[h, w] - img.epsilon_mat[h, w]
            num_active_pixels += 1
        end
    end

    @test 100 < num_active_pixels < total_pixels

    total_photons = 0.0
    for img in ea.images
        total_photons += sum(img.pixels) - sum(img.epsilon_mat)
    end

    @test num_active_photons <= total_photons  # sanity check
    @test num_active_photons > 0.9 * total_photons

    # super dim images
    for img in images
        img.pixels[:,:] = img.epsilon_mat[:,:]
    end

    # only 2 pixels per image are within 0.6 pixels of the
    # source's center (10.9, 11.5)
    Infer.load_active_pixels!(ea.images, ea.patches; min_radius_pix=0.6)

    for n in 1:ea.N
#  FIXME: is load active pixels off by (0.5, 0.5)?
#        @test sum(ea.patches[1,n].active_pixel_bitmap) == 2
    end
end

if test_long_running
    test_infer_rcf()
end

test_load_active_pixels()
test_source_division_parallelism()
test_infer_single()


text:
using UpdateIndex
using Base.Test

A = zeros(Int, 5,5)
@update! A[3,2] += 2
@test A[3,2] == 2
updateindex!(A, -, 2, 3, 2)
@test A[3,2] == 0

Asp = spzeros(5,5)
Asp[3,3] += 1.0
@update! Asp[3,3] += 2.0
@test Asp[3,3] == 3.0
updateindex!(Asp, -, 2.0, 3, 3)
@test Asp[3,3] == 1.0

met = @which(updateindex!(Asp, -, 2.0, 3, 3))
# Check overloading is woking
if VERSION <= v"0.5-dev"
    @test string(met.sig.parameters[1].name) == "SparseMatrixCSC"
else
    @test string(met.sig.parameters[2].name) == "SparseMatrixCSC"
end