text:
KL = include("./KnetLayers/src/KnetLayers.jl")
onnx_utils = include("./onnx_utils.jl")

function ONNXtoGraph(file)
    """
    Given the path of an ONNX file, construct its graph.
    """
    f = readproto(open(file), Proto.ModelProto());
    convert(f).graph
end

function PrintGraph(g)
    """
    Prints the ONNX graph.
    """
    ops = [] # operator list
    println("model inputs: ", (x->x.name).(g.input))
    println("model outputs: ", (x->x.name).(g.output))
    for (i, node) in enumerate(g.node)
        if join(node.op_type) ∉ ops
            push!(ops, join(node.op_type))
        end
        print("(op", i, ") ", node.op_type)
        println()
        for (i, input) in enumerate(node.input)
            println("\tinput", i, ": " , input)
        end
        for (i, output) in enumerate(node.output)
            println("\toutput", i, ": " , output)
        end
    end
    println("Operator list: ", ops)
end

function convert(node, g)
    """
    Given a node, calls the appropriate constructor for the corresponding (args, layer, outs)
    """
    if node.op_type == "Add"; return converter_add(node, g); end
    if node.op_type == "AveragePool"; return converter_avgpool(node,g); end
    if node.op_type == "BatchNormalization"; return converter_batchnorm(node, g); end
    if node.op_type == "Cast"; return converter_cast(node,g); end
    if node.op_type == "Concat"; return converter_concat(node,g); end
    if node.op_type == "Constant"; return converter_constant(node, g); end
    if node.op_type == "ConstantOfShape"; return converter_constantOfShape(node, g); end
    if node.op_type == "Conv"; return converter_cnn(node,g); end
    if node.op_type == "Div"; return converter_div(node,g); end
    if node.op_type == "Dropout"; return converter_dropout(node,g); end
    if node.op_type == "Flatten"; return converter_flatten(node,g); end
    if node.op_type == "Gather"; return converter_gather(node,g); end
    if node.op_type == "Gemm"; return converter_gemm(node, g); end
    if node.op_type == "GlobalAveragePool"; return converter_globalAveragePool(node, g); end
    if node.op_type == "Identity"; return converter_identity(node, g); end
    if node.op_type == "LeakyRelu"; return converter_leakyrelu(node,g); end
    if node.op_type == "MatMul"; return converter_matmul(node,g); end
    if node.op_type == "MaxPool"; return converter_maxpool(node,g); end
    if node.op_type == "Mul"; return converter_mul(node,g); end
    if node.op_type == "Relu"; return converter_relu(node, g); end
    if node.op_type == "Reshape"; return converter_reshape(node, g); end
    if node.op_type == "RNN"; return converter_rnn(node, g);
    if node.op_type == "Shape"; return converter_shape(node, g); end
    if node.op_type == "Softmax"; return converter_softmax(node, g); end
    if node.op_type == "Slice"; return converter_slice(node, g); end
    if node.op_type == "Sub"; return converter_sub(node,g); end
    if node.op_type == "Squeeze"; return converter_squeeze(node, g); end
    if node.op_type == "Transpose"; return converter_transpose(node,g); end
    if node.op_type == "Unsqueeze"; return converter_unsqueeze(node,g); end
    else; println("ONNX Operation not yet implemented: ", node.op_type);
    end
end

# Unit Layer for Testing
struct dummy_layer
end
(d::dummy_layer)(x) = x

function converter_rnn(node, g)
    return (node.input, dummy_layer(), node.output)
end

# TODO: convert all onnx params to knet params (use convert_params(w))
"""
Converters Begin Here
# A converter's inputs: graph node and the graph
# they return 3 elements:
    # - args:  the names of the tensors that will be needed for the calculations. These are just the names: strings.
    # - layer: a KnetLayer will be constructed. If the weights are in the initializer, the layer will be modified with them.
    # - outs:  the names of the tensors that are the outputs of the calculations. These are just the names: strings.
"""

# ADD
function converter_add(node, g)
    args = node.input
    outs = node.output
    layer = KL.Add()
    return (args, layer, outs)
end

# AveragePool
function converter_avgpool(node, g)
    args = node.input
    outs = node.output
    stride = 0
    padding = 0

    if :pads in keys(node.attribute); padding = node.attribute[:pads][1]; end
    if :strides in keys(node.attribute); stride = node.attribute[:strides][1]; end

    layer = KL.Pool(padding=padding, stride=stride, mode=1)
    (args, layer, outs)
end

# BATCHNORM
function converter_batchnorm(node, g)
    momentum = node.attribute[:momentum]
    epsilon = node.attribute[:epsilon]

    scale = g.initializer[node.input[2]]
    B = g.initializer[node.input[3]]
    mean = g.initializer[node.input[4]]
    variance = g.initializer[node.input[5]]

    X = node.input[1]
    outs = node.output
    layer = KL.BatchNorm(length(scale); momentum=momentum, mean=mean, var=variance)
    (X, layer, outs)
end

# CAST
function converter_cast(node, g)
    args = node.input
    outs = node.output
    to = titlecase(node.attribute[:to]) # ONNX use uppercase but Julia use capitalize on data types
    layer = KL.Cast(to)
    (args, layer, outs)
end

# CONCAT
function converter_concat(node, g)
    args = node.input
    outs = node.output
    axis = node.attribute[:axis]
    layer = KL.Concat(axis)
    (args, layer, outs)
end

# CONSTANT
function converter_constant(node, g)
    args = node.input
    outs = node.output
    onnx_constant_value = node.attribute[:value]
    julia_constant_value = onnx_utils.UInt8toFloat32(onnx_constant_value)
    layer = KL.constant_layer(julia_constant_value)
    return (args, layer, outs)
end

function converter_ConstantOfShape(node, g)
    args = node.input
    outs = node.output
    onnx_constant_value = node.attribute[:value]
    julia_constant_value = onnx_utils.UInt8toFloat32(onnx_constant_value)
    layer = KL.ConstantOfShape(julia_constant_value[1])
    (args, layer, outs)
end

# CONV
#conv1 = KnetOnnx.KnetLayers.Conv(;height=3, width=3, inout = 3=>64)
#currently treating [1,1,1,1] padding as an integer 1, same for stride
function converter_cnn(node, g)
    args = [node.input[1]]
    out = node.output

    padding = 0
    strides = 0
    if :pads in keys(node.attribute); padding = node.attribute[:pads][1]; end
    if :strides in keys(node.attribute); stride = node.attribute[:strides][1]; end

    layer = KnetOnnx.KL.Conv(height=1,width=1,inout=1=>1; padding = padding, stride = stride)

    if length(node.input) >= 2
        w_name = node.input[2]
        w = g.initializer[w_name]
        #might cause a problem later on with different convs
        layer.weight = w

    end
    if length(node.input) >= 3
        b_name = node.input[3]
        b = g.initializer[b_name]
        layer.bias = reshape(b, 1, 1, size(b)[1], 1)
    end
    (args, layer, out)
end

#Div
function converter_div(node, g)
    args = node.input
    outs = node.output
    layer = KL.Div()
    (args, layer, outs)
end

# DROPOUT
function converter_dropout(node, g)
    args = node.input
    outs = node.output
    layer = KL.Dropout(p = node.attribute[:ratio])
    (args, layer, outs)
end

# Exponential
function converter_exp(node, g)
    args = node.input
    outs = node.output
    layer = KL.Exp()
    (args, layer, outs)
end

# FLATTEN
function converter_flatten(node, g)
    args = node.input
    outs = node.output
    layer = KL.Flatten()
    (args, layer, outs)
end

# Gather
function converter_gather(node, g)
    args = node.input
    outs = node.output
    axis = node.attribute[:axis] == nothing ? 1 : node.attribute[:axis] + 1 # +1 is for Julia
    layer = KL.Gather(axis)
    (args, layer, outs)
end

# GEMM - trains bias also for gemm that is not a linear layer, fix that, write new gemm and a separate linear
function converter_gemm(node, g)
    input1 = node.input[1]

    #the layer is a Knet Layer
    layer = KnetOnnx.KnetLayers.Linear(input=1,output=1)

    # use g.initializer to modify KnetLayer
    w_name = node.input[2]
    b_name = length(node.input) >= 3 ? node.input[3] : 0
    w = g.initializer[w_name]
    w = transpose(w)
    b = g.initializer[b_name]

    w = KnetOnnx.KnetLayers.ConvertParams(w)
    b = KnetOnnx.KnetLayers.ConvertParams(b)

    layer.bias = b
    layer.mult.weight = w

    # return input tensor NAMES, it is called args: [input1, ...]
    # you can take the inputs from model.tensors using these names
    args = [input1]
    outs = [node]

    # returns these 3, use these to create ModelLayer
    (args, layer, node.output)
end

# GlobalAveragePool
function converter_globalAveragePool(node, g)
    args = node.input
    layer = KL.Pool()
    outs = node.output
    (args, layer, outs)
end

# Identity
function converter_identity(node, g)
    args = node.input
    layer = identity()
    outs = node.output
    (args, layer, outs)
end

# LEAKY RELU - done
function converter_leakyrelu(node, g)
    args = node.input
    alpha = node.attribute[:alpha]
    layer = KL.LeakyReLU(alpha)
    outs = node.output
    (args, layer, outs)
end

# MatMul
function converter_matmul(node, g)
    args = node.input
    layer = KL.MatMul()
    outs = node.output
    (args, layer, outs)
end

# MaxPool
#currently treating [1,1,1,1] padding as an integer 1, same for
function converter_maxpool(node, g)
    args = node.input
    outs = node.output
    stride = 0
    padding = 0

    if :pads in keys(node.attribute); padding = node.attribute[:pads][1]; end
    if :strides in keys(node.attribute); stride = node.attribute[:strides][1]; end

    layer = KL.Pool(padding=padding, stride=stride, mode=0)
    (args, layer, outs)
end

"""
Mul Converter
"""
function converter_mul(node,g)
    args = node.input
    layer = KL.Mul()
    outs = node.output
    (args, layer, outs)
end

# RELU - done
function converter_relu(node, g)
    args = node.input
    layer = KL.ReLU()
    outs = node.output
    (args, layer, outs)
end

"""
Reshape
"""
function converter_reshape(node, g)
    data = node.input[1]
    shape = node.input[2]
    layer = KL.Reshape(shape)
    outs = node.output
    (args, layer, outs)
end

# SHAPE
function converter_shape(node, g)
    args = node.input
    outs = node.output
    layer = KL.shape_layer()
    (args, layer, outs)
end

"""
SoftMax
The operator computes the softmax (normalized exponential) values for each layer in the batch of the given input.
"""
function converter_softmax(node, g)
    args = node.input
    axis = node.attribute[:axis] == nothing ? 1 : node.attribute[:axis] # default axis = 1
    layer = KL.SoftMax(axis)
    outs = node.output
    (args, layer, outs)
end

"""
Slice
"""
function converter_slice(node, g)
    println("Slice not implemented yet!")
    data = node.input[1]
    starts = node.input[2]
    ends = node.input[3]
    axes = node.input[4]
    steps = node.input[5]
    layer = KL.Slice()
    outs = node.output
    (args, layer, outs)
end

"""
Sub
"""
function converter_sub(node,g)
    args = node.input
    layer = KL.Sub()
    outs = node.output
    (args, layer, outs)
end

# SQUEEZE
function converter_squeeze(node, g)
    args = node.input
    outs = node.output
    layer = KL.squeeze_layer(node.attribute[:axes])
    (args, layer, outs)
end

function converter_transpose(node, g)
    args = node.input
    outs = node.output
    layer = KL.Transpose(node.attribute[:perm] .+ 1) # +1 for julia index
    (args, layer, outs)
end

# UNSQUEEZE
function converter_unsqueeze(node, g)
    args = node.input
    outs = node.output
    layer = KL.unsqueeze_layer(node.attribute[:axes])
    (args, layer, outs)
end

# Node to's

# BATCHNORM
function node_to_batchnorm(node, g)
    momentum = node.attribute[:momentum]
    epsilon = node.attribute[:epsilon]

    scale = g.initializer[node.input[2]]
    B = g.initializer[node.input[3]]
    mean = g.initializer[node.input[4]]
    variance = g.initializer[node.input[5]]

    KL.BatchNorm(length(scale); momentum=momentum, mean=mean, var=variance)
end

# RNN
function node_to_RNN(node, g)
    activations = node.attribute[:activations]
    hidden_size = node.attribute[:hidden_size]
end


text:
module GammaTest

using Test
using ReactiveMP
using Random

import SpecialFunctions: loggamma
import ReactiveMP: xtlog

@testset "Gamma" begin

    @testset "Constructor" begin
        @test GammaShapeRate <: GammaDistributionsFamily
        @test GammaShapeScale <: GammaDistributionsFamily

        @test GammaShapeRate()         == GammaShapeRate{Float64}(1.0, 1.0)
        @test GammaShapeRate(1.0)      == GammaShapeRate{Float64}(1.0, 1.0)
        @test GammaShapeRate(1.0, 2.0) == GammaShapeRate{Float64}(1.0, 2.0)
        @test GammaShapeRate(1)        == GammaShapeRate{Float64}(1.0, 1.0)
        @test GammaShapeRate(1, 2)     == GammaShapeRate{Float64}(1.0, 2.0)
        @test GammaShapeRate(1.0, 2)   == GammaShapeRate{Float64}(1.0, 2.0)
        @test GammaShapeRate(1, 2.0)   == GammaShapeRate{Float64}(1.0, 2.0)
        @test GammaShapeRate(1f0)      == GammaShapeRate{Float32}(1f0, 1f0)
        @test GammaShapeRate(1f0, 2f0) == GammaShapeRate{Float32}(1f0, 2f0)
        @test GammaShapeRate(1f0, 2)   == GammaShapeRate{Float32}(1f0, 2f0)
        @test GammaShapeRate(1f0, 2.0) == GammaShapeRate{Float64}(1.0, 2.0)

        @test GammaShapeScale()         == GammaShapeScale{Float64}(1.0, 1.0)
        @test GammaShapeScale(1.0)      == GammaShapeScale{Float64}(1.0, 1.0)
        @test GammaShapeScale(1.0, 2.0) == GammaShapeScale{Float64}(1.0, 2.0)
        @test GammaShapeScale(1)        == GammaShapeScale{Float64}(1.0, 1.0)
        @test GammaShapeScale(1, 2)     == GammaShapeScale{Float64}(1.0, 2.0)
        @test GammaShapeScale(1.0, 2)   == GammaShapeScale{Float64}(1.0, 2.0)
        @test GammaShapeScale(1, 2.0)   == GammaShapeScale{Float64}(1.0, 2.0)
        @test GammaShapeScale(1f0)      == GammaShapeScale{Float32}(1f0, 1f0)
        @test GammaShapeScale(1f0, 2f0) == GammaShapeScale{Float32}(1f0, 2f0)
        @test GammaShapeScale(1f0, 2)   == GammaShapeScale{Float32}(1f0, 2f0)
        @test GammaShapeScale(1f0, 2.0) == GammaShapeScale{Float64}(1.0, 2.0)

        @test eltype(GammaShapeRate(1.0, 2.0)) === Float64
        @test eltype(GammaShapeRate(1f0, 2f0)) === Float32
    end

    @testset "vague" begin
        vague(GammaShapeScale) == Gamma(1.0, ReactiveMP.huge)
        vague(GammaShapeRate)  == Gamma(1.0, ReactiveMP.tiny)
    end

    @testset "Stats methods for GammaShapeScale" begin

        dist1 = GammaShapeScale(1.0, 1.0)

        @test mean(dist1)     === 1.0
        @test var(dist1)      === 1.0
        @test cov(dist1)      === 1.0
        @test shape(dist1)    === 1.0
        @test scale(dist1)    === 1.0
        @test rate(dist1)     === 1.0
        @test entropy(dist1)    ≈ 1.0

        dist2 = GammaShapeScale(1.0, 2.0)

        @test mean(dist2)     === 2.0
        @test var(dist2)      === 4.0
        @test cov(dist2)      === 4.0
        @test shape(dist2)    === 1.0
        @test scale(dist2)    === 2.0
        @test rate(dist2)     === inv(2.0)
        @test entropy(dist2)    ≈ 1.6931471805599454

        dist3 = GammaShapeScale(2.0, 2.0)

        @test mean(dist3)     === 4.0
        @test var(dist3)      === 8.0
        @test cov(dist3)      === 8.0
        @test shape(dist3)    === 2.0
        @test scale(dist3)    === 2.0
        @test rate(dist3)     === inv(2.0)
        @test entropy(dist3)    ≈ 2.2703628454614764

    end

    @testset "Stats methods for GammaShapeRate" begin

        dist1 = GammaShapeRate(1.0, 1.0)

        @test mean(dist1)     === 1.0
        @test var(dist1)      === 1.0
        @test cov(dist1)      === 1.0
        @test shape(dist1)    === 1.0
        @test scale(dist1)    === 1.0
        @test rate(dist1)     === 1.0
        @test entropy(dist1)     ≈ 1.0
        @test pdf(dist1, 1.0)    ≈ 0.36787944117144233
        @test logpdf(dist1, 1.0) ≈ -1.0

        dist2 = GammaShapeRate(1.0, 2.0)

        @test mean(dist2)     === inv(2.0)
        @test var(dist2)      === inv(4.0)
        @test cov(dist2)      === inv(4.0)
        @test shape(dist2)    === 1.0
        @test scale(dist2)    === inv(2.0)
        @test rate(dist2)     === 2.0
        @test entropy(dist2)     ≈ 0.3068528194400547
        @test pdf(dist2, 1.0)    ≈ 0.2706705664732254
        @test logpdf(dist2, 1.0) ≈ -1.3068528194400546

        dist3 = GammaShapeRate(2.0, 2.0)

        @test mean(dist3)     === 1.0
        @test var(dist3)      === inv(2.0)
        @test cov(dist3)      === inv(2.0)
        @test shape(dist3)    === 2.0
        @test scale(dist3)    === inv(2.0)
        @test rate(dist3)     === 2.0
        @test entropy(dist3)     ≈ 0.8840684843415857
        @test pdf(dist3, 1.0)    ≈ 0.5413411329464508
        @test logpdf(dist3, 1.0) ≈ -0.6137056388801094

    end

    @testset "Base methods" begin
        @test convert(GammaShapeScale{Float32}, GammaShapeScale()) == GammaShapeScale{Float32}(1f0, 1f0)
        @test convert(GammaShapeScale{Float64}, GammaShapeScale(1.0, 10.0)) == GammaShapeScale{Float64}(1.0, 10.0)
        @test convert(GammaShapeScale{Float64}, GammaShapeScale(1.0, 0.1)) == GammaShapeScale{Float64}(1.0, 0.1)
        @test convert(GammaShapeScale{Float64}, 1, 1) == GammaShapeScale{Float64}(1.0, 1.0)
        @test convert(GammaShapeScale{Float64}, 1, 10) == GammaShapeScale{Float64}(1.0, 10.0)
        @test convert(GammaShapeScale{Float64}, 1.0, 0.1) == GammaShapeScale{Float64}(1.0, 0.1)

        @test convert(GammaShapeRate{Float32}, GammaShapeRate()) == GammaShapeRate{Float32}(1f0, 1f0)
        @test convert(GammaShapeRate{Float64}, GammaShapeRate(1.0, 10.0)) == GammaShapeRate{Float64}(1.0, 10.0)
        @test convert(GammaShapeRate{Float64}, GammaShapeRate(1.0, 0.1)) == GammaShapeRate{Float64}(1.0, 0.1)
        @test convert(GammaShapeRate{Float64}, 1, 1) == GammaShapeRate{Float64}(1.0, 1.0)
        @test convert(GammaShapeRate{Float64}, 1, 10) == GammaShapeRate{Float64}(1.0, 10.0)
        @test convert(GammaShapeRate{Float64}, 1.0, 0.1) == GammaShapeRate{Float64}(1.0, 0.1)

        @test convert(GammaShapeRate, GammaShapeRate(2.0, 2.0)) == GammaShapeRate{Float64}(2.0, 2.0)
        @test convert(GammaShapeScale, GammaShapeRate(2.0, 2.0)) == GammaShapeScale{Float64}(2.0, 1.0 / 2.0)

        @test convert(GammaShapeRate, GammaShapeScale(2.0, 2.0)) == GammaShapeRate{Float64}(2.0, 1.0 / 2.0)
        @test convert(GammaShapeScale, GammaShapeScale(2.0, 2.0)) == GammaShapeScale{Float64}(2.0, 2.0)

        check_basic_statistics = (left, right) -> begin
            @test mean(left)         ≈ mean(right)
            @test var(left)          ≈ var(right)
            @test cov(left)          ≈ cov(right)
            @test shape(left)        ≈ shape(right)
            @test scale(left)        ≈ scale(right)
            @test rate(left)         ≈ rate(right)
            @test entropy(left)      ≈ entropy(right)
            @test pdf(left, 1.0)     ≈ pdf(right, 1.0)
            @test pdf(left, 10.0)    ≈ pdf(right, 10.0)
            @test logpdf(left, 1.0)  ≈ logpdf(right, 1.0)
            @test logpdf(left, 10.0) ≈ logpdf(right, 10.0)
            @test mean(log, left)    ≈ mean(log, right)
            @test mean(loggamma, left) ≈ mean(loggamma, right)
            @test mean(xtlog, left)    ≈ mean(xtlog, right)
        end

        types = ReactiveMP.union_types(GammaDistributionsFamily{Float64})
        rng   = MersenneTwister(1234)

        for type in types
            left = convert(type, rand(rng, Float64), rand(rng, Float64))
            for type in types
                right = convert(type, left)
                check_basic_statistics(left, right)
            end
        end
    end

    @testset "prod" begin

        @test prod(ProdAnalytical(), GammaShapeScale(1, 1), GammaShapeScale(1, 1)) == GammaShapeScale(1, 1/2)
        @test prod(ProdAnalytical(), GammaShapeScale(1, 2), GammaShapeScale(1, 1)) == GammaShapeScale(1, 2/3)
        @test prod(ProdAnalytical(), GammaShapeScale(1, 2), GammaShapeScale(1, 2)) == GammaShapeScale(1, 1)
        @test prod(ProdAnalytical(), GammaShapeScale(2, 2), GammaShapeScale(1, 2)) == GammaShapeScale(2, 1)
        @test prod(ProdAnalytical(), GammaShapeScale(2, 2), GammaShapeScale(2, 2)) == GammaShapeScale(3, 1)

        @test prod(ProdAnalytical(), GammaShapeRate(1, 1), GammaShapeRate(1, 1)) == GammaShapeRate(1, 2)
        @test prod(ProdAnalytical(), GammaShapeRate(1, 2), GammaShapeRate(1, 1)) == GammaShapeRate(1, 3)
        @test prod(ProdAnalytical(), GammaShapeRate(1, 2), GammaShapeRate(1, 2)) == GammaShapeRate(1, 4)
        @test prod(ProdAnalytical(), GammaShapeRate(2, 2), GammaShapeRate(1, 2)) == GammaShapeRate(2, 4)
        @test prod(ProdAnalytical(), GammaShapeRate(2, 2), GammaShapeRate(2, 2)) == GammaShapeRate(3, 4)

        @test prod(ProdAnalytical(), GammaShapeScale(1, 1), GammaShapeRate(1, 1)) == GammaShapeScale(1, 1/2)
        @test prod(ProdAnalytical(), GammaShapeScale(1, 2), GammaShapeRate(1, 1)) == GammaShapeScale(1, 2/3)
        @test prod(ProdAnalytical(), GammaShapeScale(1, 2), GammaShapeRate(1, 2)) == GammaShapeScale(1, 2/5)
        @test prod(ProdAnalytical(), GammaShapeScale(2, 2), GammaShapeRate(1, 2)) == GammaShapeScale(2, 2/5)
        @test prod(ProdAnalytical(), GammaShapeScale(2, 2), GammaShapeRate(2, 2)) == GammaShapeScale(3, 2/5)

        @test prod(ProdAnalytical(), GammaShapeRate(1, 1), GammaShapeScale(1, 1)) == GammaShapeRate(1, 2)
        @test prod(ProdAnalytical(), GammaShapeRate(1, 2), GammaShapeScale(1, 1)) == GammaShapeRate(1, 3)
        @test prod(ProdAnalytical(), GammaShapeRate(1, 2), GammaShapeScale(1, 2)) == GammaShapeRate(1, 5/2)
        @test prod(ProdAnalytical(), GammaShapeRate(2, 2), GammaShapeScale(1, 2)) == GammaShapeRate(2, 5/2)
        @test prod(ProdAnalytical(), GammaShapeRate(2, 2), GammaShapeScale(2, 2)) == GammaShapeRate(3, 5/2)

    end

end

end

ext: jl


text:
mutable struct Context <: AContext
  cache::Union{IdDict{Any,Any},Nothing}
  globals::Union{Dict{GlobalRef,Any},Nothing}
end

Context() = Context(nothing, nothing)

cache(cx::Context) = cx.cache === nothing ? (cx.cache = IdDict()) : cx.cache
globals(cx::Context) = cx.globals === nothing ? (cx.globals = Dict{GlobalRef,Any}()) : cx.globals

struct Pullback{S,T}
  t::T
end

Pullback{S}(x) where S = Pullback{S,typeof(x)}(x)

struct CompileError
  T
  e
end

function Base.showerror(io::IO, e::CompileError)
  print(io, "Compiling $(e.T): ")
  showerror(io, e.e)
end

# interface2.jl

# Wrappers

_pullback(f, args...) = _pullback(Context(), f, args...)

tailmemaybe(::Nothing) = nothing
tailmemaybe(x::Tuple) = Base.tail(x)

function pullback(f, args...)
  y, back = _pullback(f, args...)
  y, Δ -> tailmemaybe(back(Δ))
end

sensitivity(y::Number) = one(y)
sensitivity(y::Complex) = error("Output is complex, so the gradient is not defined.")
sensitivity(y) = error("Output should be scalar; gradients are not defined for output $(repr(y))")

function gradient(f, args...)
  y, back = pullback(f, args...)
  return back(sensitivity(y))
end

Base.adjoint(f::Function) = x -> gradient(f, x)[1]

# Param-style wrappers

# TODO store ids only
struct Params
  order::Buffer{Any, Vector{Any}}
  params::IdSet{Any}
  Params() = new(Buffer([], false), IdSet())
end

@forward Params.order Base.iterate, Base.length

function Base.push!(ps::Params, x)
  if !(x in ps.params)
    push!(ps.order, x)
    push!(ps.params, x)
  end
  return ps
end

Base.push!(ps::Params, x...) = (foreach(x -> push!(ps, x), x); ps)

Params(xs) = push!(Params(), xs...)

function Base.show(io::IO, ps::Params)
  print(io, "Params([")
  join(io, ps.order, ", ")
  print(io, "])")
end

struct Grads
  grads::IdDict{Any,Any}
end

Base.show(io::IO, ps::Grads) = print(io, "Grads(...)")

@forward Grads.grads Base.getindex, Base.haskey

function Base.getindex(gs::Grads, x)
  isbits(x) && error("Only reference types can be differentiated with `Params`.")
  return gs.grads[x]
end

function pullback(f, ps::Params)
  cx = Context()
  y, back = _pullback(cx, f)
  y, function (Δ)
    for p in ps
      cache(cx)[p] = nothing
    end
    back(Δ)
    Grads(cx.cache) # TODO make a copy
  end
end

# Code Reflection

using InteractiveUtils
using InteractiveUtils: typesof
using Core: Typeof

function code_ir(f, T)
  m = meta(Tuple{Typeof(f),T.parameters...})
  return IR(m)
end

function code_irm(ex)
  isexpr(ex, :call) || error("@code_ir f(args...)")
  f, args = ex.args[1], ex.args[2:end]
  :(code_ir($(esc(f)), typesof($(esc.(args)...))))
end

macro code_ir(ex)
  code_irm(ex)
end

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


text:
# Thin wrapper type, mainly for dispatch
struct FullBinner{T, A <: AbstractVector{T}} <: AbstractBinner{T}
    x::A
end

FullBinner(::Type{T} = Float64) where T = FullBinner(Vector{T}(undef, 0))

@forward FullBinner.x (Base.length, Base.size, Base.lastindex, Base.ndims, Base.iterate,
                        Base.getindex, Base.setindex!, Base.view, Base.axes,
                        Base.resize!, Base.sizehint!, Base.empty!, Base.isempty)

Base.:(!=)(a::FullBinner, b::FullBinner) = a.x != b.x
Base.:(==)(a::FullBinner, b::FullBinner) = a.x == b.x
function Base.isapprox(a::FullBinner, b::FullBinner; kwargs...)
    (length(a) == length(b)) && (isempty(a) || isapprox(a.x, b.x; kwargs...))
end

# Cosmetics

function _print_header(io::IO, B::FullBinner{T,A}) where {T, A}
    print(io, "FullBinner{$(T),$(A)}")
    nothing
end

function _println_body(io::IO, B::FullBinner{T,A}) where {T, A}
    n = length(B)
    println(io)
    print(io, "| Count: ", n)
    if n > 0 && T <: Number
        print(io, "\n| Mean: ", round.(mean(B), digits=5))
        # print(io, "\n| StdError: ", round.(std_error(B), digits=5))
    end
    nothing
end

# short version (shows up in arrays etc.)
Base.show(io::IO, B::FullBinner{T,A}) where {T, A} = print(io, "FullBinner{$(T),$(A)}()")
# verbose version (shows up in the REPL)
Base.show(io::IO, m::MIME"text/plain", B::FullBinner) = (_print_header(io, B); _println_body(io, B))

# We explicitly copy here to avoid falsifying data when a pushed array is
# overwritten. (For example one might have a temporary array which accumulates
# data before pushing it to a FullBinner. If this array is reused for something
# else and we do not copy here, the pushed element of the FullBinner changes.
# This is not desirable.)
Base.push!(B::FullBinner, x) = push!(B.x, deepcopy(x))
Base.append!(B::FullBinner{<: Number}, x::AbstractArray{<: Number}) = append!(B.x, x)

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


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

## 3D tests
include("../scripts/diff_3D/diff_3D_nonlin_perf.jl")
@reset_parallel_stencil()
indsx = Int.(ceil.(LinRange(1, length(xc), 12)))
indsy = Int.(ceil.(LinRange(1, length(yc), 12)))
indsz = Int.(ceil.(LinRange(1, length(zc), 12)))
d3d   = Dict(:X=> xc[indsx], :H=>H[indsx,indsy,indsz])

@testset "Reference-tests diffusion 3D nonlin perf (noMPI)" begin
    @test_reference "reftest-files/test_diff_3D_nonlin_perf.bson" d3d by=comp
end

@reset_parallel_stencil()


text:
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Operators: Ax_ψᵃᵃᶜ, Ay_ψᵃᵃᶜ

#####
##### Momentum flux operators
#####

@inline momentum_flux_huu(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_uu(i, j, k, grid, advection, solution.uh, solution.uh) / solution.h[i, j, k]

@inline momentum_flux_huv(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_uv(i, j, k, grid, advection, solution.vh, solution.uh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)

@inline momentum_flux_hvu(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_vu(i, j, k, grid, advection, solution.uh, solution.vh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)

@inline momentum_flux_hvv(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_vv(i, j, k, grid, advection, solution.vh, solution.vh) / solution.h[i, j, k]

#####
##### Momentum flux divergence operators
#####

@inline div_hUu(i, j, k, grid, advection, solution) =
    1 / Vᵃᵃᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, momentum_flux_huu, advection, solution) +
                               δyᵃᶜᵃ(i, j, k, grid, momentum_flux_huv, advection, solution))

@inline div_hUv(i, j, k, grid, advection, solution) =
    1 / Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, momentum_flux_hvu, advection, solution) +
                               δyᵃᶠᵃ(i, j, k, grid, momentum_flux_hvv, advection, solution))

# Support for no advection
@inline div_hUu(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution) where FT = zero(FT)
@inline div_hUv(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution) where FT = zero(FT)

#####
##### Mass transport divergence operator
#####

"""
    div_uhvh(i, j, k, grid, solution)

Calculates the divergence of the mass flux into a cell,

    1/V * [δxᶜᵃᵃ(Ax * uh) + δyᵃᶜᵃ(Ay * vh)]

which will end up at the location `ccc`.
"""
@inline function div_uhvh(i, j, k, grid, solution)
    1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶜ, solution.uh) +
                             δyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶜ, solution.vh))
end


text:
"""
    AbstractInfiniteSet

An abstract type for sets that characterize infinite parameters.
"""
abstract type AbstractInfiniteSet end

"""
    InfOptParameter{T <: AbstractInfiniteSet} <: JuMP.AbstractVariable

A DataType for storing core infinite parameter information.

**Fields**
- `set::T` The infinite set that characterizes the parameter.
- `supports::Vector{<:Number}` The support points used to discretize this
                               parameter.
- `independent::Bool` Is independent of other parameters that share its group ID
                      number.
"""
struct InfOptParameter{T <: AbstractInfiniteSet} <: JuMP.AbstractVariable
    set::T
    supports::Vector{<:Number}
    independent::Bool
end

"""
    InfOptVariable <: JuMP.AbstractVariable

An abstract type for infinite, point, and hold variables.
"""
abstract type InfOptVariable <: JuMP.AbstractVariable end

"""
    AbstractReducedInfo

An abstract type reduced variable information.
"""
abstract type AbstractReducedInfo end

"""
    AbstractMeasureData

An abstract type to define data for measures to define the behavior of
[`Measure`](@ref).
"""
abstract type AbstractMeasureData end

"""
    Measure{T <: JuMP.AbstractJuMPScalar, V <: AbstractMeasureData}

A DataType for measure abstractions.

**Fields**
- `func::T` Infinite variable expression.
- `data::V` Data of the abstraction as described in a `AbstractMeasureData`
            subtype.
"""
struct Measure{T <: JuMP.AbstractJuMPScalar, V <: AbstractMeasureData}
    func::T
    data::V
end

"""
    InfiniteModel <: JuMP.AbstractModel

A DataType for storing all of the mathematical modeling information needed to
model an optmization problem with an infinite dimensional decision space.

**Fields**
- `next_meas_index::Int` Index - 1 of next measure.
- `measures::Dict{Int, Measure}` Measure indices to measure datatypes.
- `meas_to_name::Dict{Int, String}` Measure indices to names.
- `meas_to_constrs::Dict{Int, Vector{Int}}` Measure indices to dependent
                                            constraint indices.
- `meas_to_meas::Dict{Int, Vector{Int}}` Measure indices to dependent
                                         measure indices.
- `meas_in_objective::Dict{Int, Bool}` Measure indices to if used in objective.
- `integral_defaults::Dict{Symbol}` Default keyword argument settings for measures.
- `next_param_index::Int` Index - 1 of next infinite parameter.
- `next_param_id::Int` Index - 1 of the next infinite parameter group.
- `params::Dict{Int, InfOptParameter}` Infinite parameter indices to parameter
                                       datatype.
- `param_to_name::Dict{Int, String}` Infinite parameter indices to names.
- `name_to_param::Union{Dict{String, Int}, Nothing}` Names to infinite
                                                     parameters.
- `param_to_group_id::Dict{Int, Int}` Infinite parameter indices to group IDs.
- `param_to_constrs::Dict{Int, Vector{Int}}` Infinite parameter indices to list
                                             of dependent constraint indices.
- `param_to_meas::Dict{Int, Vector{Int}}` Infinite parameter indices to list
                                          of dependent measure indices.
- `param_to_vars::Dict{Int, Vector{Int}}` Infinite parameter indices to list
                                          of dependent variable indices.
- `next_var_index::Int` Index - 1 of next variable index.
- `vars::Dict{Int, Dict{Int, Union{InfOptVariable, ReducedVariable}}` Variable
                                                  indices to variable datatype.
- `var_to_name::Dict{Int, String}` Variable indices to names.
- `name_to_var::Union{Dict{String, Int}, Nothing}` Variable names to indices.
- `var_to_lower_bound::Dict{Int, Int}` Variable indices to lower bound index.
- `var_to_upper_bound::Dict{Int, Int}` Variable indices to upper bound index.
- `var_to_fix::Dict{Int, Int}` Variable indices to fix index.
- `var_to_zero_one::Dict{Int, Int}` Variable indices to binary index.
- `var_to_integrality::Dict{Int, Int}` Variable indices to integer index.
- `var_to_constrs::Dict{Int, Vector{Int}}` Variable indices to dependent
                                           constraint indices.
- `var_to_meas::Dict{Int, Vector{Int}}` Variable indices to dependent
                                        measure indices.
- `var_in_objective::Dict{Int, Bool}` Variable indices to if used in objective.
- `infinite_to_points::Dict{Int, Vector{Int}}` Infinite variable indices to
                                               dependent point variable indices.
- `infinite_to_reduced::Dict{Int, Vector{Int}}` Infinite variable indices to
                                               dependent reduced variable indices.
- `has_hold_bounds::Bool` Have hold variables with bounds been added to the model
- `reduced_to_constrs::Dict{Int, Vector{Int}}` Reduced variable indices to dependent
                                               constraint indices.
- `reduced_to_meas::Dict{Int, Vector{Int}}` Reduced variable indices to dependent
                                            measure indices.
- `reduced_info::Dict{Int, AbstractReducedInfo}` Reduced variable indices to
                                                 reduced variable information.
- `next_constr_index::Int` Index - 1 of next constraint.
- `constrs::Dict{Int, JuMP.AbstractConstraint}` Constraint indices to constraint
                                                datatypes.
- `constr_to_name::Dict{Int, String}` Constraint indices to names.
- `name_to_constr::Union{Dict{String, Int}, Nothing}` Constraint names to
                                                      indices.
- `constr_in_var_info::Dict{Int, Bool}` Constraint indices to if related to
                                        variable information constraints.
- `objective_sense::MOI.OptimizationSense` Objective sense.
- `objective_function::JuMP.AbstractJuMPScalar` Finite scalar function.
- `obj_dict::Dict{Symbol, Any}` Store Julia symbols used with `InfiniteModel`
- `optimizer_constructor` MOI optimizer constructor (e.g., Gurobi.Optimizer).
- `optimizer_model::JuMP.Model` Model used to solve `InfiniteModel`
- `ready_to_optimize::Bool` Is the optimizer_model up to date.
- `ext::Dict{Symbol, Any}` Store arbitrary extension information.
"""
mutable struct InfiniteModel <: JuMP.AbstractModel
    # Measure Data
    next_meas_index::Int
    measures::Dict{Int, Measure}
    meas_to_name::Dict{Int, String}
    meas_to_constrs::Dict{Int, Vector{Int}}
    meas_to_meas::Dict{Int, Vector{Int}}
    meas_in_objective::Dict{Int, Bool}
    integral_defaults::Dict{Symbol}

    # Parameter Data
    next_param_index::Int
    next_param_id::Int
    params::Dict{Int, InfOptParameter}
    param_to_name::Dict{Int, String}
    name_to_param::Union{Dict{String, Int}, Nothing}
    param_to_group_id::Dict{Int, Int}
    param_to_constrs::Dict{Int, Vector{Int}}
    param_to_meas::Dict{Int, Vector{Int}}
    param_to_vars::Dict{Int, Vector{Int}}

    # Variable data
    next_var_index::Int
    vars::Dict{Int, InfOptVariable}
    var_to_name::Dict{Int, String}
    name_to_var::Union{Dict{String, Int}, Nothing}
    var_to_lower_bound::Dict{Int, Int}
    var_to_upper_bound::Dict{Int, Int}
    var_to_fix::Dict{Int, Int}
    var_to_zero_one::Dict{Int, Int}
    var_to_integrality::Dict{Int, Int}
    var_to_constrs::Dict{Int, Vector{Int}}
    var_to_meas::Dict{Int, Vector{Int}}
    var_in_objective::Dict{Int, Bool}
    infinite_to_points::Dict{Int, Vector{Int}}
    infinite_to_reduced::Dict{Int, Vector{Int}}
    has_hold_bounds::Bool

    # Placeholder
    reduced_to_constrs::Dict{Int, Vector{Int}}
    reduced_to_meas::Dict{Int, Vector{Int}}
    reduced_info::Dict{Int, AbstractReducedInfo}

    # Constraint Data
    next_constr_index::Int
    constrs::Dict{Int, JuMP.AbstractConstraint}
    constr_to_name::Dict{Int, String}
    name_to_constr::Union{Dict{String, Int}, Nothing}
    constr_in_var_info::Dict{Int, Bool}

    # Objective Data
    objective_sense::MOI.OptimizationSense
    objective_function::JuMP.AbstractJuMPScalar

    # Objects
    obj_dict::Dict{Symbol, Any}

    # Optimize Data
    optimizer_constructor
    optimizer_model::JuMP.Model
    ready_to_optimize::Bool

    # Extensions
    ext::Dict{Symbol, Any}
end

"""
    InfiniteModel([optimizer_constructor; seed::Bool = false,
                  OptimizerModel::Function = TranscriptionModel,
                  caching_mode::MOIU.CachingOptimizerMode = MOIU.AUTOMATIC,
                  bridge_constraints::Bool = true, optimizer_model_kwargs...])

Return a new infinite model where an optimizer is specified if an
`optimizer_constructor` is given. The `seed` argument indicates if the stochastic
sampling used in conjunction with the model should be seeded. The optimizer
can also later be set with the [`JuMP.set_optimizer`](@ref) call. By default
the `optimizer_model` data field is initialized with a
[`TranscriptionModel`](@ref), but a different type of model can be assigned via
[`set_optimizer_model`](@ref) as can be required by extensions.

**Example**
```jldoctest
julia> using InfiniteOpt, JuMP, Ipopt;

julia> model = InfiniteModel()
An InfiniteOpt Model
Feasibility problem with:
Variables: 0
Optimizer model backend information:
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> model = InfiniteModel(Ipopt.Optimizer)
An InfiniteOpt Model
Feasibility problem with:
Variables: 0
Optimizer model backend information:
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt
```
"""
function InfiniteModel(; seed::Bool = false,
                       OptimizerModel::Function = TranscriptionModel, kwargs...)
    if seed
        Random.seed!(0)
    end
    return InfiniteModel(# Measures
                         0, Dict{Int, Measure}(), Dict{Int, String}(),
                         Dict{Int, Vector{Int}}(), Dict{Int, Vector{Int}}(),
                         Dict{Int, Bool}(),
                         Dict(:eval_method => sampling,
                              :num_supports => 10,
                              :weight_func => _w,
                              :name => "integral",
                              :use_existing_supports => false),
                         # Parameters
                         0, 0, Dict{Int, InfOptParameter}(), Dict{Int, String}(),
                         nothing, Dict{Int, Int}(), Dict{Int, Vector{Int}}(),
                         Dict{Int, Vector{Int}}(), Dict{Int, Vector{Int}}(),
                         # Variables
                         0, Dict{Int, JuMP.AbstractVariable}(),
                         Dict{Int, String}(), nothing, Dict{Int, Int}(),
                         Dict{Int, Int}(), Dict{Int, Int}(), Dict{Int, Int}(),
                         Dict{Int, Int}(), Dict{Int, Vector{Int}}(),
                         Dict{Int, Vector{Int}}(), Dict{Int, Bool}(),
                         Dict{Int, Vector{Int}}(), Dict{Int, Vector{Int}}(),
                         false,
                         # Placeholder variables
                         Dict{Int, Vector{Int}}(), Dict{Int, Vector{Int}}(),
                         Dict{Int, Tuple}(),
                         # Constraints
                         0, Dict{Int, JuMP.AbstractConstraint}(),
                         Dict{Int, String}(), nothing, Dict{Int, Bool}(),
                         # Objective
                         MOI.FEASIBILITY_SENSE,
                         zero(JuMP.GenericAffExpr{Float64, FiniteVariableRef}),
                         # Object dictionary
                         Dict{Symbol, Any}(),
                         # Optimize data
                         nothing, OptimizerModel(;kwargs...), false,
                         # Extensions
                         Dict{Symbol, Any}())
end

## Set the optimizer_constructor depending on what it is
# MOI.OptimizerWithAttributes
function _set_optimizer_constructor(model::InfiniteModel,
                                    constructor::MOI.OptimizerWithAttributes)
    model.optimizer_constructor = constructor.optimizer_constructor
    return
end

# No attributes
function _set_optimizer_constructor(model::InfiniteModel, constructor)
    model.optimizer_constructor = constructor
    return
end

# Dispatch for InfiniteModel call with optimizer constructor
function InfiniteModel(optimizer_constructor; seed::Bool = false,
                       OptimizerModel::Function = TranscriptionModel, kwargs...)
    model = InfiniteModel(seed = seed)
    model.optimizer_model = OptimizerModel(optimizer_constructor; kwargs...)
    _set_optimizer_constructor(model, optimizer_constructor)
    return model
end

# Define basic InfiniteModel extensions
Base.broadcastable(model::InfiniteModel) = Ref(model)
JuMP.object_dictionary(model::InfiniteModel) = model.obj_dict

"""
    GeneralVariableRef <: JuMP.AbstractVariableRef

An abstract type to for variable references used with infinite models.
"""
abstract type GeneralVariableRef <: JuMP.AbstractVariableRef end

"""
    MeasureFiniteVariableRef <: GeneralVariableRef

An abstract type to define finite variable and measure references.
"""
abstract type MeasureFiniteVariableRef <: GeneralVariableRef end

"""
    FiniteVariableRef <: GeneralVariableRef

An abstract type to define new finite variable references.
"""
abstract type FiniteVariableRef <: MeasureFiniteVariableRef end

"""
    HoldVariableRef <: FiniteVariableRef

A DataType for finite fixed variable references (e.g., first stage variables,
steady-state variables).

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of variable in model.
"""
struct HoldVariableRef <: FiniteVariableRef
    model::InfiniteModel
    index::Int
end

"""
    PointVariableRef <: FiniteVariableRef

A DataType for variables defined at a transcipted point (e.g., second stage
variable at a particular scenario, dynamic variable at a discretized time point).

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of variable in model.
"""
struct PointVariableRef <: FiniteVariableRef
    model::InfiniteModel
    index::Int
end

"""
    InfiniteVariableRef <: GeneralVariableRef

A DataType for untranscripted infinite dimensional variable references (e.g.,
second stage variables, time dependent variables).

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of variable in model.
"""
struct InfiniteVariableRef <: GeneralVariableRef
    model::InfiniteModel # `model` owning the variable
    index::Int           # Index in `model.variables`
end

"""
    ReducedInfiniteVariableRef <: GeneralVariableRef

A DataType for partially transcripted infinite dimensional variable references.
This is used to expand measures that contain infinite variables that are not
fully transcripted by the measure.

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of variable in model.
"""
struct ReducedInfiniteVariableRef <: GeneralVariableRef
    model::InfiniteModel
    index::Int
end

"""
    ParameterRef <: GeneralVariableRef

A DataType for untranscripted infinite parameters references that parameterize
the infinite variables.

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of variable in model.
"""
struct ParameterRef <: GeneralVariableRef
    model::InfiniteModel
    index::Int
end

"""
    MeasureRef <: FiniteVariableRef

A DataType for referring to measure abstractions.

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of variable in model.
"""
struct MeasureRef <: MeasureFiniteVariableRef
    model::InfiniteModel
    index::Int
end

"""
    IntervalSet <: AbstractInfiniteSet

A DataType that stores the lower and upper interval bounds for infinite
parameters that are continuous over a certain that interval.

**Fields**
- `lower_bound::Float64` Lower bound of the infinite parameter.
- `upper_bound::Float64` Upper bound of the infinite parameter.
"""
struct IntervalSet <: AbstractInfiniteSet
    lower_bound::Float64
    upper_bound::Float64
    function IntervalSet(lower::Float64, upper::Float64)
        if lower > upper
            error("Invalid interval set bounds, lower bound is greater than " *
                  "upper bound.")
        end
        return new(lower, upper)
    end
end

"""
    IntervalSet(lower_bound::Number, upper_bound::Number)

A constructor for [`IntervalSet`](@ref) that converts values of type `Number` to
values of type `Float64` as required by `IntervalSet`.
"""
IntervalSet(lb::Number, ub::Number) = IntervalSet(convert(Float64, lb),
                                                  convert(Float64, ub))

"""
    DistributionSet{T <: Distributions.NonMatrixDistribution} <: AbstractInfiniteSet

A DataType that stores the distribution characterizing infinite parameters that
are random.

**Fields**
- `distribution::T` Distribution of the random parameter.
"""
struct DistributionSet{T <: Distributions.NonMatrixDistribution} <: AbstractInfiniteSet
    distribution::T
end

"""
    InfiniteVariable{S, T, U, V} <: InfOptVariable
A DataType for storing core infinite variable information. Note each element of
the parameter reference tuple must contain either a single
[`ParameterRef`](@ref) or an `AbstractArray` of `ParameterRef`s where each
`ParameterRef` has the same group ID number.

**Fields**
- `info::JuMP.VariableInfo{S, T, U, V}` JuMP variable information.
- `parameter_refs::Tuple` The infinite parameters(s) that parameterize the
                          variable.
"""
struct InfiniteVariable{S, T, U, V} <: InfOptVariable
    info::JuMP.VariableInfo{S, T, U, V}
    parameter_refs::Tuple
end

"""
    PointVariable{S, T, U, V} <: InfOptVariable
A DataType for storing point variable information. Note that the elements
`parameter_values` field must match the format of the parameter reference tuple
defined in [`InfiniteVariable`](@ref)

**Fields**
- `info::JuMP.VariableInfo{S, T, U, V}` JuMP Variable information.
- `infinite_variable_ref::InfiniteVariableRef` The infinite variable associated
                                               with the point variable.
- `parameter_values::Tuple` The infinite parameter values defining the point.
"""
struct PointVariable{S, T, U, V} <: InfOptVariable
    info::JuMP.VariableInfo{S, T, U, V}
    infinite_variable_ref::InfiniteVariableRef
    parameter_values::Tuple
end

## Modify parameter dictionary to expand any multidimensional parameter keys
# Case where dictionary is already in correct form
function _expand_parameter_dict(param_bounds::Dict{ParameterRef,
                                                   IntervalSet})::Dict
    return param_bounds
end

# Case where dictionary contains vectors
function _expand_parameter_dict(param_bounds::Dict{<:Any, IntervalSet})::Dict
    # Initialize new dictionary
    new_dict = Dict{ParameterRef, IntervalSet}()
    # Find vector keys and expand
    for (key, set) in param_bounds
        # expand over the array of parameters if this is
        if isa(key, AbstractArray)
            for param in values(key)
                new_dict[param] = set
            end
        # otherwise we have parameter reference
        else
            new_dict[key] = set
        end
    end
    return new_dict
end

# Case where dictionary contains vectors
function _expand_parameter_dict(param_bounds::Dict)
    error("Invalid parameter bound dictionary format.")
end

"""
    ParameterBounds
A DataType for storing intervaled bounds of parameters. This is used to define
subdomains of [`HoldVariable`](@ref)s and [`BoundedScalarConstraint`](@ref)s.

**Fields**
- `intervals::Dict{ParameterRef, IntervalSet}` A dictionary of parameter intervals
that are tighter than those already associated those paraticular parameters.
"""
struct ParameterBounds
    intervals::Dict{ParameterRef, IntervalSet}
    function ParameterBounds(intervals::Dict)
        return new(_expand_parameter_dict(intervals))
    end
end

# Default method
function ParameterBounds()
    return ParameterBounds(Dict{ParameterRef, IntervalSet}())
end

"""
    HoldVariable{S, T, U, V} <: InfOptVariable
A DataType for storing hold variable information.

**Fields**
- `info::JuMP.VariableInfo{S, T, U, V}` JuMP variable information.
- `parameter_bounds::ParameterBounds` Valid parameter sub-domains
"""
struct HoldVariable{S, T, U, V} <: InfOptVariable
    info::JuMP.VariableInfo{S, T, U, V}
    parameter_bounds::ParameterBounds
end

"""
    ReducedInfiniteInfo <: AbstractReducedInfo

A DataType for storing reduced infinite variable information.

**Fields**
- `infinite_variable_ref::InfiniteVariableRef` The original infinite variable.
- `eval_supports::Dict{Int, Union{Number, JuMP.Containers.SparseAxisArray{<:Number}}}`
  The original parameter tuple indices to the evaluation supports.
"""
struct ReducedInfiniteInfo <: AbstractReducedInfo
    infinite_variable_ref::InfiniteVariableRef
    eval_supports::Dict{Int, Union{Number,
                                   JuMPC.SparseAxisArray{<:Number}}}
end

"""
    InfOptVariableRef

A union type for infinite, point, and hold variable references.
"""
const InfOptVariableRef = Union{InfiniteVariableRef, PointVariableRef,
                                HoldVariableRef}

# Define infinite expressions
const InfiniteExpr = Union{InfiniteVariableRef, ReducedInfiniteVariableRef,
                           JuMP.GenericAffExpr{Float64, InfiniteVariableRef},
                           JuMP.GenericAffExpr{Float64, GeneralVariableRef},
                           JuMP.GenericAffExpr{Float64, ReducedInfiniteVariableRef},
                           JuMP.GenericQuadExpr{Float64, InfiniteVariableRef},
                           JuMP.GenericQuadExpr{Float64, GeneralVariableRef},
                           JuMP.GenericQuadExpr{Float64, ReducedInfiniteVariableRef}}
const ParameterExpr = Union{ParameterRef,
                            JuMP.GenericAffExpr{Float64, ParameterRef},
                            JuMP.GenericQuadExpr{Float64, ParameterRef}}

"""
    DiscreteMeasureData <: AbstractMeasureData

A DataType for one dimensional measure abstraction data where the measure
abstraction is of the form:
``measure = \\int_{\\tau \\in T} f(\\tau) w(\\tau) d\\tau \\approx \\sum_{i = 1}^N \\alpha_i f(\\tau_i) w(\\tau_i)``.

**Fields**
- `parameter_ref::ParameterRef` The infinite parameter over which the
                                integration occurs.
- `coefficients::Vector{<:Number}` Coefficients ``\\alpha_i`` for the above
                                   measure abstraction.
- `supports::Vector{<:Number}` Support points ``\\tau_i`` for the above
                               measure abstraction.
- `name::String` Name of the measure that will be implemented.
- `weight_function::Function` Weighting function ``w`` must map support value
                              input value of type `Number` to a scalar value.
"""
struct DiscreteMeasureData <: AbstractMeasureData
    parameter_ref::ParameterRef
    coefficients::Vector{<:Number}
    supports::Vector{<:Number}
    name::String
    weight_function::Function
    # Constructor
    function DiscreteMeasureData(parameter_ref::ParameterRef,
                                 coeffs::Vector{<:Number},
                                 supports::Vector{<:Number},
                                 name::String, weight_func::Function)
        if length(coeffs) != length(supports)
            error("The amount of coefficients must match the amount of " *
                  "support points.")
        end
        set = infinite_set(parameter_ref)
        # set = parameter_ref.model.params[parameter_ref.index].set
        if !supports_in_set(supports, set)
            error("Support points violate parameter domain.")
        end
        return new(parameter_ref, coeffs, supports, name, weight_func)
    end
end

"""
    MultiDiscreteMeasureData <: AbstractMeasureData

A DataType for multi-dimensional measure abstraction data where the measure
abstraction is of the form:
``measure = \\int_{\\tau \\in T} f(\\tau) w(\\tau) d\\tau \\approx \\sum_{i = 1}^N \\alpha_i f(\\tau_i) w(\\tau_i)``.

**Fields**
- `parameter_ref::JuMP.Containers.SparseAxisArray{<:ParameterRef}` The infinite
   parameters over which the integration occurs.
- `coefficients::Vector{<:Number}` Coefficients ``\\alpha_i`` for the above
                                   measure abstraction.
- `supports::Vector{<:JuMP.Containers.SparseAxisArray{<:Number}}` Support points
   ``\\tau_i`` for the above measure abstraction.
- `name::String` Name of the measure that will be implemented.
- `weight_function::Function` Weighting function ``w`` must map a numerical
                              support of type `JuMP.Containers.SparseAxisArray`
                              to a scalar value.
"""
struct MultiDiscreteMeasureData <: AbstractMeasureData
    parameter_ref::JuMPC.SparseAxisArray{<:ParameterRef}
    coefficients::Vector{<:Number}
    supports::Vector{<:JuMPC.SparseAxisArray}
    name::String
    weight_function::Function
    # Constructor
    function MultiDiscreteMeasureData(parameter_ref::JuMPC.SparseAxisArray{<:ParameterRef},
                                      coeffs::Vector{<:Number},
                                      supports::Vector{<:JuMPC.SparseAxisArray},
                                      name::String, weight_func::Function)
        if length(coeffs) != length(supports)
            error("The amount of coefficients must match the amount of " *
                  "support points.")
        elseif !_only_one_group(parameter_ref)
            error("Parameters must all have same group ID.")
        elseif length(supports) == 0
            error("Must specify at least one support.")
        elseif keys(supports[1].data) != keys(parameter_ref.data)
            error("The keys/dimensions of the support points and parameters " *
                  "do not match.")
        end
        for key in keys(parameter_ref.data)
            set = infinite_set(parameter_ref.data[key])
            for supp in supports
                if !supports_in_set(supp.data[key], set)
                    error("Support points violate parameter domain.")
                end
            end
        end
        return new(parameter_ref, coeffs, supports, name, weight_func)
    end
end

# Define finite measure expressions (note infinite expression take precedence)
const MeasureExpr = Union{MeasureRef,
                          JuMP.GenericAffExpr{Float64, MeasureRef},
                          JuMP.GenericAffExpr{Float64, MeasureFiniteVariableRef},
                          JuMP.GenericQuadExpr{Float64, MeasureRef},
                          JuMP.GenericQuadExpr{Float64, MeasureFiniteVariableRef}}

"""
    BoundedScalarConstraint{F <: JuMP.AbstractJuMPScalar,
                            S <: MOI.AbstractScalarSet} <: JuMP.AbstractConstraint

A DataType that stores infinite constraints defined on a subset of the infinite
parameters on which they depend.

**Fields**
- `func::F` The JuMP object.
- `set::S` The MOI set.
- `bounds::ParameterBounds` Set of valid parameter sub-domains that further bound
                            constraint.
- `orig_bounds::ParameterBounds` Set of the constraint's original parameter
                                 sub-domains (not considering hold variables)
"""
struct BoundedScalarConstraint{F <: JuMP.AbstractJuMPScalar,
                               S <: MOI.AbstractScalarSet} <: JuMP.AbstractConstraint
    func::F
    set::S
    bounds::ParameterBounds
    orig_bounds::ParameterBounds
end

"""
    GeneralConstraintRef

An abstract type for constraint references unique to InfiniteOpt.
"""
abstract type GeneralConstraintRef end

"""
InfiniteConstraintRef{S <: JuMP.AbstractShape} <: GeneralConstraintRef

A DataType for constraints that contain infinite variables.

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of constraint in model.
- `shape::JuMP.AbstractShape` Shape of constraint
"""
struct InfiniteConstraintRef{S <: JuMP.AbstractShape} <: GeneralConstraintRef
    model::InfiniteModel
    index::Int
    shape::S
end

"""
    FiniteConstraintRef{S <: JuMP.AbstractShape} <: GeneralConstraintRef

A DataType for constraints that contain finite variables.

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of constraint in model.
- `shape::JuMP.AbstractShape` Shape of constraint
"""
struct FiniteConstraintRef{S <: JuMP.AbstractShape} <: GeneralConstraintRef
    model::InfiniteModel
    index::Int
    shape::S
end

"""
    MeasureConstraintRef{S <: JuMP.AbstractShape} <: GeneralConstraintRef

A DataType for constraints that contain finite variables and measures.

**Fields**
- `model::InfiniteModel` Infinite model.
- `index::Int` Index of constraint in model.
- `shape::JuMP.AbstractShape` Shape of constraint
"""
struct MeasureConstraintRef{S <: JuMP.AbstractShape} <: GeneralConstraintRef
    model::InfiniteModel
    index::Int
    shape::S
end


text:
module TestDay16

import AOC2021.Day16.part1, AOC2021.Day16.part2, AOC2021.Day16.bits, AOC2021.Day16.extract_next_packet
using Test

@testset "Day 16 basics" begin
    @testset "basic literal" begin
        literal = "D2FE28"
        literal_bitstring = bits(literal)
        packet, _ = extract_next_packet(literal_bitstring)
        @test packet.decoded == 2021
    end
    @testset "basic operator length bit type 0" begin
        operator = "38006F45291200"
        operator_bitstring = bits(operator)
        packet, _ = extract_next_packet(operator_bitstring)
        literals = [packet for packet in packet.subpackets]
        @test literals[1].decoded == 10
        @test literals[2].decoded == 20
    end
    @testset "basic operator length bit type 1" begin
        operator = "EE00D40C823060"
        operator_bitstring = bits(operator)
        packet, _ = extract_next_packet(operator_bitstring)
        literals = [packet for packet in packet.subpackets]
        @test literals[1].decoded == 1
        @test literals[2].decoded == 2
        @test literals[3].decoded == 3
    end
end

@testset "Day 16" begin
    @testset "part 1" begin
        @test part1("8A004A801A8002F478") == 16
        @test part1("620080001611562C8802118E34") == 12
        @test part1("C0015000016115A2E0802F182340") == 23
        @test part1("A0016C880162017C3686B18A3D4780") == 31
    end
    @testset "part 2" begin
        @test part2("C200B40A82") == 3
        @test part2("04005AC33890") == 54
        @test part2("880086C3E88112") == 7
        @test part2("CE00C43D881120") == 9
        @test part2("D8005AC2A8F0") == 1
        @test part2("F600BC2D8F") == 0
        @test part2("9C005AC2F8F0") == 0
        @test part2("9C0141080250320F1802104A08") == 1
    end
end

end # module


text:
type Knn
  k::Int
  x
  y
  means_train
  std_train
  norm_fx
  Knn() = new(3)
end

function train!{U<:Real,T}(model::Knn, x::Matrix{U}, y::Vector{T})
  model.x = x
  model.y = y
  nRows, nCols = size(x)
  model.means_train = Vector{Float64}(nCols)
  model.std_train = Vector{Float64}(nCols)
  for j in 1:nCols
      model.means_train[j] = mean(x[:,j])
      model.std_train[j] = std(x[:,j])
  end
end

function predict{U<:Real}(model::Knn,x::Matrix{U})
    nRows, nCols = size(x)
    x_norm = deepcopy(x)
    for j in 1:nCols
      for i in  1:nRows
         x_norm[i,j] = (x[i,j] - model.means_train[j])/(model.std_train[j])
      end
    end
    if (eltype(model.y)<:Real)
       #Regression
        resp = Vector{Float64}(nRows)
        for i in 1:size(x_norm,1)
          resp[i]= assign_value(x_norm, model.y, model.k, x_norm[i,:])
        end
        return resp
    else
      error("Classification is not implemented yet")
    end
end

function euclidean_distance(a, b)
 distance = 0.0
 for index in 1:size(a, 1)
  distance += (a[index]-b[index]) * (a[index]-b[index])
 end
 return distance
end

function get_k_nearest_neighbors(xTrain, rowTestItem, k)
   nRows, nCols = size(xTrain)
   distances = Array(Float32,nRows)
    for index in 1:nRows
     distances[index] = euclidean_distance(xTrain[index,:], rowTestItem)
    end
   sortedNeighbors = sortperm(distances)
   kNearestNeighbors = sortedNeighbors[1:k]
   return kNearestNeighbors
end

# assign_value(model.x_norm, model.y, model.k, x_norm[i,:])
function assign_value(xTrain, yTrain, k, imageI)
 kNearestNeighbors = get_k_nearest_neighbors(xTrain, imageI, k)

 soma::Float64 = 0
 highestCount = 0
 mostPopularLabel = 0
 for n in kNearestNeighbors
  soma = soma + yTrain[n]
 end
 return soma / k
end

function assign_label(xTrain, yTrain, k, imageI)
 kNearestNeighbors = get_k_nearest_neighbors(xTrain, imageI, k)
 counts = Dict{Int, Int}()
 highestCount = 0
 mostPopularLabel = 0
 for n in kNearestNeighbors
  labelOfN = yTrain[n]
  if !haskey(counts, labelOfN)
   counts[labelOfN] = 0
  end
  counts[labelOfN] += 1 #add one to the count
  if counts[labelOfN] > highestCount
   highestCount = counts[labelOfN]
   mostPopularLabel = labelOfN
  end
 end
 return mostPopularLabel
end


text:
n = 256;

data = rand(n,n);

sizeX = size(data,1);
sizeY = size(data,2);

# Set the variables for the Shearlet trasform
rows = sizeX;
cols = sizeY;
X = data;
nScales = 2;

shearletSystem = Shearlab.SLgetShearletSystem2D(rows,cols,nScales);

coeffs = Shearlab.SLsheardec2D(X,shearletSystem);

@test size(coeffs[:,:,1]) == size(data)


text:
"""
`module JAC.AngularCoefficientsRatip2013`
    ... a submodel of JAC that provides data types for angular (T and V) coefficients and methods to calculate the angular
        coefficients between pairs of relativistic configuration state functions via an external Fortran library.
"""
module AngularCoefficientsRatip2013

    using  ..Basics,  ..Defaults,  ..ManyElectron
    using  Libdl: dlopen, dlsym

    export AngularTcoeff, AngularVcoeff

    const ANCOPATH = joinpath(@__DIR__, "..", "deps", "bin", "libanco-ratip2013.so")

    """
    `AngularCoefficientsRatip2013.twice(j::AngularJ64)`  ... returns `2*j` for a given angular momentum quantum number `j`.
    """
    twice(j::AngularJ64) = Int(2*j.num//j.den)

    """
    `AngularCoefficientsRatip2013.getsubshells(csf::CsfR)`
        ... returns the array of subshells of a relativistic configuration state function.
    """
    getsubshells(csf::CsfR) = csf.useStandardSubshells ?   Defaults.GBL_STANDARD_SUBSHELL_LIST : csf.subshells

    """
    `AngularCoefficientsRatip2013.have_same_subshells(csfa::CsfR, csfb::CsfR)`
        ... checks whether two relativistic configurations state functions have the same subshells; true is returned if
            both CSF have the same subshell list, and false otherwise.
    """
    function have_same_subshells(csfa::CsfR, csfb::CsfR)
        return csfa.useStandardSubshells && csfb.useStandardSubshells || getsubshells(csfa)==getsubshells(csfb)
    end

    """
    `struct AngularCoefficientsRatip2013.Fnkappa`
        ... type that mirrors `nkappa` from the Fortran library. This type is used to communicate with the external Fortran
            library. It can be converted from and to a `Subshell` using the appropriate constructors.
    """
    struct Fnkappa      # nkappa
         n    ::Int32   # integer
         kappa::Int32   # integer
    end

    Fnkappa(s::Subshell) = Fnkappa(s.n, s.kappa)

    """
    `AngularCoefficientsRatip2013.Subshell(nk::Fnkappa)`
        ... constructor to define a Subshell (instance) from the quantum numbers n, kappa.
    """
    Subshell(nk::Fnkappa) = Subshell(nk.n, nk.kappa)

    """
    `AngularCoefficientsRatip2013.maxocc(nk::Fnkappa)`  ... return the maximal occupation number for a given relativistic subshell.
    """
    maxocc(nk::Fnkappa) = 2abs(nk.kappa)

    """
    `struct AngularCoefficientsRatip2013.FancoTcoeff`
        ... type that mirrors `anco_T_coeff` from the Fortran library. This type is used to communicate with the external
            Fortran library. It can be converted to a `AngularTcoeff` using the appropriate constructor.
    """
    struct FancoTcoeff # anco_T_coeff
        nu::Int32      # integer
        a::Fnkappa     # nkappa
        b::Fnkappa     # nkappa
        T::Float64     # real(kind=dp)
    end

    """
    `struct AngularCoefficientsRatip2013.FancoVcoeff`
        ... type that mirrors `anco_V_coeff` from the Fortran library. This type is used to communicate with the external
            Fortran library. It can be converted to an `AngularVcoeff` using the appropriate constructor.
    """
    struct FancoVcoeff # anco_V_coeff
        nu::Int32      # integer
        a::Fnkappa     # nkappa
        b::Fnkappa     # nkappa
        c::Fnkappa     # nkappa
        d::Fnkappa     # nkappa
        V::Float64     # real(kind=dp)
    end

    """
    `struct AngularCoefficientsRatip2013.Fmtcoefficient`
        ... type that mirrors `mct_coefficient` from the Fortran library. This type is used to communicate with the
            external Fortran library. It can be converted to an `AngularTcoeff` using the appropriate constructor.
    """
    struct Fmctcoefficient # mct_coefficient
        r::Int32           # integer
        s::Int32           # integer
        a::Int16           # integer(kind=i2b)
        b::Int16           # integer(kind=i2b)
        ν::Int16           # integer(kind=i2b)
        T::Float64         # real(kind=dp)
    end

    """
    `struct AngularTcoeff`  ... type for representing angular T coefficients.
    """
    struct AngularTcoeff
        nu::Int
        a::Subshell
        b::Subshell
        T::Float64
    end

    # `Base.show(io::IO, coeff::AngularTcoeff)`  ... prepares a proper printout of the variable coeff::AngularTcoeff.
    function Base.show(io::IO, coeff::AngularTcoeff)
        sa = " T^$(coeff.nu)[" * string(coeff.a) * "," * string(coeff.b) * "] = $(coeff.T)"
        print(io, sa )
    end

    """
    `AngularCoefficientsRatip2013.AngularTcoeff(t::FancoTcoeff)`
        ... constructor to define an AngularTcoeff (instance) from an internal FancoTcoeff coefficient.
    """
    AngularTcoeff(t::FancoTcoeff) = AngularTcoeff(t.nu, Subshell(t.a), Subshell(t.b), t.T)

    """
    `AngularCoefficientsRatip2013.AngularTcoeff(t::Fmctcoefficient)`
        ... constructor to define an AngularTcoeff (instance) from an internal Fmctcoefficient and the list of relativistic
            subshells the CSFs are based on.
    """
    AngularTcoeff(t::Fmctcoefficient, subshells::AbstractArray{Subshell}) = AngularTcoeff(t.ν, subshells[t.a], subshells[t.b], t.T)

    """
    `struct AngularCoefficientsRatip2013.AngularTcoeff`  ... type for representing angular V coefficients.
    """
    struct AngularVcoeff
        nu::Int
        a::Subshell
        b::Subshell
        c::Subshell
        d::Subshell
        V::Float64
    end

    # `Base.show(io::IO, coeff::AngularVcoeff)`  ... prepares a proper printout of the variable coeff::AngularVcoeff.
    function Base.show(io::IO, coeff::AngularVcoeff)
        sa = " V^$(coeff.nu)[" * string(coeff.a) * "," * string(coeff.b) * ";" * string(coeff.c) * "," * string(coeff.d) * "] = $(coeff.V)"
        print(io, sa )
    end

    """
    `AngularCoefficientsRatip2013.AngularVcoeff(t::FancoVcoeff)`
        ... constructor to define an AngularVcoeff (instance) from an internal FancoVcoeff coefficient.
    """
    AngularVcoeff(t::FancoVcoeff) = AngularVcoeff(t.nu, Subshell(t.a), Subshell(t.b), Subshell(t.c), Subshell(t.d), t.V)

    const NUMBER_OF_MCT_COEFFICIENTS_MAX = 1000

    # The following `Ref`s need to be initialized at runtime (via the __init__() function)
    # in order to enable precompilation

    const ANCOLIB  = Ref{Ptr{Nothing}}() # dlopen(ANCOPATH)

    const no_anco_t_list             = Ref{Ptr{Int32}}() # cglobal(dlsym(ANCOLIB[], :no_anco_t_list), Int32)
    const no_anco_v_list             = Ref{Ptr{Int32}}() # cglobal(dlsym(ANCOLIB[], :no_anco_v_list), Int32)
    const number_of_mct_coefficients = Ref{Ptr{Int32}}() # cglobal(dlsym(ANCOLIB[], :number_of_mct_coefficients), Int32)

    const anco_t_list = Ref{Ptr{FancoTcoeff}}()     # cglobal(dlsym(ANCOLIB[], :anco_t_list), FancoTcoeff)
    const anco_v_list = Ref{Ptr{FancoVcoeff}}()     # cglobal(dlsym(ANCOLIB[], :anco_v_list), FancoVcoeff)
    const mct_list    = Ref{Ptr{Fmctcoefficient}}() # cglobal(dlsym(ANCOLIB[], :mct_list), Fmctcoefficient)

    const t_vector   = Ref{Vector{FancoTcoeff}}()     # unsafe_wrap(Vector{FancoTcoeff}, anco_t_list, 1000)
    const v_vector   = Ref{Vector{FancoVcoeff}}()     # unsafe_wrap(Vector{FancoVcoeff}, anco_v_list, 2000)
    const mct_vector = Ref{Vector{Fmctcoefficient}}() # unsafe_wrap(Vector{Fmctcoefficient}, mct_list, NUMBER_OF_MCT_COEFFICIENTS_MAX)

    function __init__()
        ANCOLIB[] = dlopen(ANCOPATH)

        no_anco_t_list[]             = cglobal(dlsym(ANCOLIB[], :no_anco_t_list), Int32)
        no_anco_v_list[]             = cglobal(dlsym(ANCOLIB[], :no_anco_v_list), Int32)
        number_of_mct_coefficients[] = cglobal(dlsym(ANCOLIB[], :number_of_mct_coefficients), Int32)

        anco_t_list[] = cglobal(dlsym(ANCOLIB[], :anco_t_list), FancoTcoeff)
        anco_v_list[] = cglobal(dlsym(ANCOLIB[], :anco_v_list), FancoVcoeff)
        mct_list[]    = cglobal(dlsym(ANCOLIB[], :mct_list), Fmctcoefficient)

        t_vector[]   = unsafe_wrap(Vector{FancoTcoeff}, anco_t_list[], 1000)
        v_vector[]   = unsafe_wrap(Vector{FancoVcoeff}, anco_v_list[], 2000)
        mct_vector[] = unsafe_wrap(Vector{Fmctcoefficient}, mct_list[], NUMBER_OF_MCT_COEFFICIENTS_MAX)
    end

    """
    `AngularCoefficientsRatip2013.number_of_t_coeffs()`
        ... to read the number of calculated T coefficients from the external library.
    """
    number_of_t_coeffs() = unsafe_load(no_anco_t_list[])

    """
    `AngularCoefficientsRatip2013.number_of_v_coeffs()`
        ... to read the number of calculated V coefficients from the external library.
    """
    number_of_v_coeffs() = unsafe_load(no_anco_v_list[])

    """
    `AngularCoefficientsRatip2013.number_of_mct_coeffs()`
        ... to read the number of calculated MCT coefficients from the external library.
    """
    number_of_mct_coeffs() = unsafe_load(number_of_mct_coefficients[])

    """
    `AngularCoefficientsRatip2013.view_t_coeffs()`
        ... to return the array of T coefficients as a view into the memory of the external library.
    """
    view_t_coeffs() = view(t_vector[], 1:number_of_t_coeffs())

    """
    `AngularCoefficientsRatip2013.view_v_coeffs()`
        ... to return the array of V coefficients as a view into the memory of the external library.
    """
    view_v_coeffs() = view(v_vector[], 1:number_of_v_coeffs())

    """
    `AngularCoefficientsRatip2013.view_mct_coeffs()`
        ... to return the array of MCT coefficients as a view into the memory of the external library.
    """
    view_mct_coeffs() = view(mct_vector[], 1:number_of_mct_coeffs())

    """
    `AngularCoefficientsRatip2013.angular_coefficients_pair(r::Integer, s::Integer)`
        ... to calculate the angular coefficients between the `r`-th and `s`-th configuration state function of the
            CSF set that was loaded into the external library.
    """
    function angular_coefficients_pair(r::Integer, s::Integer)
        ccall(dlsym(ANCOLIB[], :jac_anco_calculate_csf_pair), Nothing, (Ref{Int32}, Ref{Int32}), r, s)
        return AngularTcoeff.(view_t_coeffs()), AngularVcoeff.(view_v_coeffs())
    end

    """
    `AngularCoefficientsRatip2013.angular_coefficients_pair_1p(ν::Integer, r::Integer, s::Integer)`
        ... to calculate the angular coefficients between the `r`-th and `s`-th configuration state function of the CSF
            set that was loaded into the external library, for a non-scalar one-particle operator of rank `ν`.
    """
    function angular_coefficients_pair_1p(ν::Integer, r::Integer, s::Integer)
        ccall(dlsym(ANCOLIB[], :jac_anco_calculate_csf_pair_1p), Nothing, (Ref{Int32}, Ref{Int32}, Ref{Int32}), ν, r, s)
        return AngularTcoeff.(view_t_coeffs())
    end

    """
    `AngularCoefficientsRatip2013.load_csl(nocsf, nwshells, nwcore, number_of_electrons, subshell,
                                           totalJs, parities, occupations, seniorities, subshellJs, subshellXs)`
        ... a wrapper function for the call to the subroutine `jac_anco_load_csl` of the external library.
    """
    load_csl(nocsf::Integer, nwshells::Integer, nwcore::Integer, number_of_electrons::Integer,
             subshell::AbstractVector{Fnkappa},
             totalJs::AbstractVector{Int8}, parities::AbstractVector{UInt8},
             occupations::AbstractMatrix{Int8}, seniorities::AbstractMatrix{Int8},
             subshellJs::AbstractMatrix{Int8},  subshellXs::AbstractMatrix{Int8}) =
        ccall(dlsym(ANCOLIB[], :jac_anco_load_csl), Nothing,
             (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Fnkappa}, Ref{Int8},
              Ref{UInt8}, Ref{Int8}, Ref{Int8}, Ref{Int8}, Ref{Int8}),
              nocsf, nwshells, nwcore, number_of_electrons, subshell, totalJs, parities,
              occupations, seniorities, subshellJs, subshellXs)

    """
    `AngularCoefficientsRatip2013.load_csl(csf::CsfR)`
        ... to load a single relativistic configuration state function into the external library. Returns the list of
            subshells the CSF is based on.
    """
    function load_csl(csf::CsfR)
        number_of_electrons = sum(csf.occupation)
        subshells = getsubshells(csf)
        nkappas = Fnkappa.(subshells)
        nwshells = length(nkappas)
        occ = reshape(Int8.(csf.occupation), nwshells, 1)
        sen = reshape(Int8.(csf.seniorityNr), nwshells, 1)
        suJ = reshape(Int8.(twice.(csf.subshellJ)), nwshells, 1)
        suX = reshape(Int8.(twice.(csf.subshellX)), nwshells, 1)
        totalJs =  Int8[twice(csf.J)]
        # parities = UInt8[Char(csf.parity)]
        if   csf.parity == Basics.plus    parities = UInt8[Char('+')]   else    parities = UInt8[Char('-')]   end
        nwcore = 0
        for i in eachindex(nkappas)
            if  occ[i]!=maxocc(nkappas[i])
                nwcore = i-1; break
            end
        end
        load_csl(1, nwshells, nwcore, number_of_electrons, nkappas, totalJs, parities, occ, sen, suJ, suX)
        subshells
    end

    """
    `AngularCoefficientsRatip2013.load_csl(csfa::CsfR, csfb::CsfR)`
        ... to load two relativistic configuration state functions into the external library. Returns a common
            list of subshells for both CSFs.
    """
    function load_csl(csfa::CsfR, csfb::CsfR)
        number_of_electrons = sum(csfa.occupation)
        sum(csfb.occupation) == number_of_electrons || error("CSFs must have the same number of electrons")
        if have_same_subshells(csfa, csfb)
            subshells = getsubshells(csfa)
            nwshells = length(subshells)
            occ = Matrix{Int8}(hcat(csfa.occupation, csfb.occupation))
            sen = Matrix{Int8}(hcat(csfa.seniorityNr, csfb.seniorityNr))
            suJ = Matrix{Int8}(twice.(hcat(csfa.subshellJ, csfb.subshellJ)))
            suX = Matrix{Int8}(twice.(hcat(csfa.subshellX, csfb.subshellX)))
        else
            subshells = sort(unique(vcat(getsubshells(csfa), getsubshells(csfb))))
            nwshells = length(subshells)
            occ = zeros(Int8, nwshells, 2)
            sen = zeros(Int8, nwshells, 2)
            suJ = zeros(Int8, nwshells, 2)
            suX = zeros(Int8, nwshells, 2)
            aindex, bindex = 1, 1
            ashells, bshells = getsubshells.((csfa, csfb))
            for i in eachindex(subshells)
                if ashells[aindex] == subshells[i]
                    occ[i,1] = csfa.occupation[aindex]
                    sen[i,1] = csfa.seniorityNr[aindex]
                    suJ[i,1] = twice(csfa.subshellJ[aindex])
                    suX[i,1] = twice(csfa.subshellX[aindex])
                    endof(ashells)<aindex && (aindex+=1)
                else
                    suX[i,1] = i==1 ? 0 : suX[i-1,1]
                end
                if bshells[bindex] == subshells[i]
                    occ[i,2] = csfb.occupation[bindex]
                    sen[i,2] = csfb.seniorityNr[bindex]
                    suJ[i,2] = twice(csfb.subshellJ[bindex])
                    suX[i,2] = twice(csfb.subshellX[bindex])
                    endof(bshells)<bindex && (bindex+=1)
                else
                    suX[i,2] = i==1 ? 0 : suX[i-1,2]
                end
            end
            aindex==endof(ashells) && bindex==endof(bshells) || error("Currently, only standard subshell ordering is supported")
        end
        totalJs =  Int8[twice(csfa.J), twice(csfb.J)]
        ## parities = UInt8[Char(csfa.parity), Char(csfb.parity)]
        parities = UInt8['a', 'a']
        if csfa.parity == Basics.plus      parities[1] = '+'     else       parities[1] = '-'      end
        if csfb.parity == Basics.plus      parities[2] = '+'     else       parities[2] = '-'      end
        nkappas = Fnkappa.(subshells)
        nwcore = 0
        for i in eachindex(subshells)
            coreocc = maxocc(nkappas[i])
            if  occ[i,1]!=coreocc || occ[i,2]!=coreocc
                nwcore = i-1; break
            end
        end
        load_csl(2, nwshells, nwcore, number_of_electrons, nkappas, totalJs, parities, occ, sen, suJ, suX)
        subshells
    end

    """
    `AngularCoefficientsRatip2013.mct_generate_coefficients(rl::Integer, ru::Integer, sl::Integer, su::Integer, iopar::Integer, rank::Integer)`
        ... to calculate the MCT coefficients between the `rl`..`ru`-th and `sl`..`su`-th configuration state functions of the CSF set
            that was loaded into the external library, for a non-scalar one-particle operator of parity `iopar` and rank `rank`.
    """
    function mct_generate_coefficients(rl::Integer, ru::Integer, sl::Integer, su::Integer, iopar::Integer, rank::Integer)
        ccall(dlsym(ANCOLIB[], :jac_mct_generate_coefficients), Nothing,
              (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}), rl, ru, sl, su, iopar, rank)
        return view_mct_coeffs()
    end

end # module


text:
"""
Julia program to implement Maximum Sum Decreasing Subsequence
In this problem, given an array we have to find the maximum sum an decreasing subsequence of that array can make.
This problem is a slight modification to the Longest Decreasing subsequence problem.
The problem can be solved using Dynamic Programming
"""

function maximum_sum_decreasing_subsequence(arr, n)
    max_sum = 0
    dp = Int[]

    # Initialize the dp array with the array values, as the maximum sum
    # at each point is atleast as the value at that point
    for i in 1:n
        push!(dp, arr[i])
    end

    for i in 1:n
        for j in 1:(i-1)
            if(arr[i] < arr[j] && dp[i] < dp[j] + arr[i])
                dp[i] = dp[j] + arr[i]
            end
        end
    end

    # Now Find the maximum element in the dp array
    max_sum = findmax(dp)[1]
    return max_sum

end

print("What is the length of the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("No numbers present in the array!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
res = maximum_sum_decreasing_subsequence(arr, n)
print("The maximum sum of an  decreasing subsequence of the given array is $res")

"""
Time Complexity: O(num ^ 2), where 'num' is the size of the given array
Space Complexity: O(num)

SAMPLE INPUT AND OUTPUT

SAMPLE 1

What is the length of the array? 10
Enter the numbers: 23 14 5 84 24 222 321 43 123 432
The maximum sum of an  decreasing subsequence of the given array is 444

SAMPLE 2

What is the length of the array? 5
Enter the numbers: 5 4 3 2 1
The maximum sum of an  decreasing subsequence of the given array is 15
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
module TestSpecies

  using GBIF
  using Test

  iver = taxon("Iris versicolor", rank=:SPECIES)
  @test iver.species == Pair("Iris versicolor", 5298019)

  iver_id = taxon(5298019)  # test using taxon id
  @test iver_id.species == iver.species

  @test iver_id == taxon(Pair("Iris versicolor", 5298019))

  i_sp = taxon(iver.genus.first; rank=:GENUS, family=iver.family.first, strict=false)
  @test ismissing(i_sp.species)

  @test typeof(occurrences(i_sp)) <: GBIFRecords

  iver_occ = occurrences(iver)
  @test typeof(iver_occ) <: GBIFRecords
  for iocc in iver_occ
    @test iocc.taxon.species == Pair("Iris versicolor", 5298019)
  end

  iver_occ_spain = occurrences(iver, "country" => "ES")
  @test typeof(iver_occ_spain) <: GBIFRecords
  for iocc in iver_occ_spain
    @test iocc.taxon.species == Pair("Iris versicolor", 5298019)
    @test iocc.country == "Spain"
  end

end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
        required = true
        arg_type = Int
        help = "seed"
    "dataset"
        required = true
        arg_type = String
        help = "dataset"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, contamination = parsed_args

#######################################################################################
################ THIS PART IS TO BE PROVIDED FOR EACH MODEL SEPARATELY ################
modelname = "if"
# sample parameters, should return a Dict of model kwargs

function set_params(data)
	return (n_estimators=100, max_samples="auto",max_features=1.0)
end
function fit(data, parameters)
	# construct model - constructor should only accept kwargs
	model = GenerativeAD.Models.IForest(;parameters...)

	# fit train data
	try
		global info, fit_t, _, _, _ = @timed fit!(model, data[1][1])
	catch e
		# return an empty array if fit fails so nothing is computed
		return (fit_t = NaN,), []
	end

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		model = nothing
		)

	# now return the different scoring functions
	training_info, [(x->predict(model, x), parameters)]
end

####################################################################
################ THIS PART IS COMMON FOR ALL MODELS ################
# set a maximum for parameter sampling retries
try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
while try_counter < max_tries
    for seed in 1:max_seed
		savepath = datadir("experiments/tabular_clean_val_default/$(modelname)/$(dataset)/seed=$(seed)")
		mkpath(savepath)

		# get data
		data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)

		# parameters
		parameters = set_params(data)

		@info "Trying to fit $modelname on $dataset with parameters $(parameters)..."
		# check if a combination of parameters and seed alread exists
		if GenerativeAD.check_params(savepath, parameters)
			# fit
			training_info, results = fit(data, parameters)
			# here define what additional info should be saved together with parameters, scores, labels and predict times
			save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, contamination = contamination))

			# now loop over all anomaly score funs
			for result in results
				GenerativeAD.experiment(result..., data, savepath; save_entries...)
			end
			global try_counter = max_tries + 1
		else
			@info "Model already present, trying new hyperparameters..."
			global try_counter += 1
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
@testset "CPReward" begin
    @testset "set_reward!(DecisionPhase)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        lh = SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, SeaPearl.CPReward, SeaPearl.FixedOutput}(agent)
        SeaPearl.update_with_cpmodel!(lh, model)

        lh.reward.value = 0
        SeaPearl.set_reward!(SeaPearl.DecisionPhase, lh, model)
        @test lh.reward.value == 0
    end

    @testset "set_reward!(StepPhase)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        lh = SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, SeaPearl.CPReward, SeaPearl.FixedOutput}(agent)
        SeaPearl.update_with_cpmodel!(lh, model)

        lh.reward.value = 0
        SeaPearl.set_reward!(SeaPearl.StepPhase, lh, model,nothing)
        @test lh.reward.value == 0
    end

    @testset "set_reward!(EndingPhase)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(3, 4, "y", trailer)
        z = SeaPearl.IntVar(4, 5, "z", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addVariable!(model, z; branchable=false)
        lh = SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, SeaPearl.CPReward, SeaPearl.FixedOutput}(agent)
        SeaPearl.update_with_cpmodel!(lh, model)

        lh.reward.value = 5
        model.statistics.numberOfNodes = 30
        SeaPearl.set_reward!(SeaPearl.EndingPhase, lh, model, :FoundSolution)
        @test lh.reward.value == 6
        SeaPearl.assign!(x,2)
        SeaPearl.set_reward!(SeaPearl.EndingPhase, lh, model, :Infeasible)
        @test lh.reward.value == 5
    end
end


text:
# Author: Dharmesh Tarapore, dharmesh@bu.edu
# Date: 12/10/2017

module GDQNACASSimulatorImpl

export
    addObserver,

    simulate,
    GDQNACASSimulator

using RiskRatio

using AbstractSimulatorImpl
using CommonInterfaces

using Base.Test

using Encounter
using PilotResponse
using DynamicModel
using WorldModel
using Sensor
using CollisionAvoidanceSystem
using RLESUtils, Observers

import CommonInterfaces.addObserver

type SimulationParameters

    em::AbstractEncounterModel
    pr::Vector{AbstractPilotResponse}
    dm::Vector{AbstractDynamicModel}
    wm::AbstractWorldModel
    sr::Vector{AbstractSensor}
    cas::Vector{AbstractCollisionAvoidanceSystem}

    number_of_aircraft::Int

    SimulationParameters() = new()
end

type GDQNACASSimulator <: AbstractSimulator

    parameters::SimulationParameters

    time_step::Float64

    observer::Observer

    function GDQNACASSimulator()

        obj = new()

        obj.parameters = SimulationParameters()
        obj.time_step = 1.

        obj.observer = Observer()

        return obj
    end
end

addObserver(sim::GDQNACASSimulator, f::Function) = add_observer(sim.observer, f)
addObserver(sim::GDQNACASSimulator, tag::AbstractString, f::Function) = add_observer(sim.observer, tag, f)

import Base.convert

convert(::Type{SimpleADMInitialState}, state::Union{CorrAEMInitialState, LLAEMInitialState}) = SimpleADMInitialState(state.t, state.x, state.y, state.h, state.v, state.psi, state.h_d)

convert(::Type{ASWMState}, state::SimpleADMOutputState) = ASWMState(state.t, state.x, state.y, state.h, state.vx, state.vy, state.vh)

convert(::Type{SimplePRCommand}, command::Union{CorrAEMCommand, LLAEMCommand}) = SimplePRCommand(command.t, command.v_d, command.h_d, command.psi_d, 0.0)

convert(::Type{SimpleADMCommand}, command::SimplePRCommand) = SimpleADMCommand(command.t, command.v_d, command.h_d, command.psi_d)

function convert(::Type{SimpleGDQNACASSensorInput}, states::Vector{ASWMState})

    states_ = Array(SimpleGDQNACASSRState, length(states))

    for i = 1:length(states)
        states_[i] = GDQNACASState(states[i].t, states[i].x, states[i].y, states[i].h, states[i].vx, states[i].vy, states[i].vh)
    end

    return GDQNACASSensorInput(states_)
end

convert(::Type{GDQNACASInput}, output::GDQNACASSensorOutput) = GDQNCASInput(output.t, output.r, output.r_d, output.a, output.a_d, output.h, output.h_d)

convert(::Type{SimplePRResolutionAdvisory}, RA::GDQNACASResolutionAdvisory) = SimplePRResolutionAdvisory(RA.h_d)
convert(::Type{SimplePRResolutionAdvisory}, RA::Void) = nothing

function simulate(sim::AbstractSimulator; bACAS = false, sample_number = 0)

    aem = sim.parameters.em
    pr = sim.parameters.pr
    adm = sim.parameters.dm
    as = sim.parameters.wm
    sr = sim.parameters.sr
    cas = sim.parameters.cas

    number_of_aircraft = sim.parameters.number_of_aircraft
    @test number_of_aircraft == 2

    if sample_number == 0
        generateEncounter(aem)
    else
        generateEncounter(aem, sample_number = sample_number)
    end

    initial_coords = zeros(number_of_aircraft, 3)

    for i = 1:number_of_aircraft
        initial = Encounter.getInitialState(aem, i)
        initial_coords[i, :] = [initial.x, initial.y, initial.h]
    end

    for i = 1:number_of_aircraft
        initial = Encounter.getInitialState(aem, i)
        state = DynamicModel.initialize(adm[i], convert(SimpleADMInitialState, initial))
        WorldModel.initialize(as, i, convert(ASWMState, state))

        Sensor.initialize(sr[i])
        CollisionAvoidanceSystem.initialize(cas[i])
        PilotResponse.initialize(pr[i])
    end

    command = nothing

    while true
        for i = 1:number_of_aircraft
            command = Encounter.update(aem, i)

            if command == nothing
                break
            end

            if bACAS
                states = WorldModel.getAll(as)
                output = Sensor.update(sr[i], convert(GDQNACASSensorInput, states))
                RA = CollisionAvoidanceSystem.update(cas[i], convert(GDQNACASInput, output))

                if RA != nothing
                    @notify_observer(sim.observer, "RA", [i, states[i].t, RA.h_d])
                end
            else
                RA = nothing
            end

            response = PilotResponse.update(pr[i], convert(SimplePRCommand, command), convert(SimplePRResolutionAdvisory, RA))
            state = DynamicModel.update(adm[i], convert(SimpleADMCommand, response))
            WorldModel.update(as, i, convert(ASWMState, state))
        end

        if command == nothing
            break
        end

        WorldModel.updateAll(as)
    end

    return 0
end

# Used for quick estimates.
function GDQNACAS_thresholds(h)
    if h < 1000
        dmod = 0
        zthr = 0
    elseif h < 2350
        dmod = 0.2 * 6076.12
        zthr = 600
    elseif h < 5000
        dmod = 0.35 * 6076.12
        zthr = 600
    elseif h < 10000
        dmod = 0.55 * 6076.12
        zthr = 600
    elseif h < 20000
        dmod = 0.80 * 6076.12
        zthr = 600
    elseif h < 42000
        dmod = 1.10 * 6076.12
        zthr = 700
    else
        dmod = 1.10 * 6076.12
        zthr = 800
    end

    return dmod, zthr
end

end


text:
"""
    sigma2(pk, R)

Variance of mass fluctuation, σ^2, with a top-hat filter with radius `R`

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - **pk times k^3 must be dimensionless**. For example, if k is in units of h/Mpc, `pk` must be in units of Mpc^3/h^3.
- `R::Real`: a top-hat smoothing scale.
   - **R times k must be dimensionless**. For example, if k is in units of h/Mpc, `R` must be in units of Mpc/h.

`sigma2` is computed as

``σ^2(R) = ∫_0^∞ k^2dk W^2(kR) pk(k) / (2π^2)``

where and W(kR) is a window function given by

``W(x) = (3/x)j_1(x) = (3/x)(sin(x)/x^2 - cos(x)/x)``

j_1(x) is the spherical Bessel function of order 1.

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function sigma2(pk, R::Real)
   function dσ2dk(k)
      x = k * R
      W = (3 / x) * (sin(x) / x^2 - cos(x) / x)
      dσ2dk = W^2 * pk(k) * k^2 / 2 / π^2
   end
   res, err = quadde(dσ2dk, 0, 20 / R)
   σ2 = res
end

"""
    dsigma2dR(pk, R)

Derivative of variance of mass fluctuation, dσ^2/dR,  with respect to a top-hat filter with radius `R`.

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - **pk times k^3 must be dimensionless**. For example, if k is in units of h/Mpc, `pk` must be in units of Mpc^3/h^3.
- `R::Real`: a top-hat smoothing scale.
   - **R times k must be dimensionless**. For example, if k is in units of h/Mpc, `R` must be in units of Mpc/h.

`dsigma2/dR` is computed as

``dσ^2/dR = 2 ∫_0^∞ k^2dk W(kR) dW/dR pk(k) / (2π^2)``

where and W(kR) is a window function given by

``W(x) = (3/x)j_1(x) = (3/x)(sin(x)/x^2 - cos(x)/x)``

and dW/dR is its derivative:

``dW/dR = k dW/dx = (-3k/x)j_2(x) = (-3k/x)[(3/x^2-1)sin(x)/x - 3cos(x)/x^2]``

j_n(x) is the spherical Bessel function of order n.

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function dsigma2dR(pk, R::Real)
   function dσ2dRdk(k)
      x = k * R
      W = (3 / x) * (sin(x) / x^2 - cos(x) / x)
      dWdx = (-3 / x) * ((3 / x^2 - 1) * sin(x) / x - 3 * cos(x) / x^2)
      dσ2dRdk = 2 * W * dWdx * pk(k) * k^3 / 2 / π^2
   end
   res, err = quadde(dσ2dRdk, 0, 20 / R)
   dσ2dR = res
end

"""
    sigma2gaus(pk, R)

Variance of mass fluctuation, σ^2(`R`), with a Gaussian filter with width `R`.

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - **pk times k^3 must be dimensionless**. For example, if k is in units of h/Mpc, `pk` must be in units of Mpc^3/h^3.
- `R::Real`: a top-hat smoothing scale.
   - **R times k must be dimensionless**. For example, if k is in units of h/Mpc, `R` must be in units of Mpc/h.

`sigma2gaus(R)` is computed as

``σ^2(R) = ∫_0^∞ k^2dk exp(-k^2 R^2) pk(k) / (2π^2)``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function sigma2gaus(pk, R::Real)
   function dσ2dk(k)
      x = k * R
      dσ2dk = exp(-x^2) * pk(k) * k^2 / 2 / π^2
   end
   res, err = quadde(dσ2dk, 0, 10 / R)
   σ2 = res
end

"""
    dsigma2gausdR(pk, R)

Derivative of variance of mass fluctuation, dσ^2/dR, with respect to a Gaussian filter with width `R`.

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - **pk times k^3 must be dimensionless**. For example, if k is in units of h/Mpc, `pk` must be in units of Mpc^3/h^3.
- `R::Real`: a top-hat smoothing scale.
   - **R times k must be dimensionless**. For example, if k is in units of h/Mpc, `R` must be in units of Mpc/h.

`dsigma2gaus/dR` is computed as

``dσ^2/dR = -2R ∫_0^∞ k^4dk exp(-k^2 R^2) pk(k) / (2π^2)``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function dsigma2gausdR(pk, R::Real)
   function dσ2dRdk(k)
      x = k * R
      dσ2dRdk = -2R * exp(-x^2) * pk(k) * k^4 / 2 / π^2
   end
   res, err = quadde(dσ2dRdk, 0, 10 / R)
   dσ2dR = res
end

"""
    d2sigma2gausdR2(pk, R)

Second derivative of variance of mass fluctuation, d^2σ^2/dR^2, with respect to a Gaussian filter with width `R`.

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - **pk times k^3 must be dimensionless**. For example, if k is in units of h/Mpc, `pk` must be in units of Mpc^3/h^3.
- `R::Real`: a top-hat smoothing scale.
   - **R times k must be dimensionless**. For example, if k is in units of h/Mpc, `R` must be in units of Mpc/h.

`d2sigma2gaus/dR2` is computed as

``d^2σ^2/dR^2 = -2 ∫_0^∞ (1 - 2x^2) k^4dk exp(-k^2 R^2) pk(k) / (2π^2)``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function d2sigma2gausdR2(pk, R::Real)
   function d2σ2dR2dk(k)
      x = k * R
      d2σ2dR2dk = (-2 + 4 * x^2) * exp(-x^2) * pk(k) * k^4 / 2 / π^2
   end
   res, err = quadde(d2σ2dR2dk, 0, 10 / R)
   d2σ2dR2 = res
end


text:
using PolarFact
using Base.Test

include("test_newton.jl")
include("test_halley.jl")
include("test_svd.jl")
include("test_hybrid.jl")
include("test_f32.jl") # test Float32 case

println("Success!")