text:
for N in [Float64, Rational{Int}, Float32]
    # Sum of 2D centered balls in norm 1 and infinity
    b1 = BallInf(N[0, 0], N(2))
    b2 = Ball1(N[0, 0], N(1))
    # Test Construction
    X = MinkowskiSum(b1, b2)
    @test X.X == b1
    @test X.Y == b2

    # swap
    ms2 = swap(X)
    @test X.X == ms2.Y && X.Y == ms2.X

    # Test Dimension
    @test dim(X) == 2
    # Test Support Vector
    d = N[1, 0]
    v = σ(d, X)
    @test v[1] == N(3)
    d = N[-1, 0]
    v = σ(d, X)
    @test v[1] == N(-3)
    d = N[0, 1]
    v = σ(d, X)
    @test v[2] == N(3)
    d = N[0, -1]
    v = σ(d, X)
    @test v[2] == N(-3)

    # Sum of not-centered 2D balls in norm 1 and infinity
    b1 = BallInf(N[-1, 3], N(2))
    b2 = Ball1(N[1, 2], N(1))
    s = b1 + b2
    # Test Support Vector
    d = N[1, 0]
    v = σ(d, s)
    @test v[1] == N(3)
    d = N[-1, 0]
    v = σ(d, s)
    @test v[1] == N(-3)
    d = N[0, 1]
    v = σ(d, s)
    @test v[2] == N(8)
    d = N[0, -1]
    v = σ(d, s)
    @test v[2] == N(2)

    # Sum of array of LazySet
    # 2-elements
    ms = MinkowskiSum(Singleton(N[1]), Singleton(N[2]))
    @test ρ(N[1], ms) == 3
    @test ρ(N[-1], ms) == -3
    # 3-elements
    ms = MinkowskiSum(Singleton(N[1]), MinkowskiSum(Singleton(N[2]), Singleton(N[3])))
    @test ρ(N[1], ms) == N(6)
    @test ρ(N[-1], ms) == N(-6)

    # boundedness
    @test isbounded(ms)
    @test !isbounded(Singleton(N[1]) + HalfSpace(N[1], N(1)))

    # isempty
    @test !isempty(ms)

    # an_element function (falls back to the default implementation
    X = MinkowskiSum(LineSegment(N[0, 0], N[1, 1]),
                     LineSegment(N[1, 0], N[0, 1]))
    v = an_element(X)
    @test v ∈ Ball1(N[1, 1], N(1))

    # constraints list
    ms = MinkowskiSum(BallInf(N[0, 0], N(1)), BallInf(N[0, 0], N(1)))
    clist = constraints_list(ms)
    P = HPolytope(clist)
    H = Hyperrectangle(N[0, 0], N[2, 2])
    @test length(clist) == 4 && P ⊆ H && H ⊆ P

    # =================
    # MinkowskiSumArray
    # =================

    # relation to base type (internal helper functions)
    @test LazySets.array_constructor(MinkowskiSum) == MinkowskiSumArray
    @test LazySets.is_array_constructor(MinkowskiSumArray)

    # array getter
    v = Vector{LazySet{N}}()
    @test array(MinkowskiSumArray(v)) ≡ v

    # constructor with size hint and type
    MinkowskiSumArray(10, N)

    # in-place modification
    msa = MinkowskiSumArray(LazySet{N}[])
    @test MinkowskiSum!(b1, b1) isa MinkowskiSum && length(array(msa)) == 0
    res = MinkowskiSum!(b1, msa)
    @test res isa MinkowskiSumArray && length(array(msa)) == 1
    res = MinkowskiSum!(msa, b1)
    @test res isa MinkowskiSumArray && length(array(msa)) == 2
    res = MinkowskiSum!(msa, msa)
    @test res isa MinkowskiSumArray && length(array(msa)) == 4

    ms = MinkowskiSum(b1, b2)
    msa = MinkowskiSumArray([b1, b2])

    # dimension
    @test dim(msa) == 2

    # support vector
    d = N[1, 1]
    @test σ(d, ms) == σ(d, msa)
    d = N[-1, 1]
    @test σ(d, ms) == σ(d, msa)

    # boundedness
    @test isbounded(msa)
    @test !isbounded(MinkowskiSumArray([Singleton(N[1]), HalfSpace(N[1], N(1))]))

    # isempty
    @test !isempty(msa)

    # =================
    # CacheMinkowskiSum
    # =================

    # caching Minkowski sum
    cms = CacheMinkowskiSum(2, N)
    x1 = BallInf(ones(N, 3), N(3))
    x2 = Ball1(ones(N, 3), N(5))
    d = ones(N, 3)
    a = array(cms)
    push!(a, x1)
    σ(d, cms)
    push!(a, x2)
    svec = σ(d, cms)
    @test σ(d, cms) == svec
    @test dim(cms) == 3
    idx = forget_sets!(cms)
    @test idx == 2 && isempty(array(cms)) && cms.cache[d][1] == 0
    push!(a, x1)
    σ(d, cms)
    push!(a, x2)
    idx = forget_sets!(cms)
    @test idx == 1 && length(array(cms)) == 1 && cms.cache[d][1] == 0
    # test issue #367: inplace modification of the direction modified the cache
    d[1] = zero(N)
    @test haskey(cms.cache, ones(N, 3))
    # getindex
    cp = LazySets.CachedPair(1, N[2])
    @test cp[1] == 1 && cp[2] == N[2]
    @test_throws ErrorException cp[3]

    # boundedness
    @test isbounded(cms)
    @test !isbounded(CacheMinkowskiSum([Singleton(N[1]), HalfSpace(N[1], N(1))]))

    # isempty
    @test !isempty(cms)

    # ================
    # common functions
    # ================

    # neutral element
    z = ZeroSet{N}(2)
    msa = MinkowskiSumArray(LazySet{N}[])
    cms = CacheMinkowskiSum(LazySet{N}[])
    @test neutral(MinkowskiSum) == neutral(MinkowskiSumArray) ==
          neutral(CacheMinkowskiSum) == ZeroSet
    @test b1 + z == z + b1 == b1
    @test msa + z == z + msa == msa
    @test cms + z == z + cms == cms
    # absorbing element
    e = EmptySet{N}()
    @test absorbing(MinkowskiSum) == absorbing(MinkowskiSumArray) ==
          absorbing(CacheMinkowskiSum) == EmptySet
    @test b1 + e == e + b1 == msa + e == e + msa == cms + e == e + cms ==
          e + e == e
    # mix of neutral and absorbing element
    @test z + e == e + z == e
end


text:
import CUDArt
using Base.Test
using Compat

#########################
# Test CUDA array types #
#########################
result = CUDArt.devices(dev->CUDArt.capability(dev)[1] >= 2, nmax=1) do devlist
    # Copying memory to and from device
    @test length(devlist) == 1
    CUDArt.device(devlist[1])
    # Test CudaArrays and CudaPitchedArrays
    for AT in (CUDArt.CudaArray, CUDArt.CudaPitchedArray)
        g = AT(Float64, (5,3))
        if AT == CUDArt.CudaPitchedArray
            p = pointer(g)
            @test p.xsize == 5
            @test p.ysize == 3
        end
        h = CUDArt.to_host(g)
        @test typeof(h) == Array{Float64,2}
        @test size(h) == (5,3)
        g2 = AT(1:5)
        g3 = similar(g2)
        @test eltype(g3) == eltype(g2)
        @test size(g3) == size(g2)
        g3 = similar(g2, Float32)
        @test eltype(g3) == Float32
        g3 = similar(g2, (7,))
        @test eltype(g3) == Int
        @test size(g3) == (7,)
        @test typeof(g3) == AT{Int, 1}
        h2 = CUDArt.to_host(g2)
        @test typeof(h2) == Array{Int,1}
        @test h2 == [1:5;]
        if AT == CUDArt.CudaPitchedArray
            p = pointer(g2)
            @test p.xsize == 5
            @test p.ysize == 1
        end
        # fill!
        fill!(g2, uint8(0))
        h2 = CUDArt.to_host(g2)
        @test h2 == zeros(Int, 5)
        fill!(g2, 0)
        h2 = CUDArt.to_host(g2)
        @test h2 == zeros(Int, 5)
        fill!(g, 0.0)
        h = CUDArt.to_host(g)
        @test all(h .== 0)
        fill!(g, 0)
        h = CUDArt.to_host(g)
        @test all(h .== 0)
        fill!(g, 8.4)
        h = CUDArt.to_host(g)
        @test all(h .== 8.4)
        fill!(g, NaN)
        h = CUDArt.to_host(g)
        @test isequal(h, fill(NaN, size(g)))
        gbig = AT(Int32, (2000,1000))
        fill!(gbig, -17)
        hbig = CUDArt.to_host(gbig)
        @test all(hbig .== -17) && size(hbig) == (2000,1000)
    end
    # Getting portions of an array
    h_src = reshape(1.0:15.0, 3, 5)
    d_src = CUDArt.CudaPitchedArray(h_src)
    d_dest = CUDArt.CudaPitchedArray(Float64, (1, 2))
    get!(d_dest, d_src, (2,2:3), NaN)
    h_dest = CUDArt.to_host(d_dest)
    @test h_dest == [5 8]
    get!(d_dest, d_src, (2,0:1), NaN)
    h_dest = CUDArt.to_host(d_dest)
    @test isequal(h_dest, [NaN 2])
    get!(d_dest, d_src, (3,5:6), 7.3)
    h_dest = CUDArt.to_host(d_dest)
    @test isequal(h_dest, [15 7.3])
    # Copies between CudaPitchedArrays and CudaArrays
    A = rand(uint16(5:11), 7, 2)
    GA = CUDArt.CudaArray(A)
    GP = CUDArt.CudaPitchedArray(eltype(A), size(A))
    copy!(GP, GA)
    H = CUDArt.to_host(GP)
    @test H == A
    A[3,2] = 0
    copy!(GP, A)
    copy!(GA, GP)
    H = CUDArt.to_host(GA)
    @test H == A
    # HostArray (pinned memory)
    A = CUDArt.HostArray(Int32, (6,4))
    A[2,2] = 7
    @test A[2,2] == 7
    G = CUDArt.CudaArray(A)
    CUDArt.device_synchronize()
    A[2,2] = 3
    @test A[2,2] == 3
    copy!(A, G)
#     @test A[2,2] == 3  # asynchronous
    CUDArt.device_synchronize()
    @test A[2,2] == 7
    S = pointer_to_array(pointer(A)+48, (6,2), false)  # a "SubArray"
    B = rand(int32(1:15), 6, 2)
    GB = CUDArt.CudaArray(B)
    copy!(S, GB)
    CUDArt.device_synchronize()
    @test A[:,3:4] == B
    GB = CUDArt.CudaPitchedArray(2B)
    copy!(S, GB)
    CUDArt.device_synchronize()
    @test A[:,3:4] == 2B
end

#####################
# Executing kernels #
#####################
result = CUDArt.devices(dev->CUDArt.capability(dev)[1] >= 2, nmax=1) do devlist
    CUDArt.device(devlist[1])
    a = rand(Float32, 3, 4)
    b = rand(Float32, 3, 4)
    c = CUDArt.CuModule("vadd.ptx") do md
        vadd = CUDArt.CuFunction(md, "vadd")
        ga = CUDArt.CudaArray(a)
        gb = CUDArt.CudaArray(b)
        gc = CUDArt.CudaArray(Float32, size(a))
        CUDArt.launch(vadd, length(a), 1, (ga, gb, gc))
        CUDArt.to_host(gc)
    end
    @test_approx_eq c (a+b)
end

gc()  # check for finalizer errors

#########################################
# Multiple devices, streams, and wait() #
#########################################
if CUDArt.devcount() > 1
    CUDArt.devices(dev->true, nmax=2) do devlist
        sleeptime = 0.5
        results = Array(Any, ceil(Int, 2.5*length(devlist)))
        streams = [(CUDArt.device(dev); CUDArt.Stream()) for dev in devlist]
        # Force one run to precompile
        CUDArt.cudasleep(sleeptime; dev=devlist[1], stream=streams[1])
        wait(streams[1])
        i = 1
        nextidx() = (idx=i; i+=1; idx)
        @sync begin
            for idev = 1:length(devlist)
                sleep(0.01)
                @async begin
                    while true
                        idx = nextidx()
                        if idx > length(results)
                            break
                        end
                        tstart = time()
                        dev = devlist[idev]
                        stream = streams[idev]
                        CUDArt.cudasleep(sleeptime; dev=dev, stream=stream)
                        wait(stream)
                        tstop = time()
                        results[idx] = (tstart, tstop, dev)
                    end
                end
            end
        end
        for i = 1:length(results)
            r = results[i]
            @test sleeptime <= (r[2]-r[1]) <= 1.1*sleeptime
            @test r[3] == devlist[mod1(i, 2)]
        end
        @test results[end][2]-results[1][1] < 1.1*sleeptime*ceil(Int, length(results)/length(devlist))
        nothing
    end
end

ext: jl


text:
#SIR-X Model

using DifferentialEquations

beta = 0.215; gamma = 0.07; N = 7000;
eta0 = 0.003; eta = 0.003

"""SIR-X Model
Extension of the basic SIR model

"""
function SIR_X!(du,u,p,t)
    du[1] = -beta* u[1]/N*u[2] - eta0*u[1]                          #S
    du[2] = beta*u[1]/N*u[2] - gamma*u[2] - eta0*u[2] - eta*u[2]    #I
    du[3] = gamma*u[2] + eta0*u[1]                                  #R
    du[4] = (eta0 + eta)*u[2]                                       #X
end

   u0 = [N-15.0;15.0;0.0;0.0]
   tspan = (0.0,100.0)
   prob = ODEProblem(SIR_X!,u0,tspan)
   sol = solve(prob)

   plot(sol)


text:
# Automatically generated test file
function super_exponential3()

g = [
1.000000000000000,
1.064494458917860,
1.284025416687741,
1.755054656960298,
2.718281828459045,
4.770733181967603,
9.487735836358526,
21.380942759123343,
54.598150033144236,
-54.598150033144236,
-21.380942759123343,
-9.487735836358526,
-4.770733181967603,
-2.718281828459045,
-1.755054656960298,
-1.284025416687741]

ζ_vec = [6649.18, 1320.94, 262.421, 52.1332, 10.3569, 2.05753, 0.408753, 0.0812037, 0.0161321, 0.00320484]

I_sols = [
[1, 16],
[1, 8, 16],
[1, 9, 10, 16],
[1, 7, 9, 10, 16],
[1, 7, 9, 10, 11, 16],
[1, 6, 8, 9, 10, 11, 16],
[1, 6, 8, 9, 10, 11, 12, 16],
[1, 5, 7, 8, 9, 10, 11, 12, 16],
[1, 5, 7, 8, 9, 10, 11, 12, 13, 16],
[1, 4, 6, 7, 8, 9, 10, 11, 12, 13, 16],
[1, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16],
[1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16],
[1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
[1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]]

f_sols = [
6906.314103950247045,
6649.182023111979106,
1780.147768604715566,
810.634794153323128,
149.682164816596014,
79.779903677661423,
16.245588913173378,
8.328684761258955,
1.929086519737211,
0.854273036134146,
0.251603511120265,
0.081666565196429,
0.041290097946330,
0.004006052640761,
-0.000000000020918]

return g, ζ_vec, I_sols, f_sols

end


text:
using Parameters

"""
Contains a distribution over a triangulated state space. `dist::Vector{Float64}`
is the distribution, and `nonzero_inds::Vector{Int}` are the indices of the
simplices with non-zero measure.
"""
@with_kw struct InvDist
    dist::Vector{Float64} = Vector{Float64}() # Distribution over the simplices
    nonzero_inds::Vector{Int} = Vector{Int}() # indices of nonzero entries
end

"""
Compute the invariant probability distribution from a square Markov matrix `M`.
This is done by repeated application of `M` on an initially random distribution
until the distribution converges.
"""
function estimate_invdist(
        M::AbstractArray{Float64, 2};
        N::Int = 100,
        tolerance::Float64 = 1/10^5,
        delta::Float64 = 1/10^5
        )

    #=
    # Start with a random distribution `Ρ` (big rho). Normalise it so that it
    # sums to 1 and forms a true probability distribution over the simplices.
    =#
    Ρ = rand(Float64, 1, size(M, 1))
    Ρ = Ρ ./ sum(Ρ, 2)

    #=
    # Start estimating the invariant distribution. We could either do this by
    # finding the left-eigenvector of M, or by repeated application of M on Ρ
    # until the distribution converges. Here, we use the latter approach,
    # meaning that we iterate until Ρ doesn't change substantially between
    # iterations.
    =#
    distribution = Ρ * M

    distance = norm(distribution - Ρ) / norm(Ρ)

    check = floor(Int, 1 / delta)
    check_pts = floor.(Int, collect(1:N).' ./ check) .* collect(1:N).'
    check_pts = check_pts[check_pts .> 0]
    num_checkpts = size(check_pts, 1)
    check_pts_counter = 1

    counter = 1
    while counter <= N && distance >= tolerance
        counter += 1
        Ρ = distribution

        # Apply the Markov matrix to the current state of the distribution
        distribution = Ρ * M

        if (check_pts_counter <= num_checkpts &&
           counter == check_pts[check_pts_counter])

            check_pts_counter += 1
            colsum_distribution = sum(distribution, 2)[1]
            if abs(colsum_distribution - 1) > delta
                distribution = distribution ./ colsum_distribution
            end
        end

        distance = norm(distribution - Ρ) / norm(Ρ)
    end

    # Do the last normalisation and check
    colsum_distribution = sum(distribution, 2)[1]

    if abs(colsum_distribution - 1) > delta
        distribution = distribution ./ colsum_distribution
    end
    # Find simplices with strictly positive measure.
    simplex_inds_nonzero = heaviside(distribution) .* collect(1:size(M, 1)).'
    simplex_inds_nonzero = round(Int, simplex_inds_nonzero)
    simplex_inds_nonzero = simplex_inds_nonzero[simplex_inds_nonzero .> 0]

    # Extract the elements of the invariant measure corresponding to these indices
    return InvDist(dist = vec(distribution), nonzero_inds = simplex_inds_nonzero)
end


text:
using JuMP
using GLPK
using ParameterJuMP
using TimerOutputs

#=
    add parameters
=#
function add_parameters2(model, N)
    x = add_parameters(model, 4.0*ones(N))
end
function add_parameter_loop(model, N)
    x = [add_parameter(model, 4.0) for i in 1:N]
end

#=
    add variables
=#
function add_variables(model, N)
    @variable(model, y[1:N])
    return y
end

#=
    add constraints
=#
function add_constraints_0(model, x, y, M)
    N_x = length(x)
    N_y = length(y)
    @constraint(model, [ctr in 1:M],
        4 * sum(3.0*x[i] for i in 1:N_x) >= 2.0)
    return nothing
end

function add_constraints_1(model, x, y, M)
    N_x = length(x)
    N_y = length(y)
    @constraint(model, [ctr in 1:M],
        sum(3.0*y[i] for i in 1:N_y) >= 2.0 - sum(7.0*x[i] for i in 1:N_x))
    return nothing
end

function add_constraints_2(model, x, y, M)
    N_x = length(x)
    N_y = length(y)
    @assert N_x == N_y
    @constraint(model, [ctr in 1:M],
        sum(3.0*y[i] - 7.0*x[i] for i in 1:N_y) >= 2.0)
    return nothing
end

function add_constraints_3(model, x, y, M)
    N_x = length(x)
    N_y = length(y)
    @assert N_x == N_y
    @constraint(model, [ctr in 1:M],
        sum(3.0*y[i] + 8.0 - 7.0*x[i] for i in 1:N_y) >= 2.0)
    return nothing
end

#=
    query duals
=#
function get_duals(model, x)
    return dual.(x)
end

#=
    Complete tests
=#

function bench_create_param(N_Parameters::Int, N_Variables::Int, N_Constraints::Int)
    model = ModelWithParams(with_optimizer(GLPK.Optimizer))
    @timeit to "create params" x = add_parameters2(model, N_Parameters)
    @timeit to "create vars" y = add_variables(model, N_Variables)

    @timeit to "ctr vars" add_constraints_0(model, y, y, N_Constraints)
    @timeit to "ctr params" add_constraints_0(model, x, x, N_Constraints)
    @timeit to "ctr 2 bl vars" add_constraints_1(model, y, y, N_Constraints)
    @timeit to "ctr 2 bl params" add_constraints_1(model, x, x, N_Constraints)
    @timeit to "ctr 2 bl both" add_constraints_1(model, x, y, N_Constraints)
    @timeit to "ctr mixed vars" add_constraints_2(model, y, y, N_Constraints)
    @timeit to "ctr mixed params" add_constraints_2(model, x, x, N_Constraints)
    @timeit to "ctr mixed both" add_constraints_2(model, x, y, N_Constraints)
    @timeit to "ctr mixed2 vars" add_constraints_3(model, y, y, N_Constraints)
    @timeit to "ctr mixed2 params" add_constraints_3(model, x, x, N_Constraints)
    @timeit to "ctr mixed2 both" add_constraints_3(model, x, y, N_Constraints)

    data = ParameterJuMP._getparamdata(model)::ParameterJuMP.ParameterData
    @timeit to "sync" ParameterJuMP.sync(data)

    @timeit to "solve" optimize!(model)

    @timeit to "duals" get_duals(model, x)
end

to = TimerOutput()
reset_timer!(to)

bench_create_param(2, 2, 2)

reset_timer!(to)

bench_create_param(1000, 1000, 1000)

show(to)


text:
#addprocs(2)
#@show nprocs()

using RoME
using Test

@error("must restore testG2oParser.jl")
# "testG2oParser.jl";  ]

testfiles = [
  "testInflation380.jl";
  "testPoint2Point2.jl";
  "testParametric.jl";
  "testTreeInitCommonMsg_IIF913.jl";
  "threeDimLinearProductTest.jl";
  "testPose3Pose3NH.jl";
  "testBeehive2D_CliqByCliq.jl";      # special case debugging
  "testhigherdimroots.jl";
  "testManifoldsPose2Equivalent.jl";
  "testDidsonFunctions.jl";
  "testPoint2Point2Init.jl";
  "testBasicPose2Stationary.jl";
  "TestPoseAndPoint2Constraints.jl";
  "testPartialRangeCrossCorrelations.jl";
  "testDynPoint2D.jl";
  "testBearingRange2D.jl";
  "testDeltaOdo.jl";
  "testFixedLagFG.jl";
  "testMultimodalRangeBearing.jl";
  "testDynPose2D.jl";
  "testPartialXYH.jl";
  "testpartialpose3.jl";
  "testpackingconverters.jl";
  "TestDefaultFGInitialization.jl";
  "testAccumulateFactors.jl";
  "testDeadReckoningTether.jl";
  "testFluxModelsPose2.jl";
]

## Tests not ready yet
# "HexagonalLightGraphs.jl"
# "testCameraFunctions.jl"
# "testmultiplefeatures.jl"
# "testPoint2Point2WorldBearing.jl";  # deprecate

for testf in testfiles
  println("[TEST] $testf")
  include(testf)
  println("[SUCCESS]")
  println()
  println()
  println()
end


text:
include("low_level.jl")
include("build_stack.jl")
include("inputs.jl")
include("github_34.jl")
include("validate.jl")

max_forks_repo_name: HolyLab/ImagineInterface


text:
using OrdinaryDiffEq, DiffEqDevTools, Sundials, ModelingToolkit, ODEInterfaceDiffEq,
      Plots, DASSL, DASKR
using LinearAlgebra

@variables t y1(t)=1.0 y2(t)=2.0 y3(t)=3.0
@parameters p1=77.27 p2=8.375e-6 p3=0.161
D = Differential(t)
eqs = [
  D(y1) ~ p1*(y2+y1*(1-p2*y1-y2))
  D(y2) ~ (y3-(1+y1)*y2)/p1
  D(y3) ~ p3*(y1-y3)
]
@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)
mmprob = ODEProblem(sys,[],(0.0,30.0))
daeprob = DAEProblem(sys,[D(y1)=>77.26935286375,
                              D(y2)=>-0.012941633234114146,
                              D(y3)=>-0.322],[],(0.0,30.0))
odaeprob = ODAEProblem(simpsys,[],(0.0,30.0))

ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14);
ode_ref_sol = solve(odaeprob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14);

probs = [mmprob,daeprob,odaeprob]
refs = [ref_sol,ref_sol,ode_ref_sol];

plot(ref_sol)

abstols = 1.0 ./ 10.0 .^ (6:9)
reltols = 1.0 ./ 10.0 .^ (2:5);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()),
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
]

wp = WorkPrecisionSet(probs,abstols,reltols,setups;print_names=true,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
          Dict(:prob_choice => 3, :alg=>TRBDF2()),
          Dict(:prob_choice => 3, :alg=>KenCarp4()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:8)
reltols = 1.0 ./ 10.0 .^ (2:4);
setups = [Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 2, :alg=>DASSL.dassl()),
          Dict(:prob_choice => 2, :alg=>DASKR.daskr()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:9)
reltols = 1.0 ./ 10.0 .^ (2:5);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()), # too slow
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          ]
gr()
wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:9)
reltols = 1.0 ./ 10.0 .^ (2:5);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
          Dict(:prob_choice => 3, :alg=>TRBDF2()),
          Dict(:prob_choice => 3, :alg=>KenCarp4()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:12)
reltols = 1.0 ./ 10.0 .^ (4:9)

setups = [Dict(:prob_choice => 1, :alg=>Rodas5()),
          Dict(:prob_choice => 3, :alg=>Rodas5()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()),
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          ]
gr()
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
#
# Conditional Sampler
#
# chooses sampler based on values of a parent choice point
#
# implementation notes:
# 	* to facilitate EDAs using model-building (i.e. selection of parents), we implicitly constrain all the underlying samplers
# 	  to be of the same type (i.e. only the parameters differ)
# 	* we permit only one parent to be specified; to support models with multiple parents, the samplers can be nested:
#     ConditionalSampler(ConditionalSampler, ...)
#   * we permit parentcpid to be passed as "nothing" (which in effect makes the sampler unconditional as it simply calls the specified default sampler)
#     since often the parent cpid may not be known when the sampler is initially constructed.
#	* parentrecency indicates which invocation of the parent choice point to consider: 1 is the most recent, 2 is the second most recent etc.
#

type ConditionalSampler <: ModifyingSampler
	defaultsampler::Sampler
	parentcpid
	parentrecency::Int
	parentrecursiondepth # if nothing, then parent's recursion depth is ignored during sampling
	restricttoancestors::Bool
	conditionalsamplers::Dict{Any,Sampler}  # TODO: type of key is really Union{Real,Void} (since parent gns can be "nothing") rather than Any - is it worth specifying?
	function ConditionalSampler(defaultsampler::Sampler, parentcpid=nothing, parentrecency::Int=1, parentrecursiondepth=nothing, restricttoancestors=false, parentgns=Vector{Real}(), conditionalparams=Vector{Vector{Float64}}())
		# default parentcpid of nothing is for ease of construction as parent cpids may not be known at construction time
		# it turns the sampler into a non-conditional sampler
		(parentcpid != nothing) || isempty(parentgns) || error("Parent Gödel number cannot be specified if no parent cpid is given ")
		(parentrecency > 0) || error("Parent recency must be an integer 1 or more")
		(parentrecursiondepth == nothing) || ((typeof(parentrecursiondepth) <: Integer) && (parentrecursiondepth >= 1)) || error("Parent recursion must be nothing, or an integer 1 or more")
		(parentrecursiondepth == nothing) || warn("for parent recursion depth to take effect, ensure distinguishconditionalitybyrecursiondepth is set to true at model level")
		length(parentgns) >= length(conditionalparams) || error("More conditional parameter settings specified than parent Gödel numbers")
		conditionalsamplers = Dict{Any,Sampler}()
		for i in 1:length(parentgns)
			!haskey(conditionalsamplers, parentgns[i]) || error("Parent Gödel number is specified more than once")
			conditionalsampler = deepcopy(defaultsampler)
			# if params are specified, then set them, otherwise use default parameters
			if i <= length(conditionalparams)
				setparams!(conditionalsampler, conditionalparams[i])
			end
			conditionalsamplers[parentgns[i]] = conditionalsampler
		end
		new(defaultsampler, parentcpid, parentrecency, parentrecursiondepth, restricttoancestors, conditionalsamplers)
	end
end

# note: in the following methods, we could make assumption that numparams for each conditional sampler is the same as default sampler
# (which is what the constructor enforces), but where possible (generally, without too much of a performance impact) we avoid this assumption
# so that the code is a little more robust

# note: the code assumes that order of samplers returned by values(conditionalsamplers) is stable

# numparams(s::ConditionalSampler) = numparams(s.defaultampler) + sum(numparams, s.conditionalsamplers)

function paramranges(s::ConditionalSampler)
	pr = paramranges(s.defaultsampler)
	for conditionalsampler in values(s.conditionalsamplers)
		pr = [pr; paramranges(conditionalsampler)]
	end
	pr
end

function setparams!(s::ConditionalSampler, params)
	nparams = numparams(s)
	length(params) == nparams || error("expected $(nparams) parameters but got $(length(params))")
	paramstart = 1
	paramcount = numparams(s.defaultsampler)
	setparams!(s.defaultsampler, params[paramstart:(paramstart+paramcount-1)])
	for conditionalsampler in values(s.conditionalsamplers)
		paramstart += paramcount
		paramcount = numparams(conditionalsampler)
		setparams!(conditionalsampler, params[paramstart:(paramstart+paramcount-1)])
	end
end

function getparams(s::ConditionalSampler)
	ps = getparams(s.defaultsampler)
	for conditionalsampler in values(s.conditionalsamplers)
		ps = [ps; getparams(conditionalsampler)]
	end
	ps
end

function sample(s::ConditionalSampler, support, cc::ChoiceContext)
	# extract Gödel number of parent choice point at its most recent invocation
	parentgn = nothing
	if s.parentcpid != nothing
		recency = 0
		# if we are restricting conditionality to choice points in the same rule or its
		# ancestors in the execution tree, then we use ancestor cp sequence numbers to move
		# quickly through the ancestors in the trace
		# otherwise we simply step backwards through trace
		i = s.restricttoancestors ? getimmediateancestorcpseqnumber(cc.derivationstate) : length(cc.derivationstate.cmtrace)
		while (parentgn == nothing) && (i != nothing) && (i >= 1)
			(earliercpid, earliertrace) = cc.derivationstate.cmtrace[i]
			if (earliercpid == s.parentcpid) &&
				((s.parentrecursiondepth == nothing) || (earliertrace[:rcd] == s.parentrecursiondepth))
			# match parent based on cpid
			# AND if we are distinguishing conditionality based on the recursion depth of the parent choice point,
			# match recursion depth of that choice point in its trace
				recency += 1
				if recency == s.parentrecency
					parentgn = earliertrace[:gdl]
					@assert parentgn == cc.derivationstate.godelsequence[i]
				end
			end
			i = s.restricttoancestors ? earliertrace[:acs] : (i - 1)
		end
	end
	# if there is a sampler specified for the parent value, then use it, otherwise revert to default
	if (s.parentcpid != nothing) && (haskey(s.conditionalsamplers, parentgn))
		x, trace = sample(s.conditionalsamplers[parentgn], support, cc)
	else
		x, trace = sample(s.defaultsampler, support, cc)
	end
	x, Dict{Symbol, Any}(:par=>parentgn, :sub=>trace)
end

# TODO for consistency with other sub-sampler, we attempt re-estimation for all
# subsamplers, if we can tell then have too few (or no) traces here: we let the sub-sampler decide what
# to do in this situation.  This may be a (small?) performance overhead
#
# NOTE: if model is changed (i.e. a different parent is identified), then traces will be invalid
# and this method must NOT be used
#
function estimateparams!(s::ConditionalSampler, traces)
	# first divide out the traces according to which subsampler was used
	defaultsamplertraces = Vector{Dict}()
	conditionalsamplertraces = Dict{Any,Vector{Dict}}()
	for parentgn in keys(s.conditionalsamplers)
		conditionalsamplertraces[parentgn] = Vector{Dict}()
	end
	for trace in traces
		parentgn = trace[:par]
		if haskey(s.conditionalsamplers, parentgn)
			push!(conditionalsamplertraces[parentgn], trace[:sub])
		else
			push!(defaultsamplertraces, trace[:sub])
		end
	end
	# then re-estimate the subsamplers
	estimateparams!(s.defaultsampler, defaultsamplertraces)
	for parentgn in keys(s.conditionalsamplers)
		estimateparams!(s.conditionalsamplers[parentgn], conditionalsamplertraces[parentgn])
	end
end

function amendtrace(s::ConditionalSampler, trace, x)
	parentgn = trace[:par]
	if haskey(s.conditionalsamplers, parentgn)
		amendtrace(s.conditionalsamplers[parentgn], trace[:sub], x)
	else
		amendtrace(s.defaultsampler, trace[:sub], x)
	end
end

#
#
# K2 Metric: consider node A which has r distinct values
# par(A) are parents of A, which have q is distinct instantions of parents
# N_ij is number of cases where A takes i^th value and parents take j^th value (N_.j is sum of this over i)
# then:
#   log K2[A|par(A)] = sum_j=1^q log[(r-1)! / (N_.j + r -1)!] + sum_j=1^q sum_i=1^r log[N_ij!]
#
function determineparentusingK2(s::ConditionalSampler, gnhistories)

	numhistories = length(gnhistories)
	# the history for the current cp (i.e. "node A") which is the last value in each history
	ahistory = map(x->x[end], gnhistories)
	# A_values are values that current cp takes
	avalues = unique(ahistory)
	# println("A takes values: $(avalues)")
	r = length(avalues)
	# println("Therefore r is: $(r)")
	lfactrminus1 = lfact(r-1)

	# calculate value for no parent (effectively a parent taking all the same value)
	bestparindices = Any[nothing]
	bestlogK2 = lfactrminus1 - lfact(numhistories + r - 1)
	for i in 1:r
		Ni = count(x->x==avalues[i], ahistory)
		bestlogK2 += lfact(Ni)
	end
	# println("No parent gives a K2 value of: $(bestlogK2)")

	### TEMP
	# warn("removing calculation of sensible minimum K2")
	# bestlogK2 = -Inf
	###

	if numhistories > 0

		# consider each possible parent in turn
		for parindex in 1:(length(gnhistories[1]) - 1) # -1 since last entry is value of current cp
			# println("Considering potential parent #$(parindex) ...")

			# the history of the i^th possible parent is the i^th entry
			parhistory = map(x->x[parindex], gnhistories)
			parvalues = unique(parhistory)
			# println("par(A) takes values: $(parvalues)")
			q = length(parvalues)
			# println("Therefore q is: $(q)")

			logK2 = 0.0

			for j in 1:q
				Ndotj = count(x->x==parvalues[j], parhistory)
				# println("For j=$(j), N_.j is: $(Ndotj)")
				logK2 += lfactrminus1 - lfact(Ndotj + r - 1)
				# lfact is log factorial approximation
				# TODO: since current always integer, could do a lookup? will this be faster than Julia's lfact?

				for i in 1:r
					Nij = count(k->((ahistory[k]==avalues[i]) && (parhistory[k]==parvalues[j])), 1:numhistories)
					# println("For j=$(j), i=$(i), N_ij is: $(Nij)")
					logK2 += lfact(Nij)
				end

			end

			# println("Potential parent $(parindex) has logK2 of: $(logK2)")
			if logK2 > bestlogK2
				# strictly greater than so simplicitly of network is encouraged
				bestparindices = [parindex]
				bestlogK2 = logK2
			elseif logK2 == bestlogK2
				push!(bestparindices, parindex)
			end

		end

	end

	(nothing in bestparindices) ? nothing : bestparindices[rand(1:length(bestparindices))]

end

function estimateconditionalmodel!(s::ConditionalSampler, cplabels, gnhistories, traces)

	# assign new parent
	parindex = determineparentusingK2(s, gnhistories)
	if parindex != nothing
		cpid, recursiondepth, recency, restricttoancestors = cplabels[parindex]
	else
		cpid, recursiondepth, recency, restricttoancestors = nothing, nothing, 1, false
	end
	s.parentcpid = cpid
	s.parentrecency = recency
	s.parentrecursiondepth = recursiondepth
	s.restricttoancestors = restricttoancestors

	# println("******* New parent cpid is index $(parindex) : $(s.parentcpid)")

	# defaultsampler
	# TODO, if there is a parent, since we will account for all values of the parents, there is nothing left
	# to estimate the default: what to do? Currently we do this for all ...

	# we use the following to filter histories passed to any sub-sampler to remove the chosen cpid
	cplabelmask = convert(Vector{Bool}, map(c->c!=parindex, 1:length(cplabels)))
	historymask = [cplabelmask; true]
	# we add a true to the cpmask to retain final value which is the current cp Gödel number

	defaulttraces = map(t->t[:sub], traces)
	if supportsconditionalmodelestimation(s.defaultsampler)
		# pass all histories and traces to default sampler
		defaultgnhistories = map(h->h[historymask], gnhistories)
		estimateconditionalmodel!(s.defaultsampler, cplabels[cplabelmask], defaultgnhistories, defaulttraces)
	else
		estimateparams!(s.defaultsampler, defaulttraces)
	end

	# remove existing conditional samplers
	s.conditionalsamplers = Dict{Any,Sampler}()

	if parindex != nothing
		# if parent identified, create and estimate a sampler for each parent values

		parhistory = map(x->x[parindex], gnhistories)
		parvalues = unique(parhistory)

		for parvalue in parvalues

			conditionalsampler = deepcopy(s.defaultsampler)
			s.conditionalsamplers[parvalue] = conditionalsampler

			parentvaluemask = convert(Vector{Bool}, map(v->v==parvalue, parhistory))
			# we use the history to determine which histories are relevant
			# note: we can't use the parent value stored with the trace since the parent may have changed

			conditionaltraces = map(x->x[:sub], traces[parentvaluemask])
			if supportsconditionalmodelestimation(conditionalsampler)
				# pass histories and traces filtered by mask where parent value matches the condition on the sampler
				# we also filter the histories to remove parent index just chosen
				conditionalgnhistories = map(h->h[historymask], gnhistories[parentvaluemask])
				# println("<<<<<<< START nested")
				estimateconditionalmodel!(conditionalsampler, cplabels[cplabelmask], conditionalgnhistories, conditionaltraces)
				# println("<<<<<<< END nested")
			else
				estimateparams!(conditionalsampler, conditionaltraces)
			end

		end

	end

end

# pretty print the model
function show(io::IO, s::ConditionalSampler, indentdepth::Int=1)
	indent = " "^(10indentdepth)
	println(io, getsamplertypename(s) * " Parent: " * ((s.parentcpid == nothing) ? "(none)" : ("$(hex(s.parentcpid)) Recency: $(s.parentrecency) Recursion Depth: "
		* ((s.parentrecursiondepth == nothing) ?  "(any)" : "$(s.parentrecursiondepth)") * " Restrict to Ancestors: " * (s.restricttoancestors ?  "yes" : "no"))))
	for parentvalue in sort(collect(keys(s.conditionalsamplers)),lt=(x,y)->(x==nothing) ? false : (y==nothing) ? true : isless(x,y))
		print(io, indent * (parentvalue == nothing ? "(none)" : "$(parentvalue)") * ": ")
		show(io, s.conditionalsamplers[parentvalue], indentdepth+1)
	end
	print(io, indent * "default: ")
	show(io, s.defaultsampler, indentdepth+1)
end


text:
using CompScienceMeshes
using BEAST

Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
X = raviartthomas(Γ)

κ, η = 1.0, 1.0
t = Maxwell3D.singlelayer(wavenumber=κ)
E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
# E = -η/(im*κ)*BEAST.CurlCurlGreen(κ, ẑ, point(2,0,0))
e = (n × E) × n

@hilbertspace j
@hilbertspace k
efie = @discretise t[k,j]==e[k]  j∈X k∈X
u = gmres(efie)

include("utils/postproc.jl")
include("utils/plotresults.jl")


text:
module Baselet

include("specialized.jl")
include("utils.jl")
include("implementations.jl")

end # module


text:
using sem, Test, SafeTestsets

@safetestset "Koch" begin include("example_koch.jl") end
@safetestset "DefVars" begin include("example_defvars.jl") end


text:
function solve!(solver::AugmentedLagrangianSolver{T,S}) where {T,S}
	initialize!(solver)
    c_max::T = Inf

	conSet = get_constraints(solver)
	solver_uncon = solver.solver_uncon::S

	# Calculate cost
    J_ = get_J(get_objective(solver))
    J = sum(J_)

    for i = 1:solver.opts.iterations
		set_tolerances!(solver, solver_uncon, i)

        step!(solver)
        J = sum(J_)
        c_max = maximum(conSet.c_max)

        record_iteration!(solver, J, c_max)

        converged = evaluate_convergence(solver)
        if converged
            break
        end

        reset!(solver_uncon, false)
    end
    return solver
end

function initialize!(solver::AugmentedLagrangianSolver)
	set_verbosity!(solver.opts)
	clear_cache!(solver.opts)

	# Reset solver
    reset!(get_constraints(solver), solver.opts)
    solver.stats.iterations = 0
	solver.stats.iterations_total = 0

	# Calculate cost
    cost!(get_objective(solver), get_trajectory(solver))
end

function step!(solver::AugmentedLagrangianSolver)

    # Solve the unconstrained problem
    solve!(solver.solver_uncon)

    # Outer loop update
    dual_update!(solver)
    penalty_update!(solver)
    max_violation!(get_constraints(solver))

	# Reset verbosity level after it's modified
	set_verbosity!(solver.opts)
end

function record_iteration!(solver::AugmentedLagrangianSolver{T,S}, J::T, c_max::T) where {T,S}

    solver.stats.iterations += 1
    i = solver.stats.iterations
    solver.stats.iterations_total += solver.solver_uncon.stats.iterations
    solver.stats.c_max[i] = c_max
	solver.stats.cost[i] = J

	conSet = get_constraints(solver)
	max_penalty!(conSet)
	solver.stats.penalty_max[i] = maximum(conSet.c_max)

	@logmsg OuterLoop :iter value=i
	@logmsg OuterLoop :total value=solver.stats.iterations_total
	@logmsg OuterLoop :cost value=J
    @logmsg OuterLoop :c_max value=c_max
	if solver.opts.verbose
		print_level(OuterLoop)
	end
end

function set_tolerances!(solver::AugmentedLagrangianSolver{T},
        solver_uncon::AbstractSolver{T},i::Int) where T
    if i != solver.opts.iterations
        solver_uncon.opts.cost_tolerance = solver.opts.cost_tolerance_intermediate
        solver_uncon.opts.gradient_norm_tolerance = solver.opts.gradient_norm_tolerance_intermediate
    else
        solver_uncon.opts.cost_tolerance = solver.opts.cost_tolerance
        solver_uncon.opts.gradient_norm_tolerance = solver.opts.gradient_norm_tolerance
    end

    return nothing
end

function evaluate_convergence(solver::AugmentedLagrangianSolver)
	i = solver.stats.iterations
    solver.stats.c_max[i] < solver.opts.constraint_tolerance ||
		solver.stats.penalty_max[i] >= solver.opts.penalty_max
end

"General Dual Update"
function dual_update!(solver::AugmentedLagrangianSolver) where {T,Q,N,M,NM}
    conSet = get_constraints(solver)
    for i in eachindex(conSet.constraints)
        dual_update!(conSet.constraints[i])
    end
end

"Dual Update for Equality Constraints"
function dual_update!(con::ConstraintVals{T,W,C}) where
		{T,W,C<:AbstractConstraint{Equality}}
	λ = con.λ
	c = con.vals
	μ = con.μ
	λ_max = con.params.λ_max
	for i in eachindex(con.inds)
		λ[i] = clamp.(λ[i] + μ[i] .* c[i], -λ_max, λ_max)
	end
end

"Dual Update for Inequality Constraints"
function dual_update!(con::ConstraintVals{T,W,C}) where
		{T,W,C<:AbstractConstraint{Inequality}}
	λ = con.λ
	c = con.vals
	μ = con.μ
	for i in eachindex(con.inds)
		λ[i] = clamp.(λ[i] + μ[i] .* c[i], 0.0, con.params.λ_max)
	end
end

"General Penalty Update"
function penalty_update!(solver::AugmentedLagrangianSolver)
    conSet = get_constraints(solver)
    for i in eachindex(conSet.constraints)
        penalty_update!(conSet.constraints[i])
    end
end

"Penalty Update for ConstraintVals"
function penalty_update!(con::ConstraintVals{T}) where T
	ϕ = con.params.ϕ
	μ = con.μ
	for i in eachindex(con.inds)
		μ[i] = clamp.(ϕ * μ[i], 0.0, con.params.μ_max)
	end
end


text:
# Benchmarks.jl
#
# Run some benchmarks for FluxArchitectures

# Settings
const RUN_GPU = true
const RUN_CPU = true

# Make sure all the required packages are available
cd(@__DIR__)
using Pkg; Pkg.activate(".")
Pkg.instantiate()

@info "Loading packages"
include("../src/FluxArchitectures.jl")
using .FluxArchitectures
using BenchmarkTools
using DataFrames
using Random

# Initialize results dataframe
results = DataFrame(Name=String[], CPU=Union{Missing,Float32}[], GPU=Union{Missing,Float32}[])

# Set up data and loss
poollength = 10
horizon = 6
datalength = 500
input_cpu, target_cpu = get_data(:solar, poollength, datalength, horizon)
input_gpu, target_gpu = gpu(input_cpu), gpu(target_cpu)
inputsize = size(input_cpu, 1)
function loss(x, y)
    Flux.reset!(model)
    return Flux.mse(model(x), permutedims(y))
end

## DARNN ----------------------------------
@info "DARNN"
push!(results,["DARNN", missing, missing])
encodersize = 10
decodersize = 10
if RUN_GPU
    model = DARNN(inputsize, encodersize, decodersize, poollength, 1) |> gpu
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_gpu, target_gpu), 5),
            ADAM(0.007))
    results[end,:GPU] = time
end
if RUN_CPU
    model = DARNN(inputsize, encodersize, decodersize, poollength, 1)
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_cpu, target_cpu), 5),
            ADAM(0.007))
    results[end,:CPU] = time
end

## DSANet ----------------------------------
@info "DSANet"
push!(results,["DSANet", missing, missing])
local_length = 3
n_kernels = 3
d_model = 4
hiddensize = 1
n_layers = 3
n_head = 2
if RUN_GPU
    Random.seed!(123)
    model = DSANet(inputsize, poollength, local_length, n_kernels, d_model,
               hiddensize, n_layers, n_head) |> gpu
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_gpu, target_gpu), 5),
            ADAM(0.007))
    results[end,:GPU] = time
end
if RUN_CPU
    Random.seed!(123)
    model = DSANet(inputsize, poollength, local_length, n_kernels, d_model,
               hiddensize, n_layers, n_head)
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_cpu, target_cpu), 5),
            ADAM(0.007))
    results[end,:CPU] = time
end

## LSTNet ----------------------------------
@info "LSTNet"
push!(results,["LSTNet", missing, missing])
convlayersize = 2
recurlayersize = 3
skiplength = 120
if RUN_GPU
    model = LSTnet(inputsize, convlayersize, recurlayersize, poollength, skiplength,
        init=Flux.zeros, initW=Flux.zeros) |> gpu
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_gpu, target_gpu), 5),
            ADAM(0.007))
    results[end,:GPU] = time
end
if RUN_CPU
    model = LSTnet(inputsize, convlayersize, recurlayersize, poollength, skiplength,
        init=Flux.zeros, initW=Flux.zeros)
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_cpu, target_cpu), 5),
            ADAM(0.007))
    results[end,:CPU] = time
end

## TPA-LSTM ----------------------------------
@info "TPA-LSTM"
push!(results,["TPA-LSTM", missing, missing])
hiddensize = 10
layers = 2
filternum = 32
filtersize = 1
if RUN_GPU
    model = TPALSTM(inputsize, hiddensize, poollength, layers, filternum, filtersize) |> gpu
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_gpu, target_gpu), 5),
            ADAM(0.007))
    results[end,:GPU] = time
end
if RUN_CPU
    model = TPALSTM(inputsize, hiddensize, poollength, layers, filternum, filtersize)
    time = @belapsed Flux.train!(loss, Flux.params(model),Iterators.repeated((input_cpu, target_cpu), 5),
            ADAM(0.007))
    results[end,:CPU] = time
end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
module MinresQlp
export minresqlp_iterable, minresqlp, minresqlp!
using Printf
import IterativeSolvers: zerox, ConvergenceHistory, reserve!, nextiter!, setconv
import IterativeSolvers: shrink!
import LinearAlgebra: BLAS.axpy!, givensAlgorithm, norm, mul!, axpy!
import Base: iterate

mutable struct MINRESQLPIterable{matT, solT, vecT <: DenseVector, shiftT<: Number, realT <: Real}
    A::matT
    shift::shiftT
    skew_hermitian::Bool
    x::solT

    # Krylov basis vectors
    r1::vecT
    r2::vecT
    r3::vecT

    #
    w::vecT
    wl::vecT
    xl2::vecT
    wl2::vecT
    v::vecT

    #flags
    flag0::Int
    flag::Int
    iters::Int
    QLPiter::Int
    beta::realT
    tau::realT
    taul::realT

    phi::realT
    beta1::realT #almost useless
    betan::realT
    gmin::realT
    cs::realT
    sn::realT
    cr1::realT
    sr1::realT
    cr2::realT
    sr2::realT
    dltan::realT
    eplnn::realT
    gama::realT
    gamal::realT
    gamal2::realT
    gama_QLP::realT
    gamal_QLP::realT
    eta::realT
    etal::realT
    etal2::realT
    vepln::realT
    veplnl::realT
    veplnl2::realT
    vepln_QLP::realT
    ul3::realT
    ul2::realT
    ul::realT
    u::realT
    ul_QLP::realT
    u_QLP::realT

    rnorm::realT
    Arnorm::realT
    xnorm::realT
    relAres::realT

    xl2norm::realT
    Axnorm::realT
    Anorm::realT
    Acond::realT
    relres::realT
    gminl::realT

    noprecon::Bool

    #lims
    TranCond::realT
    maxxnorm::realT
    Acondlim::realT

    # Bookkeeping
    mv_products::Int
    maxiter::Int
    tolerance::realT
    resnorm::realT
end

function minresqlp_iterable!(x, A, b;
    initially_zero::Bool = false,
    skew_hermitian::Bool = false,
    tol = sqrt(eps(real(eltype(b)))),
    maxiter = size(A, 2),
    shift   = 0.0,
    M       = nothing,
    TranCond=10e6,
    maxxnorm=10e6,
    Acondlim=10e14)

    T = eltype(x)
    rT= real(T)
    HessenbergT = skew_hermitian ? T : real(T)

    n = length(b)
    r1 = similar(b);    r1.=real(T)(0)
    r2 = similar(b);    copyto!(r2, b)
    r3 = similar(b);    copyto!(r3, b)
    beta1 = norm(r2)

    if isa(M, Nothing)
        noprecon = true
    else
        noprecon = false
        precond!(r3, M, r2)
        beta1 = real(r3'*r2)  #beta1 = r3.T.dot(r2) #teta
        if beta1 < 0.0
            print("Error: M is indefinite!")
        else
            beta1 = sqrt(beta1)
        end
    end
    relres = beta1/(beta1+1e-50)

    #x   = zeros(T, n)
    w   = zeros(T, n)
    wl  = zeros(T, n)
    xl2 = zeros(T, n)
    wl2 = zeros(T, n)
    v   = zeros(T, n)

    MINRESQLPIterable{typeof(A), typeof(x), typeof(b), typeof(shift), rT}(
        A, shift, skew_hermitian, x,
        r1, r2, r3,
        w, wl, xl2, wl2, v,
        -2, -2, 0, 0, zero(rT), zero(rT), zero(rT),   # flag0, flag, iters, qlpiter, beta, tau, taul
        beta1, beta1, beta1, zero(rT), -one(rT), zero(rT), -one(rT), zero(rT), #phi, beta1, betan, gmin, cs, sn, cr1, sr1
        -one(rT), zero(rT), zero(rT), zero(rT),  #cr2, sr2, dltan, eplnn,
        zero(rT), zero(rT), zero(rT), zero(rT), zero(rT), #gama, gamal, gamal2, gama_QLP, gamal_QLP
        zero(rT), zero(rT), zero(rT), zero(rT), zero(rT), zero(rT), zero(rT), #eta etal etal2 vepln veplnl veplnl2, vepln_QLP
        zero(rT), zero(rT), zero(rT), zero(rT), zero(rT), zero(rT), # ul3, ul2 ul u ul_QLP u_QLP
        beta1, zero(rT), zero(rT), zero(rT), # rnorm xnorm Arnorm relAres
        zero(rT), zero(rT), zero(rT), one(rT), #xl2norm axnorm anorm acond
        relres, zero(rT), #relres gminl,
        noprecon,
        TranCond, maxxnorm, TranCond,
        0, maxiter, tol, rT(123)
    )
end

converged(m::MINRESQLPIterable) = m.resnorm ≤ m.tolerance || m.flag != m.flag0

start(::MINRESQLPIterable) = 1

done(m::MINRESQLPIterable, iteration::Int) = iteration > m.maxiter || converged(m)

function iterate(m::MINRESQLPIterable, iteration::Int=start(m))
    if done(m, iteration) return nothing end

    m.iters += 1
    betal = m.beta
    m.beta = m.betan
    m.v .= m.r3./m.beta        # v = r3/beta
    #r3 = A*v            # r3 = Ax(A, v)
    mul!(m.r3, m.A, m.v)

    m.shift != 0.0 && axpy!(-m.shift, m.v, m.r3) # r3 .= r3 - shift*v
    #TODO Fix right mult
    m.iters > 1 && axpy!(-m.beta/betal, m.r1, m.r3) #r3 = r3 - r1*beta/betal

    alfa = real(m.r3'*m.v)                          # alfa = np.real(r3.T.dot(v))
    axpy!(-alfa/m.beta, m.r2, m.r3)   # r3 = r3 - r2*alfa/beta #TODO Fix right mult
    copyto!(m.r1, m.r2)                             # r1 = r2
    copyto!(m.r2, m.r3)                             # r2 = r3

    if m.noprecon
        m.betan = norm(m.r3)
        if m.iters == 1
            if m.betan == 0.0
                if m.alfa == 0.0
                    m.flag = 0
                    return
                else
                    m.flag = -1
                    mul!(m.x, m.b, inv(alfa)) #x = b/alfa
                    return
                end
            end
        end
    else#=
        precond!(r3, M, r2)
        m.betan = real(r2'*r3)
        if m.betan > 0.0
            m.betan = sqrt(m.betan)
        else
            print("Error: M is indefinite or singular!")
        end=#
        @error "Not Implemented"
    end
    pnorm = sqrt(betal^2 + alfa^2 + m.betan^2)

    #previous left rotation Q_{k-1}
    dbar = m.dltan
    dlta = m.cs*dbar + m.sn*alfa
    epln = m.eplnn
    gbar = m.sn*dbar - m.cs*alfa
    m.eplnn = m.sn*m.betan
    m.dltan = -m.cs*m.betan
    dlta_QLP = dlta
    #current left plane rotation Q_k
    gamal3 = m.gamal2
    m.gamal2 = m.gamal
    m.gamal = m.gama
    m.cs, m.sn, m.gama = SymGivens(gbar, m.betan)
    gama_tmp = m.gama
    taul2 = m.taul
    m.taul = m.tau
    m.tau = m.cs*m.phi
    m.Axnorm = sqrt(m.Axnorm^2 + m.tau^2)
    m.phi = m.sn*m.phi

    #previous right plane rotation P_{k-2,k}
    if m.iters > 2
        m.veplnl2 = m.veplnl
        m.etal2 = m.etal
        m.etal = m.eta
        dlta_tmp = m.sr2*m.vepln - m.cr2*dlta
        m.veplnl = m.cr2*m.vepln + m.sr2*dlta
        dlta = dlta_tmp
        m.eta = m.sr2*m.gama
        m.gama = -m.cr2 *m.gama
    end
    #current right plane rotation P{k-1,k}
    if m.iters > 1
        m.cr1, m.sr1, m.gamal = SymGivens(m.gamal, dlta)
        m.vepln = m.sr1*m.gama
        m.gama = -m.cr1*m.gama
    end

    #327
    #update xnorm
    xnorml = m.xnorm
    ul4 = m.ul3
    m.ul3 = m.ul2
    if m.iters > 2
        m.ul2 = (taul2 - m.etal2*ul4 - m.veplnl2*m.ul3)/m.gamal2
    end
    if m.iters > 1
        m.ul = (m.taul - m.etal*m.ul3 - m.veplnl * m.ul2)/m.gamal
    end
    xnorm_tmp = sqrt(m.xl2norm^2 + m.ul2^2 + m.ul^2)
    if abs(m.gama) > floatmin(real(eltype(m.A))) && xnorm_tmp < m.maxxnorm
        m.u = (m.tau - m.eta*m.ul2 - m.vepln*m.ul)/m.gama
        if sqrt(xnorm_tmp^2 + m.u^2) > m.maxxnorm
            m.u = 0
            m.flag = 6
        end
    else
        m.u = 0
        m.flag = 9
    end
    m.xl2norm = sqrt(m.xl2norm^2 + m.ul2^2)
    m.xnorm = sqrt(m.xl2norm^2 + m.ul^2 + m.u^2)

    #update w&x
    #Minres
    if (m.Acond < m.TranCond) && m.flag != m.flag0 && m.QLPiter == 0
        copyto!(m.wl2, m.wl)  # wl2 = wl
        copyto!(m.wl, m.w)    # wl = w
        m.w .= (m.v .- epln .* m.wl2 .- dlta_QLP .* m.wl)./gama_tmp
        if m.xnorm < m.maxxnorm
            m.x .+= m.tau .* m.w
        else
            m.flag = 6
        end
    #Minres-QLP
    else
        m.QLPiter += 1
        if m.QLPiter == 1
            if (m.iters > 1)    # construct w_{k-3}, w_{k-2}, w_{k-1}
                if m.iters > 3
                    m.wl2 .= gamal3.*m.wl2 .+ veplnl2.*m.wl .+ etal.*w
                end
                if m.iters > 2
                    m.wl .= gamal_QLP.*m.wl .+ vepln_QLP.*w
                end
                lmul!(gama_QLP, w)          #w = gama_QLP*w
                xl2 .= x .- m.wl.*ul_QLP .- w.*u_QLP
            end
        end

        # 369
        if m.iters == 1
            m.wl2 .= m.wl
            m.wl .=    m.v .* m.sr1
            m.w    .= .- m.v .* m.cr1
        elseif m.iters == 2
            m.wl2 .= m.wl
            m.wl  .= m.w .* m.cr1 .+ m.v .* m.sr1
            m.w   .= m.w .* m.sr1 .- m.v .* m.cr1
        else
            m.wl2 .= m.wl
            m.wl  .= m.w
            m.w   .= m.wl2 .* m.sr2 .- m.v .* m.cr2
            m.wl2 .= m.wl2 .* m.cr2 .+ m.v .* m.sr2
            m.v   .= m.wl  .* m.cr1 .+ m.w .* m.sr1
            m.w   .= m.wl  .* m.sr1 .- m.w .* m.cr1
            m.wl  .= m.v
        end
        m.xl2 .= m.xl2 .+ m.wl2 .* m.ul2
        m.x   .= m.xl2 .+ m.wl  .* m.ul .+ m.w .* m.u
    end

    # 388
    #next right plane rotation P{k-1,k+1}
    gamal_tmp = m.gamal
    m.cr2, m.sr2, m.gamal = SymGivens(m.gamal, m.eplnn)
    #transfering from Minres to Minres-QLP
    m.gamal_QLP = gamal_tmp
    m.vepln_QLP = m.vepln
    m.gama_QLP = m.gama
    m.ul_QLP = m.ul
    m.u_QLP = m.u
    ## Estimate various norms
    abs_gama = abs(m.gama)
    Anorml = m.Anorm
    m.Anorm = max(m.Anorm, pnorm, m.gamal, abs_gama)  #Anorm = maximum([Anorm, pnorm, gamal, abs_gama])
    if m.iters == 1
        m.gmin = m.gama
        m.gminl = m.gmin
    elseif m.iters > 1
        gminl2 = m.gminl
        m.gminl = m.gmin
        gmin = min(gminl2, m.gamal, abs_gama)   #gmin = minimum([gminl2, gamal, abs_gama])
    end

    #409
    Acondl = m.Acond
    m.Acond = m.Anorm / m.gmin
    rnorml = m.rnorm
    relresl = m.relres
    if m.flag != 9
        m.rnorm = m.phi
    end
    m.relres = m.rnorm / (m.Anorm * m.xnorm + m.beta1)
    rootl = sqrt(gbar ^2 + m.dltan ^ 2)
    m.Arnorm = rnorml * rootl
    m.relAres = rootl / m.Anorm
    ## See if any of the stopping criteria are satisfied.
    epsx = m.Anorm * m.xnorm * eps(typeof(m.xnorm))

    if (m.flag == m.flag0) || (m.flag == 9)
        t1 = 1 + m.relres
        t2 = 1 + m.relAres
        if m.iters >= m.maxiter
            m.flag = 8 #exit before maxit
        end
        if m.Acond >= m.Acondlim
            m.flag = 7 #Huge Acond
        end
        if m.xnorm >= m.maxxnorm
            m.flag = 6 #xnorm exceeded
        end
        if epsx >= m.beta1
            m.flag = 5 #x = eigenvector
        end
        if t2 <= 1
            m.flag = 4 #Accurate Least Square Solution
        end
        if t1 <= 1
            m.flag = 3 #Accurate Ax = b Solution
        end
        if m.relAres <= m.tolerance
            m.flag = 2 #Trustful Least Square Solution
        end
        if m.relres <= m.tolerance
            m.flag = 1 #Trustful Ax = b Solution
        end
    end
    if m.flag == 2 || m.flag == 4 || m.flag == 6 || m.flag == 7
        #possibly singular
        m.iters = m.iters - 1
        m.Acond = Acondl
        m.rnorm = rnorml
        m.relres = relresl
        #else show...
    end

    # The approximate residual is cheaply available
    m.resnorm = min(m.relres, m.relAres)

    min(m.relres, m.relAres), m.iters
end

function complete!(m::MINRESQLPIterable)
    Miter = m.iters - m.QLPiter

    mul!(m.r1, m.A, m.x)
    axpby!(m.shift, m.x, -1.0, m.r1)
    m.r1 .+= m.b

    m.rnorm = norm(m.r1)
    Arnorm = norm(m.A*m.r1 - m.shift*m.r1)
    m.xnorm = norm(m.x)
    relres = m.rnorm/(m.Anorm*m.xnorm + m.beta1)
    relAres = 0
    if rnorm > floatmin(real(eltype(A)))
        relAres = Arnorm/(Anorm*rnorm)
    end
end
"""
    minresqlp!(x, A, b; kwargs...) -> x, [history]

"""
function minresqlp!(x, A, b;
    skew_hermitian::Bool = false,
    verbose::Bool = false,
    log::Bool = false,
    tol = sqrt(eps(real(eltype(b)))),
    maxiter::Int = size(A, 2),
    initially_zero::Bool = false
)
    history = ConvergenceHistory(partial = !log)
    history[:tol] = tol
    log && reserve!(history, :resnorm, maxiter)
    log && reserve!(history, :rnorm, maxiter)
    log && reserve!(history, :Arnorm, maxiter)
    log && reserve!(history, :relres, maxiter)
    log && reserve!(history, :relAres, maxiter)
    log && reserve!(history, :Anorm, maxiter)
    log && reserve!(history, :Acond, maxiter)
    log && reserve!(history, :xnorm, maxiter)

    iterable = minresqlp_iterable!(x, A, b;
        skew_hermitian = skew_hermitian,
        tol = tol,
        maxiter = maxiter,
        initially_zero = initially_zero
    )

    if log
        history.mvps = iterable.mv_products
    end

    verbose && println("     iter     rnorm    Arnorm   relres   " *
                        "relAres    Anorm    Acond    xnorm")
    for (iteration, resnorm) = enumerate(iterable)
        m=iterable
        if log
            nextiter!(history, mvps = 1)
            push!(history, :resnorm, resnorm)
            push!(history, :rnorm, m.rnorm)
            push!(history, :Arnorm, m.Arnorm)
            push!(history, :relres, m.relres)
            push!(history, :relAres, m.relAres)
            push!(history, :Anorm, m.Anorm)
            push!(history, :Acond, m.Acond)
            push!(history, :xnorm, m.xnorm)
        end
        #verbose && @printf("%3d\t%1.2e\n", iteration, resnorm)
        verbose && @printf("%8g    %8.2e %8.2eD %8.2e %8.2eD %8.2e %8.2e  %8.2e\n",
                           m.iters, m.rnorm, m.Arnorm, m.relres, m.relAres, m.Anorm, #should be m.Anorml
                           m.Acond, m.xnorm)
    end

    verbose && println()
    log && setconv(history, converged(iterable))
    log && shrink!(history)

    log ? (iterable.x, history) : iterable.x
end

"""
    minres(A, b; kwargs...) -> x, [history]

Same as [`minres!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
minresqlp(A, b; kwargs...) = minresqlp!(zerox(A, b), A, b; initially_zero = true, kwargs...)

### ADDD
function SymGivens(a, b)
    if b == 0
        if a == 0
            c = 1
        else
            c = sign(a)
        end
        s = 0
        r = abs(a)
    elseif a == 0
        c = 0
        s = sign(b)
        r = abs(b)
    elseif abs(b) > abs(a)
        t = a / b
        s = sign(b) / sqrt(1 + t ^ 2)
        c = s * t
        r = b / s
    else
        t = b / a
        c = sign(a) / sqrt(1 + t ^ 2)
        s = c * t
        r = a / c
    end

    return c, s, r
end

end


text:
module SDDP

using Compat, Compat.LinearAlgebra

using TimerOutputs

# SDDP algorithm for StochOptInterface
using StochOptInterface
const SOI = StochOptInterface

using StructDualDynProg

# Utils
include("comp.jl")

# Path Sampler
include("sampler.jl")

include("path.jl")
include("algorithm.jl")

end


text:
# ============================================================
# Polynomial Methods
# ============================================================

include("poly_filter.jl")
include("poly_legendre.jl")
include("poly_jacobi.jl")
include("poly_lagrange.jl")
include("poly_triangle.jl")


text:
import Base.show

"""
**Show an occurrence**

    show(io::IO, o::GBIFRecord)

Displays the key, the taxon name, and the country of observation.
"""
function show(io::IO, o::GBIFRecord)
  println(io, "GBIF record $(o.key)\t$(o.taxon.name)\t($(o.country))")
end

"""
**Show several occurrences**

    show(io::IO, o::GBIFRecords)

Displays the total number, and the number of currently unmasked records.
"""
function show(io::IO, o::GBIFRecords)
  println(io, "GBIF records: viewing $(length(o)) [out of $(o.count)]")
end

"""
**Show a taxonomic record**

    show(io::IO, t::GBIFTaxon)

Displays the taxon name.
"""
function show(io::IO, t::GBIFTaxon)
  println(io, "GBIF taxon -- $(t.name)")
end


text:
"""
    AbstractMeshArray{T, N}

Subtype of AbstractArray{T, N}
"""
abstract type AbstractMeshArray{T, N} <: AbstractArray{T, N} end

"""
    gcmgrid

gcmgrid data structure. Available constructors:

```
gcmgrid(path::String, class::String, nFaces::Int,
        fSize::Array{NTuple{2, Int},1}, ioSize::Array{Int64,2},
        ioPrec::Type, read::Function, write::Function)
```
"""
struct gcmgrid
  path::String
  class::String
  nFaces::Int
  fSize::Array{NTuple{2, Int},1}
#  ioSize::NTuple{2, Int}
  ioSize::Array{Int64,2}
  ioPrec::Type
  read::Function
  write::Function
end

## concrete types and MeshArray alias:

include("gcmfaces_type.jl");
include("gcmarray_type.jl");
include("gcmvector_type.jl");

#MeshArray=gcmfaces
MeshArray=gcmarray

## Methods that apply to all AbstractMeshArray types

import Base: maximum, minimum, sum, fill, fill!

function maximum(a::AbstractMeshArray)
  c=-Inf;
  for I in eachindex(a)
    c = max(c,maximum(a[I]))
  end
  return c
end

function minimum(a::AbstractMeshArray)
  c=-Inf;
  for I in eachindex(a)
    c = min(c,minimum(a[I]))
  end
  return c
end

function sum(a::AbstractMeshArray)
  c=0.0
  for I in eachindex(a)
    c = c + sum(a[I])
  end
  return c
end

function fill(val::Any,a::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = fill(val,size(a.f[I]))
  end
  return c
end

function fill!(a::AbstractMeshArray,val::Any)
  for I in eachindex(a.f)
    fill!(a.f[I],val)
  end
  return a
end

import Base: +, -, *, /

function +(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] + b.f[I]
  end
  return c
end

function -(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] - b.f[I]
  end
  return c
end

function /(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] ./ b.f[I]
  end
  return c
end

function *(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] .* b.f[I]
  end
  return c
end


text:
using Test
using Mimi
using MimiDICE2010
using XLSX: readxlsx
using DataFrames
using CSVFiles

using MimiDICE2010: read_params, dice2010_excel_parameters

@testset "mimi-dice-2010" begin

    # ------------------------------------------------------------------------------
    #   1. Run the independent component tests
    # ------------------------------------------------------------------------------

    @testset "dice2010-components" begin

        include("test_climatedynamics.jl")
        include("test_co2cycle.jl")
        include("test_damages.jl")
        include("test_emissions.jl")
        include("test_grosseconomy.jl")
        include("test_neteconomy.jl")
        include("test_radiativeforcing.jl")
        include("test_sealevelrise.jl")
        include("test_welfare.jl")

    end # dice2010-components testset

    # ------------------------------------------------------------------------------
    #   2. Run tests on the whole model
    # ------------------------------------------------------------------------------

    @testset "dice2010-model" begin

        Precision = 1.0e-11
        T = 60
        f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

        m = MimiDICE2010.get_model()
        run(m)

        # Climate dynamics tests

        TATM    = m[:climatedynamics, :TATM]
        TOCEAN  = m[:climatedynamics, :TOCEAN]

        True_TATM   = read_params(f, "B121:BI121", T)
        True_TOCEAN = read_params(f, "B123:BI123", T)

        @test maximum(abs, TATM .- True_TATM) ≈ 0. atol = Precision
        @test maximum(abs, TOCEAN .- True_TOCEAN) ≈ 0. atol = Precision

        # CO2 Cycle tests

        MAT     = m[:co2cycle, :MAT]
        MAT_final   = m[:co2cycle, :MAT_final]
        ML      = m[:co2cycle, :ML]
        MU      = m[:co2cycle, :MU]

        True_MAT    = read_params(f, "B112:BI112", T)
        True_MAT_final  = read_params(f, "BJ112")
        True_ML     = read_params(f, "B115:BI115", T)
        True_MU     = read_params(f, "B114:BI114", T)

        @test maximum(abs, MAT .- True_MAT) ≈ 0. atol = Precision
        @test abs(MAT_final - True_MAT_final) ≈ 0. atol = Precision
        @test maximum(abs, ML .- True_ML) ≈ 0. atol = Precision
        @test maximum(abs, MU .- True_MU) ≈ 0. atol = Precision

        # Damages test

        DAMFRAC = m[:damages, :DAMFRAC]
        True_DAMFRAC    = read_params(f, "B93:BI93", T)
        @test maximum(abs, DAMFRAC .- True_DAMFRAC) ≈ 0. atol = Precision

        # Emissions tests

        CCA     = m[:emissions, :CCA]
        E       = m[:emissions, :E]
        EIND    = m[:emissions, :EIND]

        True_CCA    = read_params(f, "B117:BI117", T)
        True_E      = read_params(f, "B109:BI109", T)
        True_EIND   = read_params(f, "B110:BI110", T)

        @test maximum(abs, CCA .- True_CCA) ≈ 0. atol = Precision
        @test maximum(abs, E .- True_E) ≈ 0. atol = Precision
        @test maximum(abs, EIND .- True_EIND) ≈ 0. atol = Precision

        # Gross Economy tests

        K       = m[:grosseconomy, :K]
        YGROSS  = m[:grosseconomy, :YGROSS]

        True_K      = read_params(f, "B102:BI102", T)
        True_YGROSS = read_params(f, "B92:BI92", T)

        @test maximum(abs, K .- True_K) ≈ 0. atol = 3.0e-11 # Relax the precision just for this variable
        @test maximum(abs, YGROSS .- True_YGROSS) ≈ 0. atol = Precision

        # Net Economy tests

        ABATECOST   = m[:neteconomy, :ABATECOST]
        C           = m[:neteconomy, :C]
        CPC         = m[:neteconomy, :CPC]
        CPRICE      = m[:neteconomy, :CPRICE]
        I           = m[:neteconomy, :I]
        Y           = m[:neteconomy, :Y]
        YNET        = m[:neteconomy, :YNET]

        True_ABATECOST  = read_params(f, "B97:BI97", T)
        True_C          = read_params(f, "B125:BI125", T)
        True_CPC        = read_params(f, "B126:BI126", T)
        True_CPRICE     = read_params(f, "B134:BI134", T)
        True_I          = read_params(f, "B101:BI101", T)
        True_Y          = read_params(f, "B98:BI98", T)
        True_YNET       = read_params(f, "B95:BI95", T)

        @test maximum(abs, ABATECOST .- True_ABATECOST) ≈ 0. atol = Precision
        @test maximum(abs, C .- True_C) ≈ 0. atol = Precision
        @test maximum(abs, CPC .- True_CPC) ≈ 0. atol = Precision
        @test maximum(abs, CPRICE .- True_CPRICE) ≈ 0. atol = Precision
        @test maximum(abs, I .- True_I) ≈ 0. atol = Precision
        @test maximum(abs, Y .- True_Y) ≈ 0. atol = Precision
        @test maximum(abs, YNET .- True_YNET) ≈ 0. atol = Precision

        # Radiative Forcing test

        FORC = m[:radiativeforcing, :FORC]
        True_FORC    = read_params(f, "B122:BI122", T)
        @test maximum(abs, FORC .- True_FORC) ≈ 0. atol = Precision

        # Sea Level Rise tests

        ThermSLR    = m[:sealevelrise, :ThermSLR]
        GSICSLR     = m[:sealevelrise, :GSICSLR]
        GISSLR      = m[:sealevelrise, :GISSLR]
        AISSLR      = m[:sealevelrise, :AISSLR]
        TotSLR      = m[:sealevelrise, :TotSLR]

        True_ThermSLR    = read_params(f, "B178:BI178", T)
        True_GSICSLR    = read_params(f, "B179:BI179", T)
        True_GISSLR    = read_params(f, "B180:BI180", T)
        True_AISSLR    = read_params(f, "B181:BI181", T)
        True_TotSLR    = read_params(f, "B182:BI182", T)

        @test maximum(abs, ThermSLR .- True_ThermSLR) ≈ 0. atol = Precision
        @test maximum(abs, GSICSLR .- True_GSICSLR) ≈ 0. atol = Precision
        @test maximum(abs, GISSLR .- True_GISSLR) ≈ 0. atol = Precision
        @test maximum(abs, AISSLR .- True_AISSLR) ≈ 0. atol = Precision
        @test maximum(abs, TotSLR .- True_TotSLR) ≈ 0. atol = Precision

        # Welfare tests

        CEMUTOTPER  = m[:welfare, :CEMUTOTPER]
        PERIODU     = m[:welfare, :PERIODU]
        UTILITY     = m[:welfare, :UTILITY]

        True_CEMUTOTPER    = read_params(f, "B129:BI129", T)
        True_PERIODU    = read_params(f, "B128:BI128", T)
        True_UTILITY    = read_params(f, "B130")

        @test maximum(abs, CEMUTOTPER .- True_CEMUTOTPER) ≈ 0. atol = Precision
        @test maximum(abs, PERIODU .- True_PERIODU) ≈ 0. atol = Precision
        @test abs(UTILITY - True_UTILITY) ≈ 0. atol = Precision

    end # dice-2010 testset

    # ------------------------------------------------------------------------------
    #   3. Run tests to make sure integration version (Mimi v0.5.0)
    #   values match Mimi 0.4.0 values
    # ------------------------------------------------------------------------------

    @testset "dice2010-integration" begin

        Precision = 1.0e-11
        nullvalue = -999.999
        T = 60

        m = MimiDICE2010.get_model()
        run(m)

        for c in map(nameof, Mimi.compdefs(m)), v in Mimi.variable_names(m, c)

            # load data for comparison
            filepath = joinpath(@__DIR__, "../data/validation_data_v040/$c-$v.csv")
            results = m[c, v]

            df = load(filepath) |> DataFrame
            if typeof(results) <: Number
                validation_results = df[1,1]

            else
                validation_results = Matrix(df)

                # remove NaNs
                results[ismissing.(results)] .= nullvalue
                results[isnan.(results)] .= nullvalue
                validation_results[ismissing.(validation_results)] .= nullvalue
                validation_results[isnan.(validation_results)] .= nullvalue

                # match dimensions
                if size(validation_results, 1) == 1
                    validation_results = validation_results'
                end
            end
            @test results ≈ validation_results atol = Precision

        end # for loop
    end # dice2010-integration testset

    # ------------------------------------------------------------------------------
    #   4. Test standard api functions and SCC functions
    # ------------------------------------------------------------------------------

    @testset "Standard API functions" begin

        m = MimiDICE2010.get_model()
        run(m)

        # Test the errors
        @test_throws ErrorException MimiDICE2010.compute_scc()  # test that it errors if you don't specify a year
        @test_throws ErrorException MimiDICE2010.compute_scc(year=2020)  # test that it errors if the year isn't in the time index
        @test_throws ErrorException MimiDICE2010.compute_scc(last_year=2300)  # test that it errors if the last_year isn't in the time index
        @test_throws ErrorException MimiDICE2010.compute_scc(year=2105, last_year=2100)  # test that it errors if the year is after last_year

        # Test the SCC
        scc1 = MimiDICE2010.compute_scc(year=2015)
        @test scc1 isa Float64

        # Test that it's smaller with a shorter horizon
        scc2 = MimiDICE2010.compute_scc(year=2015, last_year=2295)
        @test scc2 < scc1

        # Test that it's smaller with a larger prtp
        scc3 = MimiDICE2010.compute_scc(year=2015, last_year=2295, prtp=0.02)
        @test scc3 < scc2

        # Test with a modified model
        m = MimiDICE2010.get_model()
        update_param!(m, :t2xco2, 5)
        scc4 = MimiDICE2010.compute_scc(m, year=2015)
        @test scc4 > scc1   # Test that a higher value of climate sensitivty makes the SCC bigger

        # Test compute_scc_mm
        result = MimiDICE2010.compute_scc_mm(year=2035)
        @test result.scc isa Float64
        @test result.mm isa Mimi.MarginalModel

    end

    # ------------------------------------------------------------------------------
    #   4. Test Deterministic SCC values
    # ------------------------------------------------------------------------------
    @testset "SCC values" begin

        atol = 1e-4 # TODO what is a reasonable tolerance given we test on a few different machines etc.

        # Test several validation configurations against the pre-saved values
        specs = Dict([
            :year => [2015, 2055],
            :eta => [0, 1.5],
            :prtp => [0.015, 0.03],
            :last_year => [2295, 2595],
            :pulse_size => [1e3, 1e7, 1e10]
        ])

        results = DataFrame(year = [], eta = [], prtp = [], last_year = [], pulse_size = [], SC = [])

        for year in specs[:year]
            for eta in specs[:eta]
                for prtp in specs[:prtp]
                    for last_year in specs[:last_year]
                        for pulse_size in specs[:pulse_size]
                            sc = MimiDICE2010.compute_scc(year=Int(year), eta=eta, prtp=prtp, last_year=Int(last_year), pulse_size=pulse_size)
                            push!(results, (year, eta, prtp, last_year, pulse_size, sc))
                        end
                    end
                end
            end
        end

        validation_results = load(joinpath(@__DIR__, "..", "data", "SC validation data", "deterministic_sc_values_v1-0-1.csv")) |> DataFrame
        # diffs = sort(results[!, :SC] - validation_results[!, :SC], rev = true)
        # println(diffs)
        @test all(isapprox.(results[!, :SC], validation_results[!, :SC], atol = atol))

    end # SCC values testset

end # mimi-dice-2010 testset

nothing


text:
abstract type EpidemicModel end

abstract type SEIR <: EpidemicModel end
abstract type SIR <: EpidemicModel end
abstract type SEI <: EpidemicModel end
abstract type SI <: EpidemicModel end


text:
#################################
# Markov chains for game theory #
#################################

"""
    game2markov(game::Dict{String,<:Array}, s::Vector{Vector{T}}) where T<:Real

Returns Markov chain of a given game and given vector of action probabilities

**Input parameters**
* `game::Dict{String,<:Array}` - dictionary of players and their payoff matrices
* `s::Vector{Vector{T}}` - collection of actions probabilities for each player
"""
function game2markov(game::Dict{String,<:Array},
                    s::Vector{Vector{T}}) where T<:Real
    # problem: 2 players
    permuts = permutations( vcat(1, fill(0, size(s[1])[1]-1 )) ) |> collect |> unique
    grid_exp = Iterators.product(collect(permuts),collect(permuts)) |> collect
    to_states = Dict(string(i) => collect(grid_exp[i]) for i in 1:length(grid_exp))
    to_states_inv = Dict(collect(grid_exp[i]) => i for i in 1:length(grid_exp))

    trans = Matrix(0I, length(to_states), length(to_states))
    for i in 1:length(to_states)
        indx = Nash._next_state(game, to_states[string(i)], to_states_inv)
        trans[i, indx] = 1/length(indx) # distributed equally when more than 1 best reply
    end

    init = Array(1:length(to_states))
    mc = MarkovChain(trans, init)
    return mc
end

"""
    plot_markov(no_steps::Int, mc::MarkovChain{Int,Array{Int,2},Array{Int,1}})

Plots a path of Markov chain simulation of a given time length and all initial values

**Input parameters**
* `no_steps::Int` - scalar of time length (number of iterations)
* `mc::MarkovChain{Int,Array{Int,2},Array{Int,1}}` - Markov chain of a game (game2markov output function)
"""
function plot_markov(no_steps::Int,
        mc::MarkovChain{Int,Array{Int,2},Array{Int,1}})
    no_states = stationary_distributions(mc)[1] |> length
    for i in 1:no_states
        line = simulate(mc, no_steps, init = i)
        i == 1 ? display(plot(1:no_steps, line, label = string(i), lw = 2)) : display(plot!(line, label = string(i), lw = 2))
    end
end


text:
using Gadfly
using Distributions

f(u) = exp(-u^2/2)/√(2pi)

const curve_plot =
    plot(f, -5, +5,
        Theme(line_width=2Gadfly.mm,
              major_label_font_size=40Gadfly.px,
              minor_label_font_size=25Gadfly.px)
    )

const u1=Uniform(-5.0,+5.0)
const u2=Uniform(0.0, 0.5)

function simulate_pt()
    x=rand(u1)
    y=rand(u2)
    y<f(x)
end

function simulate(num, running, current_approx)
     push!(running, true)
     hits=0
     for (i in 1:num)
        hits = hits + simulate_pt()
        if (i % 10000) == 0
            push!(current_approx, (i, hits/i*(0.5*10) ))
        end
    end
    estimate = hits/num*(0.5*10)
    push!(running, false)
end

function run_simulate(N, running, current_approx)
    println("Starting Simulations for N=$N")
    @async simulate(N, running, current_approx)
end

function main(window)
    push!(window.assets, "tex")
    push!(window.assets, "widgets")

    N = Signal(10000)
    btn = Signal(Any, leftbutton)

    current_approx = Signal((0, 1.0)) # current approximation
    buffered_approx = foldp((Any[0],Any[1.0]), current_approx) do prev, current
        idx, val = current
        push!(prev[1], idx), push!(prev[2], val)
    end

    running = Signal(false)
    result = map(sampleon(btn, N); typ=Any, init=nothing) do n
        run_simulate(n, running, current_approx)
    end

    vbox(
        h1("Interactive Simulations"),
        hbox(
            vbox(
                 md"""We want to estimate the following integral using a
                 *Monte Carlo* Simulation""",
                 tex("\\int_{-5}^{5} \\frac{exp(-u^2/2)}{\\sqrt{2pi}}du" , block=true),
            ),
            curve_plot
         ) |> packacross(center),
         map(buffered_approx, running) do approx, r
             xs, ys = approx
             vbox(
                 hbox("Number of runs", slider(10^3:10^3:5*10^6) >>> N, hskip(2em), hbox(button("Start", raised=true, disabled=r) >>> btn)) |> packacross(center),
                 hbox("Current value: ", hskip(1em), @sprintf("%2.4f", ys[end]) |> emph, hskip(1em)," at iteration :", hskip(1em), string(xs[end])),
                 plot(x=xs, y=ys, Geom.line) |> Escher.pad(2em)
             )
         end
    ) |> Escher.pad(2em)
end


text:
# Julia wrapper for header: /usr/include/scip/cons_countsols.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrCountsols(scip)
    ccall((:SCIPincludeConshdlrCountsols, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPdialogExecCountPresolve(scip, dialog, dialoghdlr, nextdialog)
    ccall((:SCIPdialogExecCountPresolve, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}, Ptr{SCIP_DIALOGHDLR}, Ptr{Ptr{SCIP_DIALOG}}), scip, dialog, dialoghdlr, nextdialog)
end

function SCIPdialogExecCount(scip, dialog, dialoghdlr, nextdialog)
    ccall((:SCIPdialogExecCount, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}, Ptr{SCIP_DIALOGHDLR}, Ptr{Ptr{SCIP_DIALOG}}), scip, dialog, dialoghdlr, nextdialog)
end

function SCIPdialogExecWriteAllsolutions(scip, dialog, dialoghdlr, nextdialog)
    ccall((:SCIPdialogExecWriteAllsolutions, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}, Ptr{SCIP_DIALOGHDLR}, Ptr{Ptr{SCIP_DIALOG}}), scip, dialog, dialoghdlr, nextdialog)
end

function SCIPcount(scip)
    ccall((:SCIPcount, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPgetNCountedSols(scip, valid)
    ccall((:SCIPgetNCountedSols, libscip), Clonglong, (Ptr{SCIP_}, Ptr{UInt32}), scip, valid)
end

function SCIPgetNCountedSolsstr(scip, buffer, buffersize, requiredsize)
    ccall((:SCIPgetNCountedSolsstr, libscip), Cvoid, (Ptr{SCIP_}, Ptr{Cstring}, Cint, Ptr{Cint}), scip, buffer, buffersize, requiredsize)
end

function SCIPgetNCountedFeasSubtrees(scip)
    ccall((:SCIPgetNCountedFeasSubtrees, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetCountedSparseSols(scip, vars, nvars, sols, nsols)
    ccall((:SCIPgetCountedSparseSols, libscip), Cvoid, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Cint}, Ptr{Ptr{Ptr{SCIP_SPARSESOL}}}, Ptr{Cint}), scip, vars, nvars, sols, nsols)
end

function SCIPsetParamsCountsols(scip)
    ccall((:SCIPsetParamsCountsols, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end


text:
using Interpolations

function hdr_sequence_downsample{T<:Integer}(hdrseq::HDRSequence, downsample_ratio::T)
    kernel_size = (0.75 * 1, 0.75 * 1)
    gaussian = KernelFactors.gaussian(kernel_size)

    quick_downsample_subm = [hdrseq[i][1][1:downsample_ratio÷2:size(hdrseq[i][1], 1),
                                          1:downsample_ratio÷2:size(hdrseq[i][1], 2)] for i in 1:size(hdrseq, 1)]
    filtered_image = [imfilter(RGB{N0f8}, quick_downsample_subm[i], gaussian, NA()) for i in 1:size(hdrseq, 1)]
    downsampled_seq = [(filtered_image[i], hdrseq[i][2]) for i in 1:size(hdrseq, 1)]

    return downsampled_seq
end

function inpaint_nans(v::Vector)
    y = v[!isnan(v)]
    x = find(x -> x==false, isnan(v))
    knots_x = (x, )
    interpolator = interpolate(knots_x, y, Gridded(Linear()))

    v_new = copy(v)
    for nan_idx in find(x -> x==true, isnan(v))
        v_new[nan_idx] = interpolator[nan_idx]
    end
    return v_new
end

function image_synthesis(hdr_pairs,
                         iterations = -1,
                         show_intermediate_results = false)
# description: generating a HDR image by combining multiple images of different
#              exposure time, while recovering camera's response curve by Gauss-
#              Seidel relaxation.
# input: image sequence as parsed by read_hdrgen_file
# output: one single HDR image and response function
    fw = (x) -> (x<2/255 || x>253/255) ? 0 : exp(-5 * (x - 0.5)^2 / (0.5^2)) # weighing function

    image_size = size(channelview(hdr_pairs[1][1]))
    color_depth = Integer(typemax(channelview(hdr_pairs[1][1])[1].i) + 1)
    color_channel = Integer(size(channelview(hdr_pairs[1][1]), 1))

    fi = Array(Float64, (color_depth, color_channel)) # response curve in R, G, B channel
    for i = 1:color_channel
        fi[:, i] = linspace(0, 2, color_depth)
    end

    last_ir = zeros(Float64, image_size)
    last_res = zeros(Float64, (color_depth, color_channel))

    iteration_count = 0

    while (iteration_count < iterations) || iterations == -1
        ir = zeros(Float64, image_size) #irradiance
        ir_normalizer = zeros(Float64, image_size)

        new_fi = zeros(Float64, (color_depth, color_channel))
        new_fi_normalizer = zeros(Float64, (color_depth, color_channel))

        for (image, exposure) in hdr_pairs
            image_ch = channelview(image)
            for j in 1:size(image_ch, 3), i in 1:size(image_ch, 2), k in 1:size(image_ch, 1)
                ch_i = image_ch[k, i, j] # channel intensity

                ir[k, i, j] += fw(ch_i) * exposure * fi[ch_i.i + 1]
                ir_normalizer[k, i, j] += fw(ch_i) * exposure^2
            end
        end

        ir = ir ./ ir_normalizer
        #ir[isnan(ir)] = 0 # silent NaNs

        for (image, exposure) in hdr_pairs
            image_ch = channelview(image)
            for j in 1:size(image_ch, 3), i in 1:size(image_ch, 2), k in 1:size(image_ch, 1)
                ch_i = image_ch[k, i, j]

                new_fi_normalizer[ch_i.i + 1, k] += 1
                new_fi[ch_i.i + 1, k] += exposure * ir[k, i, j]
            end
        end

        # normalizing & constraining fi(mid intensity) = 1.0
        fi = new_fi ./ new_fi_normalizer
        for ch in 1:color_channel
            fi[:, ch] = inpaint_nans(fi[:, ch]) # extrapolate NaN holes
            fi[:, ch] = fi[:, ch] / fi[color_depth ÷ 2, ch]
            #fi[:, ch] = fi_avg[:] / fi_avg[color_depth ÷ 2]

            #for intensity in 1:color_depth-1
            #    if fi[intensity, ch] > fi[intensity+1, ch]
            #        fi[intensity+1, ch] = fi[intensity, ch]
            #    end
            #end
        end

        if show_intermediate_results
            imgc, imsl = imshow(colorview(RGB, ir) / maximum(ir)) # show by simple linear mapping
            #imgc, imsl = imshow(ir / maximum(ir)) # show gray
            ImageView.annotate!(imgc, imsl,
                                ImageView.AnnotationText(40, 30, "i=$iteration_count", color = RGB(1, 0, 0),
                                fontsize = 15))
        end

        res_diff = fi - last_res
        res_diff_norm = [norm(res_diff[:, k]) for k in 1:size(res_diff, 2)]
        #println(norm(res_diff_norm)) #debug
        #println(last_res[1:10, 1])

        last_ir = copy(ir)
        last_res = copy(fi)

        iteration_count += 1
        if iterations == -1 && norm(res_diff_norm) < 2.0
            break
        end

   end #for iterations

   return (last_ir, last_res)
end

function image_synthesis_hsl(hdr_pairs,
                             iterations = -1,
                             show_intermediate_results = false)
    color_channel = size(channelview(hdr_pairs[1][1]), 1)
    num_images = size(hdr_pairs, 1)
    image_dimension = size(channelview(hdr_pairs[1][1]))

    rgb_sum = zeros(Float64, image_dimension) # TODO: add weighing func?
    for (image, exposure) in hdr_pairs
        rgb_sum += channelview(image)
    end
    rgb_avg = colorview(RGB, rgb_sum / num_images)
    hsl_avg = HSL.(rgb_avg) # we are going to use H & S channel

    hdr_light_pairs = [(N0f8.(channelview(HSL.(image))[3, :, :]), exposure)
                        for (image, exposure) in hdr_pairs]

    ir, res = image_synthesis(hdr_light_pairs, iterations, show_intermediate_results)

    # quick exponential mapping
    ir_avg = mean(ir)
    fr = (x) -> 1 - exp(-x / ir_avg)
    ir_remapped = fr.(ir)

    hsl_channel = channelview(hsl_avg)
    hsl_channel[3, :, :] = ir_remapped

    return (channelview(RGB.(colorview(HSL, hsl_channel))), res)
end


text:
module Margin

import ..Clarus

function attribution(;params...)
  return Clarus.api_request("Margin","Attribution";params...)
end

function im(;params...)
  return Clarus.api_request("Margin","IM";params...)
end

function impact(;params...)
  return Clarus.api_request("Margin","Impact";params...)
end

function optimiser(;params...)
  return Clarus.api_request("Margin","Optimiser";params...)
end

function parswapsim(;params...)
  return Clarus.api_request("Margin","ParSwapsIM";params...)
end

function plvectors(;params...)
  return Clarus.api_request("Margin","PLVectors";params...)
end

function sensitivity(;params...)
  return Clarus.api_request("Margin","Sensitivity";params...)
end

function vm(;params...)
  return Clarus.api_request("Margin","VM";params...)
end

function vmbuffer(;params...)
  return Clarus.api_request("Margin","VMBuffer";params...)
end

function vmlsoc(;params...)
  return Clarus.api_request("Margin","VMLSOC";params...)
end

end #MODULE-END


text:
function sample!(fin_el::Quadrilateral, s::Matrix{Float64} , wt::Vector{Float64})
  #
  # This subroutine returns the local coordinates and weighting coefficients
  # of the integrating points.
  #

  root3 = 1.0 / sqrt(3.0)
  r15 = 0.2 * sqrt(15.0)
  nip = size(s,1)

  w = [5.0/9.0, 8.0/9.0, 5.0/9.0]
  v = [5.0/9.0*w; 8.0/9.0*w; 5.0/9.0*w]

  if nip == 1
    s[1,1] = 0.0
    s[1,2] = 0.0
    wt[1]  = 4.0
  elseif nip == 4
    s[1,1] = -root3
    s[1,2] = root3
    s[2,1] = root3
    s[2,2] = root3
    s[3,1] = -root3
    s[3,2] = -root3
    s[4,1] = root3
    s[4,2] = -root3
    fill!(wt, 1)
  elseif nip == 9
    s[1:3:7,1] .= -r15
    s[2:3:8,2] .= 0.0
    s[3:3:9,1] .= r15
    s[1:3,2]   .= r15
    s[4:6,2]   .= 0.0
    s[7:9,2]   .= -r15
    wt[:] = v
  elseif nip == 16
    s[1:4:13,1] .= -0.861136311594053
    s[2:4:14,1] .= -0.339981043584856
    s[3:4:15,1] .= 0.339981043584856
    s[4:4:16,1] .= 0.861136311594053
    s[1:4,2]    .= 0.861136311594053
    s[5:8,2]    .= 0.339981043584856
    s[9:12,2]   .= -0.339981043584856
    s[13:16,2]  .= -0.861136311594053
    wt[1]       = 0.121002993285602
    wt[4]       = wt[1]
    wt[13]      = wt[1]
    wt[16]      = wt[1]
    wt[2]       = 0.226851851851852
    wt[3]       = wt[2]
    wt[5]       = wt[2]
    wt[8]       = wt[2]
    wt[9]       = wt[2]
    wt[12]      = wt[2]
    wt[14]      = wt[2]
    wt[15]      = wt[2]
    wt[6]       = 0.425293303010694
    wt[7]       = wt[6]
    wt[10]      = wt[6]
    wt[11]      = wt[6]
  elseif nip == 25
    s[1:5:21,1] = 0.906179845938664
    s[2:5:22,1] = 0.538469310105683
    s[3:5:23,1] = 0.0
    s[4:5:24,1] = -0.538469310105683
    s[5:5:25,1] = -0.906179845938664
    s[ 1: 5,2]  = 0.906179845938664
    s[ 6:10,2]  = 0.538469310105683
    s[11:15,2]  = 0.0
    s[16:20,2]  = -0.538469310105683
    s[21:25,2]  = -0.906179845938664
    wt[1]  = 0.056134348862429
    wt[2]  = 0.113400000000000
    wt[3]  = 0.134785072387521
    wt[4]  = 0.113400000000000
    wt[5]  = 0.056134348862429
    wt[6]  = 0.113400000000000
    wt[7]  = 0.229085404223991
    wt[8]  = 0.272286532550750
    wt[9]  = 0.229085404223991
    wt[10] = 0.113400000000000
    wt[11] = 0.134785072387521
    wt[12] = 0.272286532550750
    wt[13] = 0.323634567901235
    wt[14] = 0.272286532550750
    wt[15] = 0.134785072387521
    wt[16] = 0.113400000000000
    wt[17] = 0.229085404223991
    wt[18] = 0.272286532550750
    wt[19] = 0.229085404223991
    wt[20] = 0.113400000000000
    wt[21] = 0.056134348862429
    wt[22] = 0.113400000000000
    wt[23] = 0.134785072387521
    wt[24] = 0.113400000000000
    wt[25] = 0.056134348862429
  else
    println("Wrong number of integrating points for a quadrilateral.")
  end
end


text:
using Documenter, GeneNameGenieJ

# makedocs()
makedocs(
   modules = [GeneNameGenieJ],
   clean = false,
   format = :html,
   sitename = "GeneNameGenieJ.jl",
   doctest = false,
   pages = Any[
      "Introduction" => "index.md",
      "User Guide" => Any[
          "Getting Started" => "man/getting-started.md",
          "miRNAs" => "man/mirna-functions.md",
          "Molecular-ID handling" => "man/molecular-id-handling.md",
          "Attributes" => "man/attribute-functions.md"
      ],
      "API" => Any[
          "Functions" => "lib/functions.md"
      ]
    ]
)

# Deploy built documentation from Travis.
# =======================================

# deploydocs(
#     # options
#     repo = "github.com/JuliaData/DataFrames.jl.git",
#     target = "build",
#     julia = "0.6",
#     deps = nothing,
#     make = nothing,
# )

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


text:
"""
    FactorGraph{T<:Integer} <: AbstractIndexedGraph{T}

A type representing a sparse factor graph.

### FIELDS

- `A` -- adjacency matrix filled with `NullNumber`s. Rows are factors, columns are variables
- `X` -- square matrix for efficient access by row. `X[j,i]` points to the index of element `A[i,j]` in `A.nzval`.
"""
struct FactorGraph{T<:Integer} <: AbstractIndexedGraph{T}
    A :: SparseMatrixCSC{NullNumber, T}
    X :: SparseMatrixCSC{T, T}
end

function FactorGraph(A::AbstractMatrix{NullNumber})
    A = sparse(A)
    X = sparse(SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, collect(1:nnz(A)))')
    FactorGraph(A, X)
end

"""
    FactorGraph(A::AbstractMatrix)

Construct a factor graph from adjacency matrix `A` with the convention that
rows are factors, columns are variables
"""
function FactorGraph(A::AbstractMatrix)
    A = sparse(A)
	FactorGraph(SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, fill(NullNumber(), length(A.nzval))))
end

abstract type VariableOrFactor{T<:Integer}; end

"""
    Variable{T<:Integer}

Wraps an index to specify that it is the index of a variable.
See e.g. [`neighbors(g::FactorGraph, v::Variable)`](@ref)
"""
struct Variable{T<:Integer} <: VariableOrFactor{T}
    i::T
    function Variable(i::T) where {T<:Integer}
        i > 0 || throw(ArgumentError("Variable index must be positive, got $i"))
        new{T}(i)
    end
end

"""
    Factor{T<:Integer}

Wraps an index to specify that it is the index of a factor.
See e.g. [`neighbors(g::FactorGraph, f::Factor)`](@ref)
"""
struct Factor{T<:Integer} <: VariableOrFactor{T}
    a::T
    function Factor(a::T) where {T<:Integer}
        a > 0 || throw(ArgumentError("Factor index must be positive, got $a"))
        new{T}(a)
    end
end

"""
FactorGraphEdge{T<:Integer} <: AbstractIndexedEdge{T}

Edge type for `FactorGraph`s.
### Safe constructors
    FactorGraphEdge(f::Factor, v::Variable, idx::Integer)
    FactorGraphEdge(v::Variable, f::Factor, idx::Integer)

Example:
```@example
    FactorGraphEdge(Factor(2), Variable(4), 3)
```
"""
struct FactorGraphEdge{T<:Integer} <: AbstractIndexedEdge{T}
	a::T    # factor index
	i::T    # variable index
    idx::T  # edge index
    function FactorGraphEdge(a::T, i::T, idx::T) where {T<:Integer}
        i > 0 || throw(ArgumentError("Variable index must be positive, got $i"))
        a > 0 || throw(ArgumentError("Factor index must be positive, got $a"))
        idx > 0 || throw(ArgumentError("Edge index must be positive, got $idx"))
        new{T}(a, i, idx)
    end
end

function FactorGraphEdge(a::Integer, i::Integer, idx::Integer)
    a_, i_, idx_ = promote(a, i, idx)
    FactorGraphEdge(a_, i_, idx_)
end

function FactorGraphEdge(f::Factor{T}, v::Variable{T}, idx::T) where {T<:Integer}
    FactorGraphEdge(f.a, v.i, idx)
end
function FactorGraphEdge(v::Variable{T}, f::Factor{T}, idx::T) where {T<:Integer}
    FactorGraphEdge(f, v, idx)
end

function Base.show(io::IO, e::FactorGraphEdge)
    print(io, "Edge from factor $(e.a) to variable $(e.i) with index $(e.idx)")
end

Base.eltype(g::FactorGraph{T}) where T = T
edgetype(g::FactorGraph{T}) where T = FactorGraphEdge{T}

nvariables(g::FactorGraph) = size(g.A, 2)
nfactors(g::FactorGraph) = size(g.A, 1)

factors(g::FactorGraph) = (Factor(a) for a in 1:nfactors(g))
variables(g::FactorGraph) = (Variable(i) for i in 1:nvariables(g))

ne(g::FactorGraph) = nnz(g.A)
nv(g::FactorGraph) = nvariables(g) + nfactors(g)

edges(g::FactorGraph) = (FactorGraphEdge(g.A.rowval[k], j, k) for j=1:size(g.A,2) for k=nzrange(g.A,j))
# vertices(g::FactorGraph) = Iterators.flatten((factors(g), variables(g)))
vertices(g::FactorGraph) = 1:nv(g)

"""
    edges(g::FactorGraph, v::Variable)

Return a lazy iterator to the edges incident to variable `v`
"""
function edges(g::FactorGraph, v::Variable)
    i = v.i
    (FactorGraphEdge( g.A.rowval[k], i, k) for k in nzrange(g.A, i))
end

"""
    neighbors(g::FactorGraph, f::Factor)

Return a lazy iterator to the edges incident to factor `f`
"""
function edges(g::FactorGraph, f::Factor)
    a = f.a
    (FactorGraphEdge(a, g.X.rowval[k], k) for k in nzrange(g.X, a))
end
inedges(g::FactorGraph, x::VariableOrFactor) = edges(g, x)
outedges(g::FactorGraph, x::VariableOrFactor) = edges(g, x)

"""
    neighbors(g::FactorGraph, v::Variable)

Return a lazy iterator to the neighbors of variable `v`
"""
neighbors(g::FactorGraph, v::Variable) = @view g.A.rowval[nzrange(g.A, v.i)]

"""
    neighbors(g::FactorGraph, f::Factor)

Return a lazy iterator to the neighbors of factor `f`
"""
neighbors(g::FactorGraph, f::Factor) = @view g.X.rowval[nzrange(g.X, f.a)]
inneighbors(g::FactorGraph, x::VariableOrFactor) = neighbors(g, x)
outneighbors(g::FactorGraph, x::VariableOrFactor) = neighbors(g, x)

function neighbors(g::FactorGraph, i::Int)
    nvars = nvariables(g)
    if i > nvars
        # i is a factor
        return (j + nvars for j in neighbors(g, Factor(i-nvars)))
    else
        # i is a variable
        return neighbors(g, Variable(i))
    end
end

has_vertex(g::FactorGraph, v::Variable) = v.i ≤ nvariables(g)
has_vertex(g::FactorGraph, f::Factor) = f.a ≤ nfactors(g)
has_vertex(g::FactorGraph, i::Int) = i ≤ nv(g)
has_edge(g::FactorGraph, f::Factor, v::Variable) = g.A[f.a, v.i] != 0
has_edge(g::FactorGraph, v::Variable, f::Factor) = has_edge(g, f, v)

# return true if (src, dst) is an edge in the bipartite view of g
function has_edge(g::FactorGraph, src::Integer, dst::Integer)
    has_edge(g, Factor(src - nvariables(g)), Variable(dst))
end

Base.zero(g::FactorGraph) = FactorGraph(zero(g.A))

is_directed(g::FactorGraph) = false
is_directed(::Type{FactorGraph{T}}) where T = false
is_bipartite(g::FactorGraph) = true

edge_idx(g::FactorGraph, f::Factor, v::Variable) = nzindex(g.A, f.a, v.i)
edge_idx(g::FactorGraph, v::Variable, f::Factor) = edge_idx(g, f, v)
edge_src_dst(g::FactorGraph, id::Integer) = nzindex(g.A, id)

"""
    get_edge(g::FactorGraph, v::Variable, f::Factor)
    get_edge(g::FactorGraph, f::Factor, v::Variable)
    get_edge(g::FactorGraph, id::Integer)

Get edge given variable and factor or given edge index
"""
function get_edge(g::FactorGraph, x::VariableOrFactor, y::VariableOrFactor)
    k = edge_idx(g, x, y)
    FactorGraphEdge(x, y, k)
end

function get_edge(g::FactorGraph, id::Integer)
    a, i = edge_src_dst(g, id)
    FactorGraphEdge(a, i, id)
end

"""
    adjacency_matrix(g::FactorGraph, T::DataType=Int)

Return the symmetric adjacency matrix of size `nvariables(g) + nfactors(g)`
where no distinction is made between variable and factor nodes.
Edge indices are preserved in the bottom-left block.
"""
function adjacency_matrix(g::FactorGraph, T::DataType=Int)
    m, n = nfactors(g), nvariables(g)
    A = SparseMatrixCSC(g.A.m, g.A.n, g.A.colptr, g.A.rowval, ones(T, nnz(g.A)))
    return [ spzeros(T, n, n)  A'               ;
             A                 spzeros(T, m, m) ]
end

"""
    bipartite_view(g::FactorGraph, T::DataType=Int)

Construct an undirected bipartite `IndexedGraph` where vertices are the concatenation
of the variables and factors of the original `FactorGraph`.
Edge indices are preserved in the bottom-left block.
"""
function bipartite_view(g::FactorGraph, T::DataType=Int)
    A = adjacency_matrix(g, T)
    IndexedGraph(A)
end