text:
using PyPlot
using Random, Images, JLD2, LinearAlgebra
using JOLI, Statistics, FFTW
using Printf
using JUDI4Cloud

using ArgParse
include("../utils/parse_cmd.jl")
parsed_args = parse_commandline()
L = parsed_args["nv"]
vm = parsed_args["vm"]
niter = parsed_args["niter"]
nth = parsed_args["nth"]
nsrc = parsed_args["nsrc"]
batchsize = parsed_args["bs"]

Random.seed!(1234);

creds=joinpath(pwd(),"..","credentials.json")
init_culsterless(batchsize*L; credentials=creds, vm_size=vm, pool_name="JRM", verbose=1, nthreads=nth, auto_scale=false)

JLD2.@load "../models/Compass_tti_625m.jld2"
JLD2.@load "../models/timelapsevrho$(L)vint.jld2" vp_stack rho_stack
JLD2.@load "../data/dobs$(L)vint$(nsrc)nsrc.jld2" dobs_stack q_stack

idx_wb = find_water_bottom(rho.-rho[1,1])

nsrc = dobs_stack[1].nsrc

v0_stack = deepcopy(vp_stack./1f3)
for i = 1:L
    v0_stack[i][:,maximum(idx_wb)+1:end] .= 1f0./convert(Array{Float32,2},imfilter(1f3./vp_stack[i][:,maximum(idx_wb)+1:end], Kernel.gaussian(10)))
end
m0_stack = [1f0./v0_stack[i].^2f0 for i = 1:L]

model0_stack = [Model(n,d,o,m0_stack[i]; nb=80) for i = 1:L]

opt = JUDI.Options(isic=true)

# Preconditioner

Tm = judiTopmute(model0_stack[1].n, maximum(idx_wb), 3)  # Mute water column
S = judiDepthScaling(model0_stack[1])
Mr = Tm*S

# Linearized Bregman parameters
nn = prod(model0_stack[1].n)
x = zeros(Float32, nn)
z = zeros(Float32, nn)

fval = zeros(Float32, niter)

# Soft thresholding functions and Curvelet transform
soft_thresholding(x::Array{Float64}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float64, lambda), 0.0)
soft_thresholding(x::Array{Float32}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float32, lambda), 0f0)

C0 = joCurvelet2D(2*n[1], 2*n[2]; zero_finest = false, DDT = Float32, RDT = Float64)

function C_fwd(im, C, n)
	im = hcat(reshape(im, n), reshape(im, n)[:, end:-1:1])
    im = vcat(im, im[end:-1:1,:])
	coeffs = C*vec(im)/2f0
	return coeffs
end

function C_adj(coeffs, C, n)
	im = reshape(C'*coeffs, 2*n[1], 2*n[2])
	return vec(im[1:n[1], 1:n[2]] + im[1:n[1], end:-1:n[2]+1] + im[end:-1:n[1]+1, 1:n[2]] + im[end:-1:n[1]+1, end:-1:n[2]+1])/2f0
end

C = joLinearFunctionFwd_T(size(C0, 1), n[1]*n[2],
                          x -> C_fwd(x, C0, n),
                          b -> C_adj(b, C0, n),
                          Float32,Float64, name="Cmirrorext")

src_list = [collect(1:nsrc) for i = 1:L]

ps = 0
γ  = L/2f0 # hyperparameter to tune

x = [zeros(Float32, size(C,2)) for i=1:L+1];
z = [zeros(Float64, size(C,1)) for i=1:L+1];

tau = [zeros(Float32,size(C,1)) for i=1:L+1];
flag = [BitArray(undef,size(C,1)) for i=1:L+1];
sumsign = [zeros(Float32,size(C,1)) for i=1:L+1]

lambda = zeros(Float64,L+1)

# Main loop

for  j=1:niter

	# Main loop
    @printf("JRM Iteration: %d \n", j)
    flush(Base.stdout)

    xdm = [1f0/γ*x[1]+x[i+1] for i=1:L];
    global inds = [zeros(Int, batchsize) for i=1:L]

    for i = 1:L
        length(src_list[i]) < batchsize && (global src_list[i] = collect(1:nsrc))
		src_list[i] = src_list[i][randperm(MersenneTwister(i+2000*j),length(src_list[i]))]
		global inds[i] = [pop!(src_list[i]) for b=1:batchsize]
		println("Vintage $(i) Imaging source $(inds[i])")
    end

    source = [q_stack[i][inds[i]] for i = 1:L]
    dObs = [dobs_stack[i][inds[i]] for i = 1:L]
    dmx = [Mr*xdm[i] for i = 1:L]

    phi, g1 = lsrtm_objective(model0_stack, source, dObs, dmx; options=opt, nlind=false)
    g = [1f0/γ*C*Mr'*vec(sum(g1)), [C*Mr'*vec(g1[i]) for i=1:L]...]

	@printf("At iteration %d function value is %2.2e \n", j, phi)
	flush(Base.stdout)
	# Step size and update variable

	global t = γ^2f0*2f-5/(L+γ^2f0)/12.5 # fixed step

    # anti-chatter
	for i = 1:L+1
		global sumsign[i] = sumsign[i] + sign.(g[i])
		global tau[i] .= t
		global tau[i][findall(flag[i])] = deepcopy((t*abs.(sumsign[i])/j)[findall(flag[i])])
		global z[i] -= tau[i] .* g[i]
	end

	(j==1) && global lambda = [quantile(abs.(vec(z[i])), .9) for i = 1:L+1]	# estimate thresholding parameter at 1st iteration
    lambda1 = maximum(lambda[2:end])
    for i = 2:L+1
        global lambda[i] = lambda1
    end

	# Update variables and save snapshot
	for i = 1:L+1
		global x[i] = adjoint(C)*soft_thresholding(z[i], lambda[i])
		global flag[i] = flag[i] .| (abs.(z[i]).>=lambda[i])     # check if ever pass the threshold
	end
    JLD2.@save "../results/JRM$(j)Iter$(L)vintages$(nsrc)nsrc.jld2" x z g lambda phi
end


text:
using Dates

"""
Runs a Monte Carlo simulation with the FUND model over it's distirbutional parameters, and calculates
the Social Cost of Carbon for the specified `years` and discount `rates`.
`trials`: the number of trials to run.
`ntimesteps`: how many timesteps to run.
`output_dir`: an output directory; if none provided, will create and use "output/SCC yyyy-mm-dd HH-MM-SS MCtrials".
`save_trials`: whether or not to generate and save the MC trial values up front to a file.
"""
function run_fund_scc_mcs(trials = 10000; years = [2020], rates = [0.03], ntimesteps = MimiFUND.default_nsteps + 1, output_dir = nothing, save_trials = false)

    # Set up output directories
    output_dir = output_dir === nothing ? joinpath(@__DIR__, "../../output/", "SCC $(Dates.format(now(), "yyyy-mm-dd HH-MM-SS")) MC$trials") : output_dir
    mkpath("$output_dir/results")
    save_trials ? trials_output_filename = joinpath(@__DIR__, "$output_dir/trials.csv") : trials_output_filename = nothing

    # Write header in SCC output file
    scc_file = joinpath(output_dir, "scc.csv")
    open(scc_file, "w") do f
        write(f, "trial, rate, year, scc\n")
    end

    # Scenario set up
    scenario_args = Any[
        :rate           => rates
        :emissionyear   => years
    ]

    # get models and sim
    mcs = getmcs()
    mm = get_marginal_model()

    # Define scenario function
    function _scenario_func(mcs::SimulationInstance, tup::Tuple)

        # Unpack scenario tuple argument
        (rate, emissionyear) = tup

        # Get models
        mm = mcs.models[1]

        # Perturb emissons in the marginal model
        perturb_marginal_emissions!(mm.marginal, emissionyear)

    end

    # Define post trial function
    function scc_calculation(mcs::SimulationInstance, trialnum::Int, ntimesteps::Int, tup::Tuple)

        # Unpack scenario tuple argument
        (rate, emissionyear) = tup

        # Get marginal damages
        mm = mcs.models[1]
        marginaldamages = mm[:impactaggregation, :loss]

        # Calculate discount factor
        T = ntimesteps == typemax(Int) ? length(Mimi.dimension(mm.base, :time)) : ntimesteps
        discount_factor = zeros(T)
        idx = getindexfromyear(emissionyear)
        discount_factor[idx:T] = [1 / ((1 + rate) ^ t) for t in 0:T-idx]

        # Sum discounted global damages to scc
        scc = sum(sum(marginaldamages, dims = 2)[2:T] .* discount_factor[2:end])

        # Write output
        open(scc_file, "a") do f
            write(f, "$trialnum, $rate, $emissionyear, $scc\n")
        end

        return nothing
    end

    # Run monte carlo trials
    res = run(mcs, mm, trials;
        ntimesteps = ntimesteps,
        trials_output_filename = trials_output_filename,
        results_output_dir = "$output_dir/results",
        scenario_args = scenario_args,
        scenario_func = _scenario_func,
        post_trial_func = scc_calculation
        )

    return res
end

ext: jl


text:
facts("constructors") do
  ticker_str = "TICKER"
#   context("inner constructor from TimeArray") do
#     fts = FinancialSeries.FinancialTimeSeries(
#       ohlc, FinancialSeries.Stock(FinancialSeries.Ticker(ticker_str)))
#     @fact values(fts.series["Open"])[1] => roughly(104.88)
#   end
#   context("from TimeArray components and instrument object") do
#     @fact_throws FinancialSeries.FinancialTimeSeries(
#       [DateTime(1795,10,31)], [0.0 0.0], ["col-1" "col-2"], "just-a-string")
#     fts = FinancialSeries.FinancialTimeSeries(
#       [DateTime(1795,10,31)], [0.0 0.0], ["col-1", "col-2"],
#       FinancialSeries.Stock(FinancialSeries.Ticker(ticker_str)))
#     @fact isa(fts.instrument, FinancialSeries.Stock) => true
#     @fact string(fts.instrument.ticker) => ticker_str
#   end

  ### TODO: redo these tests (or delete?) now that TimeArray has meta
  @fact 0 => 1
end

facts("element and time series method wrappers") do
  inst = FinancialSeries.Stock(FinancialSeries.Ticker("AAPL"))
  fts_cl = FinancialSeries.FinancialTimeSeries(cl, inst)
  fts_ohlc = FinancialSeries.FinancialTimeSeries(ohlc, inst)
#   context("element wrappers isolate elements") do
#     # correct types
#     @fact isa(timestamp(fts_cl), Array{Date,1}) => true
#     @fact isa(values(fts_cl), Array{Float64,1}) => true
#     @fact isa(values(fts_ohlc), Array{Float64,2}) => true
#     @fact isa(colnames(cl), Array{UTF8String, 1}) => true
#     @fact isa(instrument(fts_cl), FinancialSeries.Stock) => true
#     # proper indexing
#     @fact values(fts_cl)[2] => roughly(102.5)
#     @fact timestamp(fts_cl)[3] => DateTime.Date(2000,1,5)
#   end

  ### TODO: fix wrapper functions
  @fact 0 => 1
end

facts("time series method wrappers") do
#   inst = FinancialSeries.Stock(FinancialSeries.Ticker("AAPL"))
#   fts_cl = FinancialSeries.FinancialTimeSeries(cl, inst)
#   context("correct simple return value") do
#     @fact percentchange(fts_cl).values[1] => roughly((102.5-111.94)/111.94)
#   end
#   context("correct log return value") do
#     @fact percentchange(fts_cl, method="log").values[1] => roughly(log(102.5) - log(111.94))
#   end

  ### TODO: fix wrapper functions

  ## TODO: add tests for preserving meta after percentchange etc. calls
  @fact 0 => 1
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Random: randstring

@testset "ifelse" begin
    @test ifelse(true, 1, 2) == 1
    @test ifelse(false, 1, 2) == 2

    let s = Set()
        ifelse(true, push!(s, 1), push!(s, 2))
        @test s == Set([1, 2])
    end

    let s = Set()
        true ? push!(s, 1) : push!(s, 2)
        false ? push!(s, 3) : push!(s, 4)
        @test s == Set([1, 4])
    end

    let B = [true true false]
        @test ifelse.(B, 1, 2) == [1 1 2]
        @test ifelse.(B, 1, [2 3 4]) == [1 1 4]
        @test ifelse.(B, [2 3 4], 1) == [2 3 1]
        @test ifelse.(B, [2 3 4], [5 6 7]) == [2 3 7]
    end
end

@testset "operations on Pairs" begin
    @test reverse(Pair(1,2)) == Pair(2,1)
    @test reverse(Pair("13","24")) == Pair("24","13")
    @test typeof(reverse(Pair{String,Int64}("a",1))) == Pair{Int64,String}
    @test convert(Pair{Float64,Float64}, 17 => 4711) === (17.0 => 4711.0)
    @test convert(Pair{Int,Float64}, 17 => 4711) === (17 => 4711.0)
    @test convert(Pair{Float64,Int}, 17 => 4711) === (17.0 => 4711)
    @test convert(Pair{Any,Any}, 17 => 4711) === Pair{Any,Any}(17, 4711)
    @test convert(Pair{Number,Number}, 17 => 4711) === Pair{Number,Number}(17, 4711)
    @test promote(1=>1, 2=>2.0) === (1=>1.0, 2=>2.0)
    @test promote(1=>1, 2.0=>2) === (1.0=>1, 2.0=>2)
    @test promote(1=>1.0, 2.0=>2) === (1.0=>1.0, 2.0=>2.0)
    @test promote(1=>1, :b=>2.0) === (Pair{Any,Float64}(1,1.0),Pair{Any,Float64}(:b,2.0))
    @test isa([:a=>1, :b=>2], Vector{Pair{Symbol,Int}})
    @test isa([:a=>1, :b=>2.0], Vector{Pair{Symbol,Float64}})
    @test isa(["a"=>1, :b=>2.0], Vector{Pair{Any,Float64}})

    p = 1=>:foo
    @test first(p) == 1
    @test last(p)  == :foo
    @test first(reverse(p)) == :foo
    @test last(reverse(p))  == 1
    @test lastindex(p) == 2
    @test p[lastindex(p)] == p[end] == p[2] == :foo
end

# Infix `isa`
@test 1 isa Integer

@test (|)(2) == 2
@test xor(2) == 2
@test (⊻)(2) == 2

@test_throws MethodError min(Set([1]), Set([2]))
@test_throws MethodError max(Set([1]), Set([2]))
@test_throws MethodError minmax(Set([1]), Set([2]))

# Test if isless (not <) is used by min, max, minmax
# and commutativity.
struct TO23094
    x::Int
end
Base.isless(a::TO23094, b::TO23094) = isless(a.x, b.x)
Base.isequal(a::TO23094, b::TO23094) = isequal(a.x, b.x)
import Base.<
<(a::TO23094, b::TO23094) = error("< should not be called")

@test isequal(min(TO23094(1), TO23094(2)), TO23094(1))
@test isequal(min(TO23094(2), TO23094(1)), TO23094(1))
@test isequal(max(TO23094(1), TO23094(2)), TO23094(2))
@test isequal(max(TO23094(2), TO23094(1)), TO23094(2))
@test isequal(minmax(TO23094(1), TO23094(2))[1], TO23094(1))
@test isequal(minmax(TO23094(1), TO23094(2))[2], TO23094(2))
@test isequal(minmax(TO23094(2), TO23094(1))[1], TO23094(1))
@test isequal(minmax(TO23094(2), TO23094(1))[2], TO23094(2))

@test isless('a','b')

@testset "isgreater" begin
    # isgreater should be compatible with min.
    min1(a, b) = Base.isgreater(a, b) ? b : a
    # min promotes numerical arguments to the same type, but our quick min1
    # doesn't, so use float test values instead of ints.
    values = (1.0, 5.0, NaN, missing, Inf)
    for a in values, b in values
        @test min(a, b) === min1(a, b)
        @test min((a,), (b,)) === min1((a,), (b,))
        @test all(min([a], [b]) .=== min1([a], [b]))
    end
end

@testset "isunordered" begin
    @test  isunordered(NaN)
    @test  isunordered(NaN32)
    @test  isunordered(missing)
    @test !isunordered(1)
    @test !isunordered([NaN, 1])
    @test !isunordered([1.0, missing])
end

@testset "vectorized comparisons between numbers" begin
    @test 1 .!= 2
    @test 1 .== 1
    @test 1 .< 2
    @test 1 .<= 2
end

# issue #13144: max() with 4 or more array arguments
let xs = [[i:i+4;] for i in 1:10]
    for n in 2:10
        @test max.(xs[1:n]...) == [n:n+4;]
    end
end

# issue #19714
struct T19714 <: Integer end
Base.float(::T19714) = 19714.0
Base.:/(::T19714, ::T19714) = T19714()
Base.convert(::Type{T19714}, ::Int) = T19714()
Base.promote_rule(::Type{T19714}, ::Type{Int}) = T19714
@test T19714()/1 === 1/T19714() === T19714()

# pr #17155 and #33568
@testset "function composition" begin
    @test (uppercase∘(x->string(x,base=16)))(239487) == "3A77F"
    @test ∘(x -> x-2, x -> x-3, x -> x+5)(7) == 7
    fs = [x -> x[1:2], uppercase, lowercase]
    @test ∘(fs...)("ABC") == "AB"

    # Like +() and *() we leave ∘() undefined.
    # While `∘() = identity` is a reasonable definition for functions, this
    # would cause headaches for composition of user defined morphisms.
    # See also #34251
    @test_throws(MethodError, ∘())

    @test ∘(x -> (x, 1))(0) === (0, 1)
    @test ∘(x -> (x, 2), x -> (x, 1))(0) === ((0, 1), 2)
    @test ∘(x -> (x, 3), x -> (x, 2), x->(x,1))(0) === (((0, 1), 2), 3)
    @test ∘(x -> (x, 4), x -> (x, 3), x->(x,2), x-> (x, 1))(0) === ((((0, 1), 2), 3), 4)

    # test that user defined functors only need to overload the two arg version
    struct FreeMagma
        word
    end
    Base.:(∘)(a::FreeMagma, b::FreeMagma) = FreeMagma((a.word, b.word))

    @test ∘(FreeMagma(1)) === FreeMagma(1)
    @test ∘(FreeMagma(1), FreeMagma(2)) === FreeMagma((1,2))
    @test ∘(FreeMagma(1), FreeMagma(2), FreeMagma(3)) === FreeMagma(((1,2), 3))
    @test ∘(FreeMagma(1), FreeMagma(2), FreeMagma(3), FreeMagma(4)) === FreeMagma((((1,2), 3), 4))

    @test fieldtypes(typeof(Float64 ∘ Int)) == (Type{Float64}, Type{Int})

    @test repr(uppercase ∘ first) == "uppercase ∘ first"
    @test sprint(show, "text/plain", uppercase ∘ first) == "uppercase ∘ first"

    # test keyword ags in composition
    function kwf(a;b,c); a + b + c; end
    @test (abs2 ∘ kwf)(1,b=2,c=3) == 36

end

@testset "function negation" begin
    str = randstring(20)
    @test filter(!isuppercase, str) == replace(str, r"[A-Z]" => "")
    @test filter(!islowercase, str) == replace(str, r"[a-z]" => "")
end

# issue #19891
@testset "chained comparison" begin
    B = 0 .< [1 -1 5] .< 3
    @test B == [true false false]
    B = 3 .> [1 -1 5] .> 0
    @test B == [true false false]
end

struct TypeWrapper
    t::Type
end
Base.:(<)(x::TypeWrapper, y::TypeWrapper) = (x.t <: y.t) & (x.t != y.t)
@testset "poset" begin
    #   Real
    #  /    \
    # Int  Float64
    #  \    /
    #  Union{}
    @test TypeWrapper(Int) <= TypeWrapper(Int)
    @test TypeWrapper(Int) <= TypeWrapper(Real)
    @test !(TypeWrapper(Int) <= TypeWrapper(Float64))
end

# issue #20355
@testset "mod1, fld1" begin
    for T in [Int8, Int16, Int32, Int64],
        x in T[typemin(T); typemin(T) + 1; -10:10; typemax(T)-1; typemax(T)],
        y in T[typemin(T); typemin(T) + 1; -10:-1; 1:10; typemax(T)-1; typemax(T)]

        m = mod1(x, y)
        @test mod(x, y) == mod(m, y)
        if y > 0
            @test 0 < m <= y
        else
            @test y <= m < 0
        end
        if x == typemin(T) && y == -1
            @test_throws DivideError fld1(x, y)
        else
            f = fld1(x, y)
            @test (f - 1) * y + m == x
        end
    end

    for T in [UInt8, UInt16, UInt32, UInt64],
        x in T[0:10; typemax(T)-1; typemax(T)],
        y in T[1:10; typemax(T)-1; typemax(T)]

        m = mod1(x, y)
        @test mod(x, y) == mod(m, y)
        @test 0 < m <= y
        f = fld1(x, y)
        @test (f - 1) * y + m == x
    end

    for T in [Float32, Float64, Rational{Int64}],
        x in T[k // 4 for k in -10:10],
        y in T[k // 4 for k in [-10:-1; 1:10]]

        m = mod1(x, y)
        @test mod(x, y) == mod(m, y)
        if y > 0
            @test 0 < m <= y
        else
            @test y <= m < 0
        end
        f = fld1(x, y)
        @test (f - 1) * y + m == x
    end

    @test fldmod1(4.0, 3) == fldmod1(4, 3)
end

@testset "Fix12" begin
    x = 9
    y = 7.0
    fx = Base.Fix1(/, x)
    fy = Base.Fix2(/, y)
    @test fx(y) == x / y
    @test fy(x) == x / y
end

@testset "curried comparisons" begin
    eql5 = (==)(5)
    neq5 = (!=)(5)
    gte5 = (>=)(5)
    lte5 = (<=)(5)
    gt5  = (>)(5)
    lt5  = (<)(5)

    @test eql5(5) && !eql5(0)
    @test neq5(6) && !neq5(5)
    @test gte5(5) && gte5(6)
    @test lte5(5) && lte5(4)
    @test gt5(6) && !gt5(5)
    @test lt5(4) && !lt5(5)
end

@testset "ni" begin
    @test ∋([1,5,10,11], 5)
    @test !∋([1,10,11], 5)
    @test ∋(5)([5,1])
    @test !∋(42)([0,1,100])
    @test ∌(0)(1:10)
    @test ∋(0)(-2:2)
end

a = rand(3, 3)
@test transpose(a) === a'ᵀ

@test [Base.afoldl(+, 1:i...) for i = 1:40] == [i * (i + 1) ÷ 2 for i = 1:40]


text:
# bearing.jl
# adding noise to bearing modality

struct NoisyBearing <: FEBOL.Sensor
    bo::BearingOnly
end

NoisyBearing(ns::Real) = NoisyBearing( BearingOnly(ns) )

function FEBOL.observe(theta::FEBOL.LocTuple, s::NoisyBearing, p::Pose)
    truth = true_bearing(p, theta)
    noise = s.bo.noise_sigma * randn()
    if p[1] < 200.0

        # determine line to radio source
        dx = theta[1] - p[1]        # x-distance to jammer
        dy = theta[2] - p[2]        # y-distance to jammer

        y_at_200 = p[2] + (dy / dx) * (200.0 - p[1])

        # we get the wall...
        if y_at_200 > 150.0

            y1 = p[2]
            y2 = theta[2]
            d = theta[1] - p[1]
            x2 = d / (1.0 + y2/y1)

            r_x = theta[1] - x2         # location of reflection
            truth = true_bearing(p, (r_x,0.0))
            noise = 15.0 * randn()
        end
    end
    return mod(truth + noise, 360.0)
end

# This shouldn't matter, because we only observe
#function O(bo::NoisyBearing, theta::TargetTuple, p::Pose, o)
#
#    # Calculate true bearing, and find distance to bin edges
#    ang_deg = true_bearing(p, theta)
#    o_diff = fit_180(o - ang_deg)
#
#    # now look at probability
#    d = Normal(0, bo.noise_sigma)
#    p = pdf(d, o_diff)
#    return p
#end


text:
using AccuracyAtTop: hinge, quadratic

# Surrogate functions
ϑ = rand()

@testset "surrrogate functions" begin
    @test hinge.([-1/ϑ - 1, 0, -1/ϑ + 1], ϑ) ≈ [0, 1, ϑ]
    @test quadratic.([-1/ϑ - 1, 0, -1/ϑ + 1], ϑ) ≈ [0, 1, ϑ^2]
end


text:
"""
after all processing we need to establish weather the mask is full  (in case of first iteration also is it empty)
"""
module IsBlockToStayActive
using CUDA, ..CUDAGpuUtils, Logging,StaticArrays,..HFUtils
export isActiveForFirstPass, isActiveForNormalPass
"""
We check is block still active in first pass - we need to check is empty or full and set the metadata
    for reference below te fourth dimension of metadata
        1) isActiveOrFullForPasssegm - true if other  mask is acive for modifications
        2) isActiveOrFullForPassgold -  true if gold standard mask is acive for modifications
        3) isFullSegm - true if other mask is full (only ones)
        4) isFullGold - true if gold standard mask is full (only ones)

        shmem- shared memory
        metaData - 3 dim array with metadata of data blocks
        currMatadataBlockX,currMatadataBlockY, currMatadataBlockZ - cartesian coordinates of current block in metadaa!!!
        mainQuesCounter - counter that we will update atomically and will be usefull to populate the work queue
        mainWorkQueue - the list of the indicies of  data blocks in metadata with additional information is it referencing the goldpass or second one

"""
function isActiveForFirstPass(isMaskFull::MVector{1,Bool}
                            , isMaskEmpty::MVector{1,Bool}
                            ,shmem
                            ,currMatadataBlockX::UInt8
                            ,currMatadataBlockY::UInt8
                            ,currMatadataBlockZ::UInt8
                            ,isPassGold::Bool
                            ,metaData
                            ,mainQuesCounter
                            ,mainWorkQueue )
    @inbounds shmem[1,threadIdxY(),20]=  reduce_warp_and(isMaskFull, UInt8(32))
    @inbounds shmem[1,threadIdxY(),21]=  reduce_warp_and(isMaskEmpty, UInt8(32))
    #so now we have 32 booleans in shared memory so we need to reduce it one more time using single warp
    #TODO() check weather warps are column or row wise this below also I am not sure is it good

    if(threadIdxY()==1 && threadIdxX()==1)
        @inbounds  shmem[1,1,22]= reduce_warp_and(shmem[1,threadIdxY(),20], UInt8(32))  #true if mask is full
    end
    if(threadIdxY()==2 && threadIdxX()==2)
        @inbounds  shmem[1,2,22]= reduce_warp_and(shmem[1,threadIdxY(),21], UInt8(32))  #true if mask is empty
    end
    sync_threads()
    #we have needed data in shmem

    # we update matedata concurrently
    if(threadIdxY()==1 && threadIdxX()==1 && (shmem[1,2,22] || shmem[1,1,22]))
        metaData[currMatadataBlockX,currMatadataBlockY,currMatadataBlockZ,isPassGold+1]=false # we set is inactive
    end#if
    if(threadIdxY()==1 && threadIdxX()==1 && (shmem[1,2,22] || shmem[1,1,22]))
        metaData[currMatadataBlockX,currMatadataBlockY,currMatadataBlockZ,isPassGold+3]=true # we set is as full
    end#if
    #so in case it not empty and not full we need to put it into the work queue and increment appropriate counter
    if(threadIdxY()==3&& threadIdxX()==3 && !(shmem[1,2,22] || shmem[1,1,22]))
        mainWorkQueue[CUDA.atomic_inc!(pointer(mainQuesCounter), UInt16(1))+1]=[currMatadataBlockX,currMatadataBlockY,currMatadataBlockZ,UInt8(isPassGold)] #x dim of block in metadata
    end#if

end#isActiveForFirstPass

"""
We check is block still active in first pass - we need to check is full (if it would be empty we would not start to analyze it after first pass)
shmem- shared memory
metaData - 3 dim array with metadata of data blocks
currMatadataBlockX,currMatadataBlockY, currMatadataBlockZ - cartesian coordinates of current block in metadaa!!!
mainQuesCounter - counter that we will update atomically and will be usefull to populate the work queue
mainWorkQueue - the list of the indicies of  data blocks in metadata with additional information is it referencing the goldpass or second one
workingQueuePosition - points out to the index of currently analyzed spot in work queue that we are in order to

"""
function isActiveForNormalPass(isMaskFull::MVector{1,Bool}
                                ,shmem
                                ,currMatadataBlockX::UInt8
                                ,currMatadataBlockY::UInt8
                                ,currMatadataBlockZ::UInt8
                                ,isPassGold::Bool
                                ,metaData
                                ,mainQuesCounter
                                ,mainWorkQueue )
    @inbounds shmem[1,threadIdxY(),20]=  reduce_warp_and(isMaskFull, UInt8(32))
    #so now we have 32 booleans in shared memory so we need to reduce it one more time using single warp
    #TODO() check weather warps are column or row wise this below also I am not sure is it good
    if(threadIdxY()==1 && threadIdxX()==1)
        @inbounds  shmem[1,1,22]= reduce_warp_and(shmem[1,threadIdxY(),20], UInt8(32))  #true if mask is full
    end
    sync_threads()
    #we have needed data in shmem
    if(threadIdxY()==1 && threadIdxX()==1 && (shmem[1,2,22]))
        metaData[currMatadataBlockX,currMatadataBlockY,currMatadataBlockZ,isPassGold+1]=false # we set is inactive
    end#if
    if(threadIdxY()==2 && threadIdxX()==1 && (shmem[1,2,22]))
        metaData[currMatadataBlockX,currMatadataBlockY,currMatadataBlockZ,isPassGold+3]=true # we set is as full
    end#if

    #so in case it not empty and not full we need to put it into the work queue and increment appropriate counter
    if(threadIdxY()==3&& threadIdxX()==3 && !(shmem[1,2,22]))
        mainWorkQueue[CUDA.atomic_inc!(pointer(mainQuesCounter), UInt16(1))+1]=[currMatadataBlockX,currMatadataBlockY,currMatadataBlockZ,UInt8(isPassGold)] #x dim of block in metadata
    end#if

end#isActiveForNormalPass

end#IsBlockToStayActive


text:
# Condition and ALU flags for conditional execution, etc.
#
# Make this a struct with an API to possibly reimplement as a single integer
# with bit-twiddling instead of a struct of Bools. Probably not needed but
# may save some space in the future.
struct ALUFlags
    carry       :: Bool
    negative    :: Bool
    overflow    :: Bool
    zero        :: Bool
end
ALUFlags() = ALUFlags(false, false, false, false)

# ------------------------- #
# --- Address Generator --- #
# ------------------------- #

# Things to think about:
#
# * Maybe make the "number of bits" for the address generator a parameter if
# we ever want have generators outside of the 8-bit range of Asap4
#
# * Think about write-back timings.

@with_kw_noshow mutable struct AddressGenerator
    start       ::UInt8 = 0
    stop        ::UInt8 = 0
    current     ::UInt8 = 0
    # Old value of the pointer. For rolling back in case of branch mispredict.
    old         ::UInt8 = 0
    stride      ::UInt8 = 0
    # Flag to indicate this needs an increment. Used because multiple operands
    # in an instruction can read from the same address generator.
    needs_increment :: Bool = false
end

function set!(ag::AddressGenerator, start, stop)
    ag.start = start
    ag.stop  = stop
    ag.current = ag.start
    return nothing
end

function set_stride!(ag::AddressGenerator, stride)
    ag.stride = stride
    return nothing
end

Base.read(ag::AddressGenerator) = ag.current
mark(ag::AddressGenerator) = ag.needs_increment = true

function increment!(ag::AddressGenerator)
    # If the address generator does not need to be incremented, do nothing and
    # just return.
    ag.needs_increment || return nothing
    # Clear increment flag.
    ag.needs_increment = false

    # Save the old pointer value.
    ag.old = ag.current

    # Take advantage of the binary arithmetic of julia - pointer will
    # automatically wrap around at 256.
    ag.current = (ag.current == ag.stop) ? ag.start : ag.current + ag.stride
    return nothing
end
rollback!(ag::AddressGenerator) = ag.current = ag.old

mutable struct CondExec
    flag        :: Bool
    mask        :: Int16
    # One of :OR, :AND, :XOR
    unary_op    :: Symbol
    early_kill  :: Bool
end
CondExec() = CondExec(false, 0, :OR, false)

function set!(c::CondExec, mask)
    # Get bit 13 for the early kill.
    c.early_kill = isbitset(mask, 13)

    # Get bits 12 and 11 for the reduction operator.
    #
    # If OO or 01, use OR
    # If 10, use AND
    # If 11, use XOR
    if isbitset(mask, 12)
        if isbitset(mask, 11)
            c.unary_op = :XOR
        else
            c.unary_op = :AND
        end
    else
        c.unary_op = :OR
    end
    c.mask = mask
    return nothing
end

# For keeping track of instructions through the pipeline.
# Essentially just a wrapper for the instruction with an extra slot for the
# result.
struct PipelineEntry
    instruction :: AsapInstruction# = NOP()
    # The actual source and result values
    src1_value  :: Int16# = 0
    src2_value  :: Int16# = 0
    result      :: Int16# = 0

    # Data for rolling-back branches. This doesn't really need to be recorded
    # at every stage, but it will be cleaner to do it this way. It doesn't take
    # that much memory anyways.
    old_return_address :: Int64# = 0
    old_alternate_pc   :: Int64# = 0
    old_repeat_count   :: Int64# = 0
end

PipelineEntry(instruction = InstNOP()) = PipelineEntry(instruction,0,0,0,0,0,0)

# Convenience types for changing fields of the PipelineEntry immutable
# struct.
struct SRC <: Mutator{PipelineEntry}
    src1_value :: Int16
    src2_value :: Int16
end
struct RESULT <: Mutator{PipelineEntry}
    result::Int16
end
struct INSTRUCTION <: Mutator{PipelineEntry}
    instruction::AsapInstruction
end
struct RETURNS <: Mutator{PipelineEntry}
    old_return_address :: Int64
    old_alternate_pc   :: Int64
    old_repeat_count   :: Int64
end

# Asap Pipeline.
#
# Note: The data in each stage is the END of each stage.
mutable struct AsapPipeline
    stage1 :: PipelineEntry
    stage2 :: PipelineEntry
    stage3 :: PipelineEntry
    stage4 :: PipelineEntry
    stage5 :: PipelineEntry
    stage6 :: PipelineEntry
    stage7 :: PipelineEntry
end

# Initialize to an empty pipeline.
AsapPipeline() = AsapPipeline(
    PipelineEntry(),
    PipelineEntry(),
    PipelineEntry(),
    PipelineEntry(),
    PipelineEntry(),
    PipelineEntry(),
    PipelineEntry(),
)

mutable struct SmallAssociative{K,V}
    items::Vector{Pair{K,V}}
end
SmallAssociative{K,V}() where {K,V} = SmallAssociative(Vector{Pair{K,V}}())
Base.length(s::SmallAssociative) = length(s.items)

function Base.getindex(s::SmallAssociative{K,V}, key :: K) where {K,V}
    # Iterate through the wrapped container - try to find the given key.
    for (k, v) in s.items
        if k == key
            return v
        end
    end
    throw(KeyError(key))
end

function Base.setindex!(s::SmallAssociative, value, key)
    # Check to see if this key already exists
    foundkey = false
    for (i, (k, v)) in enumerate(s.items)
        if k == key
            s.items[i] = value
            foundkey = true
            break
        end
    end

    # If key was not found, make a new entry
    if !foundkey
        push!(s.items, key => value)
    end
    return nothing
end

# Iterator protocol
function Base.iterate(s::SmallAssociative, state = 1)
    state > length(s.items) ? nothing : (s.items[state], state + 1)
end

# Key and Value iterators
struct SmallKeyIterator{K,V}
    s :: SmallAssociative{K,V}
end

Base.keys(s::SmallAssociative) = SmallKeyIterator(s)

function Base.iterate(i::SmallKeyIterator, state = 1)
    state > length(i.s.items) ? nothing : (first(i.s.items[state]), state + 1)
end

struct SmallValueIterator{K,V}
    s :: SmallAssociative{K,V}
end

Base.values(s::SmallAssociative) = SmallValueIterator(s)
function Base.iterate(i::SmallValueIterator, state = 1)
    state > length(i.s.items) ? nothing : (last(i.s.items[state]), state + 1)
end

# Basic AsapCore. Parameterize based on
#
# 1. The type of FIFO - allows for injection of TestFifos to help with tests
@with_kw mutable struct AsapCore{F <: AbstractFifo}
    # --- Timing Control --- #

    # Clock period - for registering updates.
    # This number is inherently meaningless, but should usually be interpreted
    # as PS.
    #
    # Note that if interpreting this as PS ... a default value of 1 does not
    # really make a whole lot of sense.
    clockperiod ::Int64 = 1

    # --- Program --- #

    # Stored program and program counter.
    program :: AsapProgram = AsapProgram()

    # Note on Program Counter (pc) - In actual hardware, it is base 0 ...
    # because hardware is base 0. However, Julia is base 1, so we'll have to
    # keep that in mind.
    pc :: Int64 = 0

    # Hardware registers holding information about repeat values.
    repeat_count  :: Int64 = 0
    repeat_start  :: Int64 = 0
    repeat_end    :: Int64 = 0

    # Hardware return address buffer.
    return_address :: Int64 = 0

    # Number of pending NOPs to inser in stage 3 of the pipeline.
    # This is set by the options in the assembler.
    pending_nops :: Int16 = 0

    # Mispredicted branch - done in S4
    branch_mispredict :: Bool = false

    # ---------------------- #
    # Input and Output Fifos #
    # ---------------------- #

    # Input fifo - default the element types of the fifo to Int16s.
    fifos::Vector{F} = [
        DualClockFifo(Int16, 32),
        DualClockFifo(Int16, 32),
    ]

    # Set this up as a Dict to handle cases where directions do not have
    # connected outputs.
    #
    # NOTE: In Julia 0.7+, small unions are faster than in 0.6+, so it may be
    # possible to convert this to a Union{Void,DualClockFifo{Int16}} without
    # incurring any nasty runtime penalties.
    #outputs::Dict{Int,F} = Dict{Int,DualClockFifo{Int16}}()
    outputs :: SmallAssociative{Int,F} = SmallAssociative{Int,DualClockFifo{Int16}}()

    # Flags
    aluflags :: ALUFlags = ALUFlags()

    # 40-bit accumulator
    # TODO: Think about wrapping this in a type that auto-truncates to
    # a set number of bits.
    accumulator :: Int64 = 0

    # ---------------------------- #
    # Dynamic Configuration Memory #
    # ---------------------------- #

    # Address generators.
    address_generators::Vector{AddressGenerator} = [
        AddressGenerator(),
        AddressGenerator(),
        AddressGenerator(),
    ]

    # Hardware pointers
    pointers :: Vector{Int16} = zeros(Int16, 4)

    # Conditional execution blocks
    condexec :: Vector{CondExec} = [
        CondExec(),
        CondExec(),
    ]

    # Mask for output
    obuf_mask :: BitVector = falses(8)

    # ----------- #
    # Data memory #
    # ----------- #

    dmem::Vector{Int16} = zeros(Int16, 256)

    # -------- #
    # Pipeline #
    # -------- #
    pipeline::AsapPipeline = AsapPipeline()
end

clockperiod(core::AsapCore) = core.clockperiod

function setobufmask!(core, mask)
    for index in 1:length(core.obuf_mask)
        core.obuf_mask[index] = isbitset(mask, index - 1)
    end
end

# Convert from index 0 to index 1.
dmem(core, address) = core.dmem[address + 1]

################################################################################
# Stall Detection
################################################################################

"""
    StallSignals

Collection of stall signals that can occur in the Asap4 processor.

Parameters:

* `stall_01` - Stall stages 0 and 1. Happens if stalled on a fifo or if stage
    2 has pending NOPS to insert.
* `nop_2` - Insert NOPs into Stage 2. Happens when there are pending NOPs.
* `stall_234` - Stall stages 2, 3, and 4. Happens if stalled on fifo.
* `stall_5678` - Stall stages 5, 6, and 7. Happens if stalled on a fifo and none
    of the destinations in stages 5, 6, and 7 is an output.
* `nop_5` - Insert NOPS into stage 5. Happens if stalled on an output but
    one of the destinations in stages 5, 6, or 7 is an output.
"""
struct StallSignals
    stall_01   :: Bool
    nop_2      :: Bool
    stall_234  :: Bool
    stall_5678 :: Bool
    nop_5      :: Bool
end

# Use an ENUM to encode the stall reason.
@enum StallReason NoStall StallOp EmptyIbuf FullObuf FullObufMask

function stall_check(core::AsapCore)
    # Check if stalled on a FIFO.
    stallreason = stall_fifo_check(core)

    fifo_stall = stallreason != NoStall

    # Default stall signals to "false"
    stall_01 = false
    nop_2 = false
    stall_234 = false
    stall_5678 = false
    nop_5 = false

    # If we're stalled on a fifo, stall stages 0, 1, 2, 3, and 4.
    #
    # Check stages 5, 6, and 7 to see if they are writing to an output.
    # If not, stall stages 5, 6, and 7 also.
    if fifo_stall
        stall_01  = true
        stall_234 = true

        if core.pipeline.stage5.instruction.dest_is_output ||
            core.pipeline.stage6.instruction.dest_is_output ||
            core.pipeline.stage7.instruction.dest_is_output

            nop_5 = true
        else
            stall_5678 = true
        end

    # Otherwise, stall stages 0 and 1 if stage 2 has pending NOPs or will
    # request NOPs this cycle.
    else
        # Like in the RTL for asap4, if pending nops = 1, we're on the last
        # NOP, so we can stop the NOP signal then.
        if core.pending_nops > 0
            stall_01 = true
            nop_2 = true
        elseif core.pipeline.stage2.instruction.nops > 0
            stall_01 = true
        end
    end

    # Return a collection of stall signals.
    return StallSignals(
        stall_01,
        nop_2,
        stall_234,
        stall_5678,
        nop_5,
    )
end

"""
    stall_fifo_check(core::AsapCore)

Return `true` if the core should stall, `false` otherwise. A core will stall if:

* The instruction in stage 4 is a STALL instruction. In this case, the core will
    only stall if all masked inputs are empty and all masked outputs are full.

* The instruction in stage 4 is NOT a STALL instruction. Then accesses to inputs
    and outputs will be checked. A stall will be generated if any read input is
    empty or any written-to output is full.
"""
function stall_fifo_check(core::AsapCore) :: StallReason
    # Early termination of stage 3 is a NOP - NOPs can't stall.
    core.pipeline.stage3.instruction.op == NOP && (return NoStall)

    # Get the pipeline stage 4 instruction.
    stage3 = core.pipeline.stage3

    # If this instruction is a STALL instruction, pass its mask is found in
    # the src1_value field
    if stage3.instruction.op == STALL
        return stall_check_stall_op(core, stage3.src1_value)

    # Otherwise, do a normal stall check.
    else
        status = stall_check_io(core, stage3.instruction)
        return status
    end
end

"""
    stall_check_stall_op(core::AsapCore, mask::Integer)

Return `true` if `core` should stall as a result of a STALL instruction.
Argument `mask` is the stall bit mask indicating which hardware resources
should be checked.

Core will stall only if all hardware resources checked are stalling. That is,
all inputs are empty and all outputs are full.

The bit encoding is:

* Bit 0: Stall on empty input fifo 0
* Bit 1: Stall on empty input fifo 1
* Bit 2: Stall on empty packet in (TODO)
* Bit 3: Stall on empty memory in (TODO)
* Bit 4: Stall on all obuf_mask direction full
* Bit 5: Stall on full packet out (TODO)
* Bit 6: Stall on full memory out (TODO)
"""
function stall_check_stall_op(core::AsapCore, mask)
    # Set default stall to be true. If any further evaluations break the default
    # We can terminate early.
    default = true

    # Start checking all of the bits. This is a little painful, but I can't
    # really think of a more elegant way to do this.

    # --- Check inputs ---
    if isbitset(mask, 0)
        # Note: Converting from index 0 to index 1
        stall_check_ibuf(core, 1, default) != default && return NoStall
    end

    if isbitset(mask, 1)
        # Note: Converting from index 0 tio index 1
        stall_check_ibuf(core, 2, default) != default && return NoStall
    end

    # --- Check obuf mask ---
    if isbitset(mask, 4)
        # Iterate through the OBUF mask. If a flag is set, check that output.
        for (index, flag) in enumerate(core.obuf_mask)
            flag || continue
            stall_check_obuf(core, index, default) != default && return NoStall
        end
    end

    # Give "unimplemented error" for unimplemented flags.
    if mask & (0b1101100) != 0
        error("""
        Bits 2, 3, 5, and 6 of the STALL instruction mask are not yet
        implemented.
        """)
    end

    return StallOp
end

function stall_check_io(core::AsapCore, instruction::AsapInstruction)
    # Default to not stall. Exit early if any IO interaction should cause a
    # stall to happen.
    default = false

    # --- Check inputs ---

    # NOTE: Must do conversion from index 0 to index 1
    if sym(instruction.src1) == IBUF || sym(instruction.src1) == IBUF_NEXT
        # The fifo to check should be in the src1_index field of the instruction.
        stall_check_ibuf(core.fifos[ind(instruction.src1) + 1], default) != default && return EmptyIbuf
    end

    if sym(instruction.src2) == IBUF || sym(instruction.src2) == IBUF_NEXT
        stall_check_ibuf(core.fifos[ind(instruction.src2) + 1], default) != default && return EmptyIbuf
    end

    # Check writes to an output
    if sym(instruction.dest) == OUTPUT
        stall_check_obuf(core.outputs[ind(instruction.dest) + 1], default) != default && return FullObuf
    # Check if doing a broadcast. Then, stall if any of the outputs selected
    # by the obuf_mask are stalling.
    elseif sym(instruction.dest) == OBUF
        for (index, flag) in enumerate(core.obuf_mask)
            flag || continue
            stall_check_obuf(core.outputs[index], default) != default && return FullObufMask
        end
    end

    # All checks passed, don't stall!
    return NoStall
end

# Break out the various stall checks into a bunch of little functions.

@inline function stall_check_ibuf(fifo, default) :: Bool
    # If the default value is to stall, return false if the given fifo is
    # not empty to go against the stall.
    if default == true && isreadready(fifo)
        return false

    # Otherwise, if the default is not to stall, return "true" if the fifo is
    # empty to indicate the core should stall
    elseif default == false && !isreadready(fifo)
        return true

    # If the above conditions are not met, just return the default.
    else
        return default
    end
end

@inline function stall_check_obuf(fifo, default :: Bool) :: Bool
    # Logic is similar to the ibuf check, but looks for full fifos instead of
    # empty fifos.
    if default == true && iswriteready(fifo)
        return false
    elseif default == false && !iswriteready(fifo)
        return true
    else
        return default
    end
end

"""
    checkoutputs(core::AsapCore)

Return `true` if any of the output buffers selected by `core`'s obuf_mask is
full. Otherwise, return `false`.
"""
function checkoutputs(core::AsapCore)
    for (index, bit) in enumerate(core.obuf_mask)
        bit || continue
        if !iswriteready(core.outputs[index])
            return true
        end
    end
    return false
end


text:
#=
    Decomposition of global constraint matrix_element in Julia ConstraintSolver.jl

    matrix_element(model,x,i,j,val)

    ensures that  x[i,j] = val

    Model created by Hakan Kjellerstrand, hakank@gmail.com
    See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# matrix_element is in constraints_utils.jl
#
function matrix_element_test(n=4,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            "branch_split"=>:Biggest, # fastest for this simple model!
                                                            # "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>3,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    m = n
    @variable(model, 1 <= x[1:n,1:m] <= n, Int)

    # @constraint(model, x[:] in CS.AllDifferent())
    @variable(model, 1 <= ii <= n, Int)
    @variable(model, 1 <= jj <= m, Int)
    @variable(model, 1 <= val <= n, Int)

    for i in 1:n
        @constraint(model, x[i,:] in CS.AllDifferent())
        @constraint(model, x[:,i] in CS.AllDifferent())
    end

    @constraint(model, val == n)
    # @constraint(model, ii == 2)
    # @constraint(model, jj == 3)

    matrix_element(model,x,ii,jj,val)

    # @objective(model,Min,x[n,m])

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                iix = convert.(Integer,JuMP.value.(ii; result=sol))
                jjx = convert.(Integer,JuMP.value.(jj; result=sol))
                valx = convert.(Integer,JuMP.value.(val; result=sol))

                # print_grid(xx)
                display(xx)
                println("ii:$iix jj:$jjx val:$valx")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

print_solutions=true
all_solutions=false
for n in 2:20
    println("\nn:$n")
    @time status = matrix_element_test(n,print_solutions,all_solutions)
    if status != MOI.OPTIMAL
        break
    end
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# on vérifie simplement que le constructeur par défaut fait son travail
test_meta = StoppingMeta()

@test test_meta.optimality0         == 1.0
@test test_meta.unbounded_threshold == 1.0e50
@test test_meta.unbounded_x         == 1.0e50
@test test_meta.max_f               == 9223372036854775807
@test test_meta.max_eval            == 20_000
@test test_meta.max_iter            == 5_000
@test test_meta.max_time            == 300.0
@test isnan(test_meta.start_time)
@test test_meta.fail_sub_pb         == false
@test test_meta.unbounded           == false
@test test_meta.unbounded_pb        == false
@test test_meta.tired               == false
@test test_meta.stalled             == false
@test test_meta.iteration_limit     == false
@test test_meta.optimal             == false
@test test_meta.suboptimal          == false
@test test_meta.main_pb             == false
@test test_meta.nb_of_stop          == 0

try
 StoppingMeta(tol_check = x -> x)
 @test false #tol_check must have 3 arguments and not one
catch
 @test true
end


text:
import Libdl, LibGit2, LinearAlgebra

# Use this version of Elemental
Elsha = "79987d38b04838acf6b6195be1967177521ee908"

if Sys.iswindows()
    error("Elemental only works on Unix Platforms")
end

depdir = dirname(@__FILE__)

if !isdir(joinpath(depdir, "src"))
    mkdir(joinpath(depdir, "src"))
end
srcdir = joinpath(depdir, "src", "Elemental")

if !isdir(joinpath(depdir, "usr"))
    mkdir(joinpath(depdir, "usr"))
end
prefix = joinpath(depdir, "usr")

if !isdir(srcdir)
    LibGit2.clone("https://github.com/elemental/Elemental.git", "$srcdir")
end
cd(srcdir) do
    LibGit2.checkout!(LibGit2.GitRepo("."), "$Elsha")
end

LinearAlgebra.BLAS.check()
blas = LinearAlgebra.BLAS.vendor()
mathlib = Libdl.dlpath(LinearAlgebra.BLAS.libblas)
blas64 = LinearAlgebra.USE_BLAS64 ? "ON" : "OFF"
blas_suffix = blas === :openblas64 ? "_64_" : "_"
build_procs = (haskey(ENV, "CI") && ENV["CI"] == "true") ? 2 : Sys.CPU_THREADS

builddir = joinpath(depdir, "builds")
if isdir(builddir)
    rm(builddir, recursive=true)
end
mkdir(builddir)

cd(builddir) do
    run(`cmake -D CMAKE_INSTALL_PREFIX=$prefix
               -D INSTALL_PYTHON_PACKAGE=OFF
               -D PYTHON_EXECUTABLE=""
               -D PYTHON_SITE_PACKAGES=""
               -D EL_USE_64BIT_INTS=$blas64
               -D EL_USE_64BIT_BLAS_INTS=$blas64
               -D MATH_LIBS=$mathlib
               -D EL_BLAS_SUFFIX=$blas_suffix
               -D EL_LAPACK_SUFFIX=$blas_suffix
               -D CMAKE_INSTALL_RPATH=$prefix/lib
               $srcdir`)
    run(`make -j $build_procs`)
    run(`make install`)
end
GC.gc() # work-around for https://github.com/JuliaLang/julia/issues/28306


text:
using ReactionNetworks, Distributions
using Distances
using StatsBase
using DecisionTree
using StructArrays

k1, k2, k3 = 50.0, 50.0, 1.0
s0, e0, c0, p0 = 100, 5, 0, 0

tspan = (0.0, 1000.0)

const DATA = [
    1.73,
    3.80,
    5.95,
    8.10,
    11.17,
    12.92,
    15.50,
    17.75,
    20.17,
    23.67,
]

function f_lo(y_lo, θ)
    k1 = θ[1]
    k2 = θ[2]
    k3 = θ[3]
    _lo = EnzymeModels.EnzymeMM(k1, k2, k3, e0)
    _lo_out = @timed simulate(_lo, [s0, p0], tspan)

    t_lo, x_lo, coupling_lo = _lo_out.value
    c_lo = _lo_out.time

    for i in 1:10
        idx = findfirst(==(10*i), selectdim(x_lo, 1, 2))
        y_lo[i] = t_lo[idx]
    end
    return c_lo, coupling_lo
end

function f_hi(y_hi, θ)
    k1 = θ[1]
    k2 = θ[2]
    k3 = θ[3]
    _hi = EnzymeModels.Enzyme(k1, k2, k3)
    _hi_out = @timed simulate(_hi, [s0, e0, c0, p0], tspan)

    t_hi, x_hi, coupling_hi = _hi_out.value
    c_hi = _hi_out.time

    for i in 1:10
        idx = findfirst(==(10*i), selectdim(x_hi, 1, 4))
        y_hi[i] = t_hi[idx]
    end

    return c_hi, coupling_hi
end

function f_hi(y_hi, θ, coupling_lo)
    k1 = θ[1]
    k2 = θ[2]
    k3 = θ[3]
    _hi = EnzymeModels.Enzyme(k1, k2, k3)

    _hi_out = @timed simulate(_hi, [s0, e0, c0, p0], tspan, (nothing, nothing, coupling_lo[1]))

    t_hi, x_hi, coupling_hi = _hi_out.value
    c_hi = _hi_out.time

    for i in 1:10
        idx = findfirst(==(10*i), selectdim(x_hi, 1, 4))
        y_hi[i] = t_hi[idx]
    end

    return c_hi, coupling_hi
end

function make_abc(ϵ, data=DATA)
    ℓ_ABC = function (θ, y)
        d = euclidean(y, data)
        return float(d<ϵ)
    end
    return ℓ_ABC
end

function make_bsl(data=DATA)
    ℓ_bsl = function (θ, y)
        mu = mean.(eachrow(y))
        sigma = cov(y')
        return pdf(MvNormal(mu, sigma), data)
    end
    return ℓ_bsl
end

q = product_distribution([
    Uniform(10.0, 100.0),
    Uniform(10.0, 100.0),
    Uniform(0.1, 10.0),
])
prior = q

G = θ -> θ[3]

##############
# BURN IN
##############

struct MFIteration
    θ::Vector{Float64}
    # Low fidelity
    y_lo::Vector{Float64}
    c_lo::Float64
    L_lo::Float64
    # Choose num high fidelity
    m::Int64
    μ::Float64
    # High fidelity
    Σc_hi::Float64
    ΣL_hi::Float64
    ΣL_hi2::Float64
    # Correction terms
    ΣL_mf::Float64
    ΣL_mf2::Float64
end
MFSample = StructVector{MFIteration}
MFWeight(p::MFIteration) = p.L_lo + (p.ΣL_mf)/p.μ
is_output(p::MFIteration) = !iszero(MFWeight(p))
ImportanceWeight(p::MFIteration; prior=prior, q=q) = exp(logpdf(prior, p.θ) - logpdf(q, p.θ))

function StatsBase.mean(G, S::MFSample)
    Gθ = G.(S.θ)
    w_mf = MFWeight.(S)
    return mean(Gθ, Weights(w_mf))
end
ImportanceWeight(S::MFSample; prior=prior, q=q) = ImportanceWeight.(S; prior=prior, q=q)
function Δ(G, S::MFSample; kwargs...)
    Ḡ = mean(G, S)
    Delta = function (p::MFIteration)
        ImportanceWeight(p; kwargs...)*abs(G(p.θ) - Ḡ)
    end
    return Delta
end

function MFABCIteration(ℓ_lo, ℓ_hi, μfun=(θ, y_lo)->1.0, M=μ->Poisson(μ))
    θ = rand(q)
    y_lo = zeros(10)
    y_hi = zeros(10)
    Σc_hi = 0.0
    ΣL_hi = 0.0
    ΣL_hi2 = 0.0
    ΣL_mf = 0.0
    ΣL_mf2 = 0.0

    # Low fidelity
    c_lo, cpl = f_lo(y_lo, θ)
    L_lo = ℓ_lo(θ, y_lo)

    # Choose num high fidelity
    μ = μfun(θ, y_lo)
    m::Int64 = rand(M(μ))

    # Do high fidelity
    for _i in 1:m
        c_hi_i, _ = f_hi(y_hi, θ, cpl)
        L_hi_i = ℓ_hi(θ, y_hi)

        Σc_hi += c_hi_i
        ΣL_hi += L_hi_i
        ΣL_hi2 += L_hi_i^2
        ΣL_mf += (L_hi_i - L_lo)
        ΣL_mf2 += (L_hi_i - L_lo)^2
    end

    return MFIteration(θ, y_lo, c_lo, L_lo, m, μ, Σc_hi, ΣL_hi, ΣL_hi2, ΣL_mf, ΣL_mf2)
end

function MFABCBatch(N::Int; ϵ::Real, μfun=(θ, y_lo)->1.0, M=μ->Poisson(μ))
    ℓ_lo = make_abc(ϵ)
    ℓ_hi = make_abc(ϵ)
    return StructVector(map(1:N) do i
        MFABCIteration(ℓ_lo, ℓ_hi, μfun, M)
    end)
end
MFABCBurnin(N; ϵ) = MFABCBatch(N; ϵ=ϵ, μfun=(θ, y_lo)->1, M=μ -> DiscreteUniform(μ,μ))

function MFBSLIteration(ℓ_lo, ℓ_hi, μfun=(θ, y_lo)->1.0, M=μ->Poisson(μ); n_bsl = 100)
    θ = rand(q)
    y_lo = zeros(10, n_bsl)
    y_hi = zeros(10, n_bsl)
    Σc_hi = 0.0
    ΣL_hi = 0.0
    ΣL_hi2 = 0.0
    ΣL_mf = 0.0
    ΣL_mf2 = 0.0

    # Low fidelity
    c_lo_batch = map(eachcol(y_lo)) do y_lo_i
        f_lo(y_lo_i, θ)
    end
    L_lo = ℓ_lo(θ, y_lo)
    c_lo = sum(tup->tup[1], c_lo_batch)
    cpl = (tup[2] for tup in c_lo_batch)

    # Choose num high fidelity
    μ = μfun(θ, vec(y_lo))
    m::Int64 = rand(M(μ))

    # Do high fidelity
    for _i in 1:m
        c_hi_i = 0.0
        for (y_hi_c, cpl_c) in zip(eachcol(y_hi), cpl)
            c_hi_i_c, _ = f_hi(y_hi_c, θ, cpl_c)
            c_hi_i += c_hi_i_c
        end
        L_hi_i = ℓ_hi(θ, y_hi)

        Σc_hi += c_hi_i
        ΣL_hi += L_hi_i
        ΣL_hi2 += L_hi_i^2
        ΣL_mf += (L_hi_i - L_lo)
        ΣL_mf2 += (L_hi_i - L_lo)^2
    end

    return MFIteration(θ, vec(y_lo), c_lo, L_lo, m, μ, Σc_hi, ΣL_hi, ΣL_hi2, ΣL_mf, ΣL_mf2)
end

function MFBSLBatch(N::Int; μfun=(θ, y_lo)->1.0, M=μ->Poisson(μ), n_bsl=100)
    ℓ_lo = make_bsl()
    ℓ_hi = make_bsl()
    return StructVector(map(1:N) do i
        MFBSLIteration(ℓ_lo, ℓ_hi, μfun, M; n_bsl=n_bsl)
    end)
end

struct MFNode{L,R}
    featid::Int64
    featval::Float64
    left::L
    right::R
end
mutable struct MFLeaf
    logμ::Float64
    c_k::Float64
    V_k_G0::Float64
    V_k_G1::Float64
    V_k_G2::Float64
    function MFLeaf()
        return new(0.0, 0.0, 0.0, 0.0, 0.0)
    end
end
mutable struct MFMean{R<:MFNode}
    root::R
    r::Int64
    c̄_lo::Float64
    V_mf_G0::Float64
    V_mf_G1::Float64
    V_mf_G2::Float64
    Ḡ_num::Float64
    Ḡ_den::Float64
    function MFMean(decision_tree::DecisionTree.Node)
        root = make_MF_vertex(decision_tree)
        R = typeof(root)
        μ = new{R}(root)
        μ.r = 0
        μ.c̄_lo = 0.0
        μ.V_mf_G0 = 0.0
        μ.V_mf_G1 = 0.0
        μ.V_mf_G2 = 0.0
        μ.Ḡ_num = 0.0
        μ.Ḡ_den = 0.0
        return μ
    end
end

include("mean_fun.jl")

function MFMean(S::MFSample, G=G)
    Ḡ = mean(G, S)

    idx = findall(p->p.m>0, S)
    N = length(idx)
    dim_features = length(S.θ[1]) + length(S.y_lo[1])

    features = zeros(N, dim_features)
    for (_n, _i) in enumerate(idx)
        view(features, _n, 1:3) .= S.θ[_i]
        view(features, _n, 4:dim_features) .= S.y_lo[_i]
    end

    # Create a tree from the burn-in sample
    Stree = S[idx]
    Δ = ImportanceWeight(Stree) .* abs.(G.(Stree.θ) .- Ḡ)

    mu_star = @. Δ * sqrt(Stree.ΣL_mf2/Stree.Σc_hi)
    mu_star ./= maximum(mu_star)
    model = build_tree(mu_star, features)

    μ = MFMean(model)
    increment!(μ, S, G)
    return μ
end

function increment!(μ::MFMean, S::MFSample, G=G)
    for p in S
        increment!(μ, p, G)
    end
    return nothing
end
function increment!(μ::MFMean, p::MFIteration, G=G)
    μ.r += 1

    # Increment cost
    μ.c̄_lo += p.c_lo

    # Increment mean
    Gθ = G(p.θ)
    w_mf = p.L_lo + (p.ΣL_mf / p.μ)
    μ.Ḡ_num += w_mf*Gθ
    μ.Ḡ_den += w_mf

    # Increment V_mf components
    V_mf_inc = (p.ΣL_hi^2 - p.ΣL_hi2) / (p.μ^2)
    μ.V_mf_G0 += V_mf_inc
    μ.V_mf_G1 += V_mf_inc * Gθ
    μ.V_mf_G2 += V_mf_inc * Gθ * Gθ

    # Find leaf
    leaf = get_MF_leaf(μ, p)

    # Increment cost
    leaf.c_k += (p.Σc_hi / p.μ)

    # Increment V_k components
    V_k_inc = p.ΣL_mf2 / p.μ
    leaf.V_k_G0 += V_k_inc
    leaf.V_k_G1 += V_k_inc * Gθ
    leaf.V_k_G2 += V_k_inc * Gθ * Gθ

    return nothing
end

function descend_gradient!(μ::MFMean, n=1, δ=1.0)
    Ḡ = μ.Ḡ_num / μ.Ḡ_den
    r = μ.r

    f_c = leaf -> exp(leaf.logμ) * leaf.c_k / r
    f_V = leaf -> exp(-leaf.logμ) * ((leaf.V_k_G0 * Ḡ^2) - 2*(leaf.V_k_G1 * Ḡ) + (leaf.V_k_G2)) / r

    c̄_lo = μ.c̄_lo / r
    V_mf = ((μ.V_mf_G0 * Ḡ^2) - 2*(μ.V_mf_G1 * Ḡ) + (μ.V_mf_G2)) / r

    for _ in 1:n
        c = c̄_lo + sum(f_c, leaves(μ))
        V = V_mf + sum(f_V, leaves(μ))

        for leaf in leaves(μ)
            leaf.logμ -= δ * (
                f_c(leaf)*V - f_V(leaf)*c
            )
        end
    end

    ν = map(leaves(μ)) do leaf
        sqrt((((leaf.V_k_G0 * Ḡ^2) - 2*(leaf.V_k_G1 * Ḡ) + (leaf.V_k_G2))/V_mf) / (leaf.c_k / c̄_lo))
    end
    return ν
end

######

function do_abc(N=10_000; ϵ=5, G=G, N_batch=1000)
    ℓ = make_abc(ϵ)

    θ = rand(q, N_batch)
    y = zeros(10, N_batch)
    n = 0
    c = 0.0
    w = Float64[]
    Gθ = Float64[]
    while n<N
        c_batch = f_hi.(eachcol(y), eachcol(θ))
        w_batch = ℓ.(eachcol(θ), eachcol(y))
        Gθ_batch = G.(eachcol(θ))

        c += sum(tup->tup[1], c_batch)
        n += N_batch
        append!(w, w_batch)
        append!(Gθ, Gθ_batch)

        Distributions.rand!(q, θ)
    end

    Ḡ = sum(Gθ.*w)/sum(w)
    σ2 = sum((w .* (Gθ .- Ḡ)).^2) / sum(w)^2

    return Ḡ, σ2, c
end

function do_bsl(N=10_000; n_bsl=100, G=G, N_batch=1000)
    ℓ = make_bsl()

    θ = rand(q, N_batch)
    y = zeros(10, n_bsl, N_batch)
    n = 0
    c = 0.0
    w = Float64[]
    Gθ = Float64[]

    w_batch = zeros(N_batch)
    Gθ_batch = zeros(N_batch)
    while n<N
        for i in 1:N_batch
            for j in 1:n_bsl
                _t, _ = f_hi(view(y, :, j, i), view(θ, :, i))
                c += _t
            end

            w_batch[i] = ℓ(view(θ, :, i), view(y, :, :, i))
            Gθ_batch[i] = G(view(θ, :, i))
        end
        n += N_batch
        append!(w, w_batch)
        append!(Gθ, Gθ_batch)

        Distributions.rand!(q, θ)
    end

    Ḡ = sum(Gθ.*w)/sum(w)
    σ2 = sum((w.*(Gθ .- Ḡ)).^2) / sum(w)^2

    return Ḡ, σ2, c
end

function do_mfabc(N=100_000; ϵ=5, N₀=10000, G=G, N_batch=1000, GD_batch=1000, δ=0.01)
    S = MFABCBatch(N₀, ϵ=ϵ)
    μ = MFMean(S, G)
    leaf_set = leaves(μ)
    leaf_traj = zeros(N, length(μ.root))
    ν_traj = zeros(N, length(μ.root))

    ctr = N₀
    while length(S)<N
        ν = descend_gradient!(μ, GD_batch, δ)
        dS = MFABCBatch(N_batch, ϵ=ϵ, μfun = μ)
        increment!(μ, dS, G)
        append!(S, dS)
        leaf_traj_slice = selectdim(leaf_traj, 1, ctr .+ (1:N_batch))
        for (leaf_traj_slice_col, leaf) in zip(eachcol(leaf_traj_slice), leaf_set)
            leaf_traj_slice_col .= leaf.logμ
        end
        for i in 1:N_batch
            selectdim(ν_traj, 1, ctr+i) .= ν
        end
        ctr += N_batch
    end
    return S, μ, leaf_traj, ν_traj
end

function do_mfbsl(N=10_000; n_bsl=100, N₀=1000, G=G, N_batch=1000, GD_batch=1000, δ=0.01)
    S = MFBSLBatch(N₀, n_bsl=n_bsl)
    μ = MFMean(S, G)
    leaf_set = leaves(μ)
    leaf_traj = zeros(N, length(μ.root))
    ν_traj = zeros(N, length(μ.root))

    ctr = N₀
    while length(S)<N
        ν = descend_gradient!(μ, GD_batch, δ)
        dS = MFBSLBatch(N_batch, μfun = μ, n_bsl=n_bsl)
        increment!(μ, dS, G)
        append!(S, dS)
        leaf_traj_slice = selectdim(leaf_traj, 1, ctr .+ (1:N_batch))
        for (leaf_traj_slice_col, leaf) in zip(eachcol(leaf_traj_slice), leaf_set)
            leaf_traj_slice_col .= leaf.logμ
        end
        for i in 1:N_batch
            selectdim(ν_traj, 1, ctr+i) .= ν
        end
        ctr += N_batch
    end
    return S, μ, leaf_traj, ν_traj
end

max_stars_repo_name: tpprescott/mf-lf

max_issues_repo_name: tpprescott/mf-lf

max_forks_repo_name: tpprescott/mf-lf


text:
using Vofinit
using Test

function test_getcelltype()
    x0 = Cdouble.((0, 0, 0))
    h0 = Cdouble.((1, 1, 1))

    getcelltype(x0, h0) do x...
        sum(x .* x) - 0.25
    end
end

function test_getcc()
    x0 = Cdouble.((0, 0, 0))
    h0 = Cdouble.((1, 1, 1))
    xex = Cdouble[0, 0, 0, 0]

    getcc(x0, h0, xex) do x...
        sum(x .* x) - 0.25
    end
end

@testset "Vofinit.jl" begin
    @test test_getcelltype() == -1
    @test isapprox(test_getcc(), π / 48, rtol = 5e-3)
end


text:
# Julia wrapper for header: /Users/dls/jl/openssl-1.0.1e/include/openssl/sha.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 SHA_Init (Ptr{SHA_CTX},) libcrypto
@c Int32 SHA_Update (Ptr{SHA_CTX}, Ptr{None}, Csize_t) libcrypto
@c Int32 SHA_Final (Ptr{Uint8}, Ptr{SHA_CTX}) libcrypto
@c Ptr{Uint8} SHA (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto
@c None SHA_Transform (Ptr{SHA_CTX}, Ptr{Uint8}) libcrypto
@c Int32 SHA1_Init (Ptr{SHA_CTX},) libcrypto
@c Int32 SHA1_Update (Ptr{SHA_CTX}, Ptr{None}, Csize_t) libcrypto
@c Int32 SHA1_Final (Ptr{Uint8}, Ptr{SHA_CTX}) libcrypto
@c Ptr{Uint8} SHA1 (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto
@c None SHA1_Transform (Ptr{SHA_CTX}, Ptr{Uint8}) libcrypto
@c Int32 SHA224_Init (Ptr{SHA256_CTX},) libcrypto
@c Int32 SHA224_Update (Ptr{SHA256_CTX}, Ptr{None}, Csize_t) libcrypto
@c Int32 SHA224_Final (Ptr{Uint8}, Ptr{SHA256_CTX}) libcrypto
@c Ptr{Uint8} SHA224 (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto
@c Int32 SHA256_Init (Ptr{SHA256_CTX},) libcrypto
@c Int32 SHA256_Update (Ptr{SHA256_CTX}, Ptr{None}, Csize_t) libcrypto
@c Int32 SHA256_Final (Ptr{Uint8}, Ptr{SHA256_CTX}) libcrypto
@c Ptr{Uint8} SHA256 (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto
@c None SHA256_Transform (Ptr{SHA256_CTX}, Ptr{Uint8}) libcrypto
@c Int32 SHA384_Init (Ptr{SHA512_CTX},) libcrypto
@c Int32 SHA384_Update (Ptr{SHA512_CTX}, Ptr{None}, Csize_t) libcrypto
@c Int32 SHA384_Final (Ptr{Uint8}, Ptr{SHA512_CTX}) libcrypto
@c Ptr{Uint8} SHA384 (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto
@c Int32 SHA512_Init (Ptr{SHA512_CTX},) libcrypto
@c Int32 SHA512_Update (Ptr{SHA512_CTX}, Ptr{None}, Csize_t) libcrypto
@c Int32 SHA512_Final (Ptr{Uint8}, Ptr{SHA512_CTX}) libcrypto
@c Ptr{Uint8} SHA512 (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto
@c None SHA512_Transform (Ptr{SHA512_CTX}, Ptr{Uint8}) libcrypto


text:
using ApproxFun, LazyArrays, FillArrays, LinearAlgebra, SpecialFunctions, Test
    import ApproxFun: interlace, Multiplication, ConstantSpace, PointSpace,
                        ArraySpace, testblockbandedoperator, blocklengths, ∞,
                        testraggedbelowoperator, Vec

@testset "Vector" begin
    @testset "Construction" begin
        f = Fun(x->[1.,0.])
        @test f(0.) ≈ [1.,0.]
        @test f[end] ≈ 0
        @test f[firstindex(f)] ≈ 1
    end

    @testset "Broadcast" begin
        x = Fun()

        @test x .+ [1,2] ≈ [x+1,x+2]
        @test [1,2] .+ x ≈ [x+1,x+2]

        @test x.*[1,2] ≈ [x,2x]
        @test [1,2].*x ≈ [x,2x]
        @test (x*[1,2])(0.1) ≈ [0.1,0.2]
        @test ([1,2]*x)(0.1) ≈ [0.1,0.2]

        @test (x*[1,2])(0.1) ≈ [0.1,0.2]
        @test ([1,2]*x)(0.1) ≈ [0.1,0.2]

        @test (x*Fun([1,2]))(0.1) ≈ [0.1,0.2]
        @test (Fun([1,2])*x)(0.1) ≈ [0.1,0.2]

        @test (Fun(x->[1., 2.]) + [2, 2])(0.) ≈ [3., 4.]
    end

    @testset "Vector*Vector{Fun}" begin
        f = Fun(x->[exp(x),cos(x)])
        @test f(0.1) ≈ [exp(0.1),cos(0.1)]
        @test ([1 2]*f)(0.1) ≈ [1 2]*f(0.1)
        @test (f*3)(0.1) ≈ f(0.1)*3
        @test (3*f)(0.1) ≈ f(0.1)*3

        @test_broken transpose(f)*[1,2] ≈ transpose(f(0.1))*[1,2]

        @test norm(f) ≈ sqrt(sinh(2)+1+cos(1)sin(1))
        @test norm(f,2) ≈ sqrt(sinh(2)+1+cos(1)sin(1))
        @test norm(f,1) ≈ sqrt(sinh(2))+sqrt(1+cos(1)sin(1))
    end

    @testset "Chebyshev Vector" begin
        a = [1 2; 3 4]
        f = Fun(x->[exp(x),cos(x)])
        @test (a*f)(0.1) ≈ [exp(0.1)+2cos(0.1); 3exp(0.1)+4cos(0.1)]
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f
    end

    @testset "Chebyshev Matrix" begin
        a = [1 2; 3 4]
        m = Fun(x->[exp(x) cos(x); sin(x) airyai(x)])
        @test (a*m)(0.1) ≈ [exp(0.1)+2sin(0.1) cos(0.1)+2airyai(0.1);
                            3exp(0.1)+4sin(0.1) 3cos(0.1)+4airyai(0.1)]
        @test (a*m)(0.1) ≈ a*m(0.1)
        @test (m*a)(0.1) ≈ m(0.1)*a
        @test Fun(a)*m   ≈ a*m
        @test m*Fun(a)   ≈ m*a
        @test Fun(a*Array(m))   ≈ a*m
        @test Fun(Array(m)*a)   ≈ m*a

        @test (a+m)(0.1) ≈ a+m(0.1)
        @test (m+a)(0.1) ≈ m(0.1)+a

        @test (m+I)(0.1) ≈ m(0.1)+I
    end

    @testset "CosSpace Vector" begin
        a = [1 2; 3 4]
        f = Fun(θ->[1,cos(θ)],CosSpace())
        @test (a*f)(0.1) ≈ [1+2cos(0.1); 3+4cos(0.1)]
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f
    end

    @testset "CosSpace Matrix" begin
        a = [1 2; 3 4]
        m = Fun(θ->[1 cos(θ); cos(2θ) cos(cos(θ))],CosSpace())
        @test (a*m)(0.1) ≈ a*m(0.1)
        @test (m*a)(0.1) ≈ m(0.1)*a
        @test Fun(a)*m   ≈ a*m
        @test Fun(a*Array(m))   ≈ a*m

        @test (a+m)(0.1) ≈ a+m(0.1)
        @test (m+a)(0.1) ≈ m(0.1)+a

        @test (m+I)(0.1) ≈ m(0.1)+I
    end

    @testset "SinSpace Vector" begin
        a = [1 2; 3 4]
        f = Fun(θ->[sin(θ),sin(2θ)],SinSpace())
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f

        @test all(sp -> sp isa SinSpace, space(a*f).spaces)
    end

    @testset "CosSpace Matrix" begin
        a = [1 2; 3 4]
        m = Fun(θ->[sin(3θ) sin(θ); sin(2θ) sin(sin(θ))],SinSpace())
        @test (a*m)(0.1) ≈ a*m(0.1)
        @test (m*a)(0.1) ≈ m(0.1)*a
        @test Fun(a)*m   ≈ a*m
        @test Fun(a*Array(m))   ≈ a*m

        @test all(sp -> sp isa SinSpace, space(a*m).spaces)

        @test (a+m)(0.1) ≈ a+m(0.1)
        @test (m+a)(0.1) ≈ m(0.1)+a

        @test (m+I)(0.1) ≈ m(0.1)+I
    end

    @testset "Fourier Vector" begin
        a = [1 2; 3 4]
        f = Fun(θ->[sin(θ),sin(2θ)],Fourier())
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f
        @test norm(f) ≈ sqrt(2π)
        @test norm(f,2) ≈ sqrt(2π)
    end

    @testset "Matrix{Fun}*Matrix{Fun}" begin
    # note that 2x2 and 3x3 mult are special cases
        x=Fun()
        A = [x x; x x]

        @test A isa Fun

        @test norm(map(norm,A*A-[2x^2 2x^2; 2x^2 2x^2])) <eps()
        @test norm(map(norm,A*Array(A)-[2x^2 2x^2; 2x^2 2x^2])) < eps()
        @test norm(map(norm,Array(A)*A-[2x^2 2x^2; 2x^2 2x^2])) < eps()
        @test norm(map(norm,Array(A)*Array(A)-Array([2x^2 2x^2; 2x^2 2x^2]))) < eps()

        @test [x;x] == A[:,1]
        @test [[x;2x] [3x;4x]] == [x 3x; 2x 4x]
        @test [x x] == A[1:1,:]

        A = fill(x,3,3)
        @test norm(map(norm,A^3-fill(9x^3,3,3))) <eps()
        @test norm((A^2*Fun(A)-fill(9x^3,3,3)).coefficients) < eps()
        @test norm((A*Fun(A)^2-fill(9x^3,3,3)).coefficients) < eps()
        @test norm((Fun(A)^3-fill(9x^3,3,3)).coefficients) < eps()

        A = fill(x,4,4)
        @test norm(map(norm,A^2-fill(4x^2,4,4))) <eps()
        @test norm((A*Fun(A)-fill(4x^2,4,4)).coefficients) < eps()
        @test norm((Fun(A)*A-fill(4x^2,4,4)).coefficients) < eps()
        @test norm((Fun(A)^2-fill(4x^2,4,4)).coefficients) < eps()
    end

    @testset "Vector ODE" begin
        d=ChebyshevInterval()
        D=Derivative(d);
        B=ldirichlet();
        Bn=lneumann();
        A=[B 0;
           0 B;
           D-I 2I;
           0I D+I];

        f=Fun(x->[exp(x),cos(x)],d)

        b=Any[0.,0.,f...]
        f1,f2=vec(f)
        @time u=A\b
        u1=vec(u)[1];u2=vec(u)[2];

        @test norm(u1'-u1+2u2-f1)<10eps()
        @test norm(u2'+u2-f2)<10eps()

        A=[B 0;
           Bn 0;
           0 B;
           D^2-I 2I;
           0 D+I];

        b=Any[0.,0.,0.,f...]

        @time u=A\b
        u1=vec(u)[1];u2=vec(u)[2];

        @test norm(differentiate(u1,2)-u1+2u2-f1) < 4eps()
        @test norm(u2'+u2-f2) < 4eps()
    end

    @testset "Matrix exponential" begin
        n=4
        d=fill(Interval(0.,1.),n)
        B=Evaluation(d,0.)
        D=Derivative(d)
        A=rand(n,n)
        L=[B;D-A]
        @time u=L\Matrix(I,2n,n)
        @test norm(u(1.)-exp(A))<eps(1000.)

        n=4
        d=fill(Interval(0.,1.),n)
        B=Evaluation(d,0.)
        D=Derivative(d)
        A=rand(n,n)
        L=[B;D-A]
        @time u=L\Matrix(I,2n,2)
        @test norm(u(1.)-exp(A)[:,1:2])<eps(1000.)
    end

    @testset "Multiplication" begin
        d = ChebyshevInterval()
        t=Fun(identity,d)
        f = Fun([t^2, sin(t)])
        @test norm(((Derivative(space(f))*f)-Fun(t->[2t,cos(t)])).coefficients)<100eps()
        @test norm((([1 2;3 4]*f)-Fun(t->[t^2+2sin(t),3t^2+4sin(t)])).coefficients)<100eps()
    end

    @testset "Two circles" begin
        Γ = Circle() ∪ Circle(0.5)

        f = Fun(z -> in(z,component(Γ,2)) ? 1 : z,Γ)
        @test f(exp(0.1im)) ≈ exp(0.1im)
        @test f(0.5exp(0.1im)) ≈ 1

        G = Fun(z -> in(z,component(Γ,2)) ? [1 -z^(-1); 0 1] :
                                            [z 0; 0 z^(-1)], Γ);

        @test G(exp(0.1im)) ≈ [exp(0.1im) 0 ; 0 exp(-0.1im)]
        @test G(0.5exp(0.1im)) ≈ [1 -2exp(-0.1im) ; 0 1]

        G1=Fun(Array(G)[:,1])

        @test G1(exp(0.1im)) ≈ [exp(0.1im),0.]
        @test G1(0.5exp(0.1im)) ≈ [1,0.]

        M = Multiplication(G, space(G1))
        testblockbandedoperator(M)

        for z in (0.5exp(0.1im),exp(0.2im))
            @test G[1,1](z) ≈ G[1](z)
            @test (M.op.ops[1,1]*G1[1])(z) ≈ M.f[1,1](z)*G1[1](z)
            @test (M.op.ops[2,1]*G1[1])(z) ≈ M.f[2,1](z)*G1[1](z)
            @test (M.op.ops[1,2]*G1[2])(z) ≈ M.f[1,2](z)*G1[2](z)
            @test (M.op.ops[2,2]*G1[2])(z) ≈ M.f[2,2](z)*G1[2](z)
        end

        u = M*G1
        @test norm(u(exp(.1im))-[exp(.2im),0])<100eps()
        @test norm(u(.5exp(.1im))-[1,0])<100eps()
    end

    @testset "Circle" begin
        G = Fun(z->[-1 -3; -3 -1]/z +
                   [ 2  2;  1 -3] +
                   [ 2 -1;  1  2]*z, Circle())

        @test G[1,1](exp(0.1im)) == G(exp(0.1im))[1,1]

        F̃ = Array((G-I)[:,1])
        F = (G-I)[:,1]

        @test Fun(F) ≡ F

        @test F(exp(0.1im)) ≈ [-exp(-0.1im)+1+2exp(0.1im);-3exp(-0.1im)+1+1exp(0.1im)]
        @test Fun(F̃,space(F))(exp(0.1im)) ≈ [-exp(-0.1im)+1+2exp(0.1im);-3exp(-0.1im)+1+1exp(0.1im)]

        @test coefficients(F̃,space(F)) == F.coefficients
        @test Fun(F̃,space(F)) == F

        @test F == Fun(vec(F),space(F))

        @test inv(G(exp(0.1im))) ≈ inv(G)(exp(0.1im))

        @test Fun(Matrix(I,2,2),space(G))(exp(0.1im)) ≈ Matrix(I,2,2)
        @test Fun(I,space(G))(exp(0.1im)) ≈ Matrix(I,2,2)
    end

    @testset "Conversion" begin
        f=Fun(t->[cos(t) 0;sin(t) 1],-π..π)
        g=Fun(f,Space(PeriodicSegment(-π,π)))
        @test g(.1) ≈ f(.1)

        a = ArraySpace(JacobiWeight(1/2,1/2, Chebyshev()), 2)
        b = ArraySpace(JacobiWeight(1/2,1/2, Ultraspherical(1)), 2)
        C = Conversion(a, b)

        f = Fun(a, rand(10))
        @test f(0.1) ≈ (C*f)(0.1)

        a = ArraySpace(JacobiWeight(1/2,1/2, Chebyshev()), 2,3)
        b = ArraySpace(JacobiWeight(1/2,1/2, Ultraspherical(1)), 2,3)
        C = Conversion(a, b)

        f = Fun(a, rand(10))
        @test f(0.1) ≈ (C*f)(0.1)
    end

    @testset "Interlace" begin
        S1 = Chebyshev()^2
        S2 = Chebyshev()
        TS = ArraySpace([ConstantSpace(),S1,ConstantSpace(),S2,PointSpace([1.,2.])])
        f = Fun(TS,collect(1:13))
        @test f[1] == Fun(TS[1],[1.])
        @test f[2] == Fun(TS[2],[2.,7.,8.,10.,11.,12.])
        @test f[3] == Fun(TS[3],[3.])
        @test f[4] == Fun(TS[4],[4.,9.,13.])
        @test f[5] == Fun(TS[5],[5.,6.])
    end

    @testset "Operator * Matrix" begin
        f = [Fun(exp) Fun(cos)]
        @test f(0.1) ≈ [exp(0.1) cos(0.1)]

        D=Derivative()
        u=D*f
        @test u(0.1) ≈ [exp(0.1) -sin(0.1)]
        u=D*Array(f)
        @test u(0.1) ≈ [exp(0.1) -sin(0.1)]
    end

    @testset "Floquet" begin
        # TODO: fix when IntervalSets.jl is fixed
        T = Float64(π, RoundUp); a=0.15
        t = Fun(identity,0..T)
        d=domain(t)
        D=Derivative(d)
        B=ldirichlet(d)

        B_row = [D             -I  0I            0I]

        f=Fun(exp,d)
        @test norm((B_row*[f;f;f;f])[1]) ≤ 1000eps()
        @test B_row isa ApproxFun.MatrixInterlaceOperator

        @test size([B_row;B_row].ops) == (2,4)

        @test norm(([B_row;B_row]*[f;f;f;f])[1]) ≤ 1000eps()

        n=4
        Dg = Operator(diagm(0 => fill(ldirichlet(d),n)))

        @test Dg isa ApproxFun.MatrixInterlaceOperator
        @test size([Dg; B_row].ops) == (5,4)

        L = [Dg; B_row]
        testraggedbelowoperator(L)
        @test ([Dg; B_row]*[f;f;f;f])(0.1) ≈ [fill(1.0,4);0]

        @test hcat(Dg).ops == Dg.ops

        B_row2 = [(2+a*cos(2t))   D  -I            0I]
        @test ([B_row;B_row2]*[f;f;f;f])(0.1) ≈ [0.,(2+a*cos(2*0.1))*f(0.1) + f'(0.1) - f(0.1)]

        A=[ Operator(diagm(0 => fill(ldirichlet(d),n)));
            D             -I  0I            0I;
           (2+a*cos(2t))   D  -I            0I;
           0I             0I   D            -I;
           -I             0I  (2+a*cos(2t))  D]

        Φ = A\Matrix(I,2n,n);

        @test Φ(π) ≈ [-0.170879 -0.148885 -0.836059 0.265569;
                 0.732284 -0.170879 -0.612945 -0.836059;
                 -0.836059 0.265569 -0.170879 -0.148885;
                 -0.612945 -0.836059 0.732284 -0.170879] atol=1E-3
    end

    @testset "diagonal" begin
        x = Fun(0.0..1)
        @test isambiguous(domain(zero(typeof(x))))
        @test diagm(0 => [x,x]) * [1,2] == [x,2x]
    end
end


text:
module Traffic

using LightGraphs, SimpleWeightedGraphs
using Distributions
using JuMP
import Ipopt

export qʰALLTime, qᵍALLTime, solve_TA, Nₕ

##########################################
# SiouxFalls traffic info
##########################################
include("readTrafficInfo.jl")
const network_data_file = "./data/SiouxFalls_net.tntp";
const trip_table_file = "./data/SiouxFalls_trips.tntp";
td = load_ta_network("SiouxFalls", network_data_file, trip_table_file);
t₀ = td.free_flow_time;
cₐ = td.capacity;
n = td.number_of_nodes;
b = td.number_of_links;
roads = [start_node => end_node for (start_node, end_node) in zip(td.init_node, td.term_node)];
road2roadidx = Dict(roads[i]=>i  for i in 1:length(roads))

##########################################
# additional travel info
##########################################
BRP_int(xₐ, t₀, cₐ) = t₀ * (xₐ + 0.03 * xₐ^5  / cₐ^4) #  integral of BRP

# traffic graph
g = SimpleWeightedDiGraph(n);
for (start_node, end_node, link_length) in zip(td.init_node, td.term_node, td.link_length)
    add_edge!(g, start_node, end_node, link_length)
end
# Yen's algorithm to find k shortest paths
# - there is a bug in `rem_edge!` of SimpleWeightedGraphs.jl
# - so convert SimpleWeightedGraphs to a SimpleDiGraph
# https://github.com/JuliaGraphs/LightGraphs.jl/issues/1505
# https://github.com/JuliaGraphs/SimpleWeightedGraphs.jl/issues/66
K = 10 # K shortest path in each OD
available_path_GV = Dict()
for start_node in 1:n
    for end_node in 1:n
        if start_node ≠ end_node
            path_state = yen_k_shortest_paths(SimpleDiGraph(g), start_node, end_node, weights(g), K)
            temp_via_roads_idx = []
            for path in path_state.paths
                via_roads = [path[i-1] => path[i] for i in 2:length(path)]
                via_roads_idx = [road2roadidx[road] for road in via_roads]
                push!(temp_via_roads_idx, via_roads_idx)
            end
            available_path_GV[start_node=>end_node] = temp_via_roads_idx
        end
    end
end

##########################################
# additional hydrogen refueling station info
##########################################
Eₕ = 10; # hydrogen required per vehicle: 10kg
Nₕ = 8; # Number of hydrogen refueling stations

# ρⱼ = ones(Nₕ) * 10; # price of hydrogen kg/$
ω = 0.5 # price of time $/minute
#  hydrogen stations are on these roads

station_placement = [16=>8, 3=>4, 6=>5, 9=>10, 22=>21, 12=>11, 17=>19, 23=>14];
station2roadidx = [road2roadidx[road] for road in station_placement];

# find K paths that contain at least one HRS
available_path_FCEV = Dict()
available_path_via_stationidx_FCEV = Dict()
for start_node in 1:n
    for end_node in 1:n
        if start_node ≠ end_node
            path_state = yen_k_shortest_paths(SimpleDiGraph(g), start_node, end_node, weights(g), 200)
            temp_via_roads_idx = []
            temp_via_stationidx = []
            # find K paths
            for path in path_state.paths
                via_roads = [path[i-1] => path[i] for i in 2:length(path)]
                if (via_roads ∩ station_placement) ≠ []
                    via_roads_idx = [road2roadidx[road] for road in via_roads]
                    push!(temp_via_roads_idx, via_roads_idx)
                    push!(temp_via_stationidx, [findfirst(x->x==station, station_placement)
                        for station in via_roads ∩ station_placement])
                    if length(temp_via_roads_idx) >= K
                        break
                    end
                end
            end
            available_path_FCEV[start_node=>end_node] = temp_via_roads_idx
            available_path_via_stationidx_FCEV[start_node=>end_node] = temp_via_stationidx
        end
    end
end

# from here, all dict keys are converted to the "index" in OD
# convenient to use in solve_TA
OD = [start_node=>end_node for start_node in 1:n for end_node in 1:n if start_node ≠ end_node]
# shape (OD_idx, path_idx, road)
pathʰ = [available_path_FCEV[OD_pairs] for OD_pairs in OD]
pathᵍ = [available_path_GV[OD_pairs] for OD_pairs in OD]
pathʰ_via_stationidx = [available_path_via_stationidx_FCEV[OD_pairs] for OD_pairs in OD]

# from road index to know all (OD, path) through it
road_idx2pathʰidx = [[] for _ in roads]
road_idx2pathᵍidx = [[] for _ in roads]
for (OD_idx, paths) in enumerate(pathʰ)
    for (path_idx, path) in enumerate(paths)
        for road_idx in path
            push!(road_idx2pathʰidx[road_idx], (OD=OD_idx, path=path_idx))
        end
    end
end
for (OD_idx, paths) in enumerate(pathᵍ)
    for (path_idx, path) in enumerate(paths)
        for road_idx in path
            push!(road_idx2pathᵍidx[road_idx], (OD=OD_idx, path=path_idx))
        end
    end
end

##########################################
# additional traffic demand info
##########################################
# assume the traffic demand is normal distribution
# assume 8:00 and 20:00 is the peak of the distribution
# its value: base_travel_demand * adjust_factor
# 2 a.m. to 2 p.m. : Normal(8, 3) * base_travel_demand * adjust_factor
# 2 p.m. to 2 a.m. : - Normal(20, 3) * base_travel_demand * adjust_factor
# the reverse traffic demand is the transpose of demand matrix
base_travel_demand = td.travel_demand
adjust_factor = 1
# travel demand for all FCEVs and GVs
travel_demand = [] # Vector{Matrix} (time, start_node, end_node)
T = 24
for t in 1:T
    if t < 2
        push!(travel_demand, pdf(Normal(20, 5), 24+t) / pdf(Normal(20, 5), 20) * base_travel_demand' * adjust_factor)
    elseif t >= 2 && t < 14
        push!(travel_demand, pdf(Normal(8, 5), t) / pdf(Normal(8, 5), 8)  * base_travel_demand * adjust_factor)
    else
        push!(travel_demand, pdf(Normal(20, 5), t) / pdf(Normal(20, 5), 20) * base_travel_demand' * adjust_factor)
    end
end
# FCEV retio, we think the FCEVs that don't require hydrogen as GVs
GV_ratio = 0.9998
FCEV_ration = 0.0002
# qʰALLTime shape: (T, OD)
qʰALLTime = [
    [ FCEV_ration * travel_demand[t][OD_pair.first, OD_pair.second] for OD_pair in OD]
        for t in 1:T
]

qᵍALLTime = [
    [GV_ratio * travel_demand[t][OD_pair.first, OD_pair.second]  for OD_pair in OD]
        for t in 1:T
]

"""
- @description: solve a TA problem
- @param {ρⱼ} hydrogen price for each HRS
- @param {qᵍ} GV traffic demand for each OD
- @param {qʰ} FCEV traffic demand for each OD
- @return {hydrogen_demand}  FCEV refuling demand for each HRS
"""
function solve_TA(ρⱼ, qᵍ, qʰ)
    # pick the cheapest hydrogen station for each FCEV path
    pathʰ2stationidx = [ [ reduce((x, y) -> (ρⱼ[x] < ρⱼ[y]) ? x : y, stations_idx)
        for stations_idx in paths] for paths in pathʰ_via_stationidx]

    # nonlinear optimization
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    register(model, :BRP_int, 3, BRP_int; autodiff = true)
    # JuMP.set_optimizer_attribute(model, "hessian_approximation", "limited-memory")

    ## decision variables
    @variable(model, xₐ[1:length(roads)] >= 0) # traffic flow for each road
    # traffic flow for each path
    fᵍ = [@variable(model, [1:length(pathᵍ[i])], lower_bound=0) for i in 1:length(OD)] # (OD, path)
    fʰ = [@variable(model, [1:length(pathʰ[i])], lower_bound=0) for i in 1:length(OD)] # (OD, path)

    ## constraints
    # road flow is sum of path flows through the road
    @constraint(model, [a=1:length(roads)],
        xₐ[a] == sum(fᵍ[idx.OD][idx.path] for idx in road_idx2pathᵍidx[a])
        + sum(fʰ[idx.OD][idx.path] for idx in road_idx2pathʰidx[a]))

    # OD demand is sum of their path flows
    @constraint(model, [i=1:length(OD)], qʰ[i] == sum(fʰ[i]))
    @constraint(model, [i=1:length(OD)], qᵍ[i] == sum(fᵍ[i]))

    ## objective
    @NLobjective(model, Min,
        ω * sum(BRP_int(xₐ[i], t₀[i], cₐ[i]) for i in 1:length(roads)) # cost of time
        + Eₕ * sum( ρⱼ[station_idx] * fʰ[OD_idx][path_idx] # cost of hydrogen
            for (OD_idx, paths) in enumerate(pathʰ2stationidx)
            for (path_idx, station_idx) in enumerate(paths) )
    )

    optimize!(model)
    println("【Traffic Result】:", raw_status(model))
    # if primal_status(model) == MOI.FEASIBLE_POINT
    #     println("Solution is a FEASIBLE POINT")
    # elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
    #     println("Solution is suboptimal due to a time limit, but a primal solution is available")
    #     println(ρⱼ)
    # else
    #     println("The model was not solved correctly.")
    #     println(ρⱼ)
    # end

    ## calculate the hydrogen demand for each HRS
    hydrogen_demand = zeros(Nₕ)
    for (OD_idx, paths) in enumerate(pathʰ2stationidx)
        for (path_idx, station_idx) in enumerate(paths)
            hydrogen_demand[station_idx] += Eₕ * value(fʰ[OD_idx][path_idx])
        end
    end

    return hydrogen_demand
end

## test
# solve_TA(ones(Nₕ), qᵍALLTime[1], qʰALLTime[1])

end


text:
# #######################################################################################################################################################
# cmaps: #Accent, Accent_r, Blues, Blues_r, BrBG, BrBG_r, BuGn, BuGn_r, BuPu, BuPu_r, CMRmap, CMRmap_r, Dark2, Dark2_r, GnBu, GnBu_r, Greens, Greens_r,
# Greys, Greys_r, OrRd, OrRd_r, Oranges, Oranges_r, PRGn, PRGn_r, Paired, Paired_r, Pastel1, Pastel1_r, Pastel2, Pastel2_r,
# PiYG, PiYG_r, PuBu, PuBuGn, PuBuGn_r, PuBu_r, PuOr, PuOr_r, PuRd, PuRd_r, Purples, Purples_r, RdBu, RdBu_r, RdGy, RdGy_r,
# RdPu, RdPu_r, RdYlBu, RdYlBu_r, RdYlGn, RdYlGn_r, Reds, Reds_r, Set1, Set1_r, Set2, Set2_r, Set3, Set3_r,
# Spectral, Spectral_r, Wistia, Wistia_r, YlGn, YlGnBu, YlGnBu_r, YlGn_r, YlOrBr, YlOrBr_r, YlOrRd, YlOrRd_r, afmhot, afmhot_r,
# autumn, autumn_r, binary, binary_r, bone, bone_r, brg, brg_r, bwr, bwr_r, cividis, cividis_r, cool, cool_r, coolwarm, coolwarm_r, copper, copper_r,
# cubehelix, cubehelix_r, flag, flag_r, gist_earth, gist_earth_r, gist_gray, gist_gray_r, gist_heat, gist_heat_r, gist_ncar, gist_ncar_r,
# gist_rainbow, gist_rainbow_r, gist_stern, gist_stern_r, gist_yarg, gist_yarg_r, gnuplot, gnuplot2, gnuplot2_r, gnuplot_r, gray, gray_r, hot, hot_r, hsv, hsv_r,
# inferno, inferno_r, jet, jet_r, magma, magma_r, nipy_spectral, nipy_spectral_r, ocean, ocean_r, pink, pink_r, plasma, plasma_r, prism, prism_r,
# rainbow, rainbow_r, seismic, seismic_r, spring, spring_r, summer, summer_r, tab10, tab10_r, tab20, tab20_r, tab20b, tab20b_r, tab20c, tab20c_r,
# terrain, terrain_r, twilight, twilight_r, twilight_shifted, twilight_shifted_r, viridis, viridis_r, winter, winter_r'
# #######################################################################################################################################################

using Dates
using PyPlot, Colors

include("startup.jl")
include("manualMassLibrary.jl")

# *****************************    DEFINE PLOT TIME    *********************************************************************
selectedTime = true
toi = [Dates.DateTime(2018,5,10), Dates.DateTime(2018,5,22)]	# have to restart julia for new toi #
interpolateTraces = true
avMinutes = 30
intpolT = collect(Dates.datetime2unix(toi[1]):(60*avMinutes):Dates.datetime2unix(toi[2])) 	# 5min averages
intpolDT = Dates.unix2datetime.(intpolT)

# ********************************    LOAD VARS    *************************************************************************
plotSMPS = true
plotNO3cims = true
plotACSM = true
plotMeteo = true
plotSRR = true
plotBC = false
#plotBLinfluence = false

if !(@isdefined createTotalPTR3data)
	println("load variables:")
	include("Bolivia-DefineVariables.jl")
else println("variables already defined!")
end

# ********************************    PREPARE FIGURE    ********************************************************************
figure(figsize=(20,20))
subplots_adjust(left=0.03, right=0.95, bottom = 0.03, top = 0.97, hspace=0.01)
N = length(findall([plotSMPS, plotNO3cims, plotACSM, plotMeteo, plotSRR])) + 1 		# nr of Plots

# ********************************    LOAD PTR3 DATA    ********************************************************************

if !(@isdefined ptr3Traces)
	println("  -> load PTR3 Data")
	ptr3Peaks, ptr3Times, ptr3Traces = createTotalPTR3data()
	ptr3Masses = round.(ptr3Peaks[:,1]; digits = 3)
	ptr3PeaksNoO, ptr3TimesNoO, ptr3TracesNoO = createTotalPTR3dataNonOxidized()
	ptr3MassesNoO = round.(ptr3PeaksNoO[:,1]; digits = 3)
end

massesToPlot = round.([
		H3O[1]
		# DMS[1]
		#MasslistFunctions.createMassList2(; C=8, O=2, H=20, N=0, S=0, nHplus=1)[1]	# butanol
		#MasslistFunctions.createMassList2(; C=8, O=2, H=18, N=0, S=0, nHplus=1)[1]	# butanol
		#MasslistFunctions.createMassList2(; C=4, O=1, H=10, N=0, S=0, nHplus=1)[1]	# butanol
		#MasslistFunctions.createMassList2(; C=4, O=0, H=8, N=0, S=0, nHplus=1)[1]	# butanol
		#MasslistFunctions.createMassList2(; C=4, O=1:2, H=6:2:8, N=0, S=0, nHplus=1)[1]
		# MasslistFunctions.createMassList2(; C=3, O=0:20, H=8:2:12, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=4, O=1, H=6, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=4, O=1, H=8, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=4, O=2, H=8, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=4, O=3:1:5, H=8, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=4, O=0, H=8, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=4, O=1, H=10, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=8, O=2, H=20, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=8, O=1, H=18, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=6, O=3, H=5, N=1, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=1, O=0:4, H=4:2:6, N=0, S=1, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=2, O=0:1, H=4:2:6, N=0, S=1, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=2, O=3:5, H=4:2:6, N=0, S=1, nHplus=1)[1]
	        #MasslistFunctions.createMassList2(; C=6, O=1:5, H=5, N=1, S=0, nHplus=1)[1]
		# MasslistFunctions.createMassList2(; C=6, O=2:3, H=5:1:7, N=0:1, S=0, nHplus=1)[1]
		# MasslistFunctions.createMassList2(; C=5, O=0:20, H=6:2:8, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=5, O=0:20, H=10:12, N=0, S=0, nHplus=1)[1]
		#MasslistFunctions.createMassList2(; C=4, O=1, H=8, N=0, S=0, nHplus=1)[1]
		# MasslistFunctions.createMassList2(; C=4, O=1, H=6, N=0, S=0, nHplus=1)[1]
		# always exclude C10H12O, C10H13O, C9H11ON, C6H12O4, C10H12O2, when looking at 12.5. - 15.5. !!!
		# always exclude C4H8 and C8H20O2 - Butanol !!!
		#MasslistFunctions.createMassList2(; C=1, O=1:5, H=1:10, N=0, S=0, nHplus=1, allowRadicals=true)[1]
		#MasslistFunctions.createMassList2(; C=9, O=1, H=10:8:18, N=0, S=0, nHplus=1, allowRadicals=true)[1]
		MasslistFunctions.createMassList2(; C=7, O=0:10, H=8, N=0, S=0, nHplus=1)[1]
		]; digits = 3)
massesToPlotnoO = round.([
		H3O[1]
		MasslistFunctions.createMassList2(; C=7, O=0, H=2:20, N=0:1, S=0:1, nHplus=1)[1]
		]; digits = 3)
wantedIndices = findall(x -> x in massesToPlot, ptr3Masses)
println("  -> interpolate chosen PTR3 data")
# !!!!!!!!! #
ptr3Traces2plot = InterpolationFunctions.interpolateMatrix_1D(intpolT, ptr3Times, ptr3Traces[:,wantedIndices]) ./ collect(range(0.101,stop = 0.70, step = 0.6/577))
createTotalPTR3dataNonOxidized()

# ********************************    PLOTTING    *************************************************************************

# PTR3 plot ##############################
n = 1
println("Plot Nr. $(n)")
PTR3plot = PyPlot.subplot(N,1,n, label="PTR3")
PTR3plot.patch.set_alpha(0)
PTR3plot.semilogy(intpolDT, ptr3Traces2plot) # !!!!!!!!!!!!!!!!!!!
#PTR3plot.semilogy(intpolDT, DatasetFunctions.nansum(ptr3Traces2plot, 2))
labels = ptr3Peaks[wantedIndices,2]
#PTR3plot.semilogy(intpolDT, ptr3Traces2plot[:,1]./ptr3Traces2plot[:,2])
#PTR3plot.semilogy(intpolDT, ptr3Traces2plot[:,3]./ptr3Traces2plot[:,4])
#labels = [string(ptr3Peaks[wantedIndices,2][1], "/", ptr3Peaks[wantedIndices,2][2]), string(ptr3Peaks[wantedIndices,2][3], "/", ptr3Peaks[wantedIndices,2][4])]
#append!(labels,"sum")
PTR3plot.legend(labels)
#PTR3plot.set_ylim(1e-2,2*DatasetFunctions.nanmax(DatasetFunctions.nansum(ptr3Traces2plot, 2)))
PTR3plot.set_xlim(toi[1], toi[2])
PTR3plot.set_yscale("linear")

n = n+1

# SRR plot ###############################
if plotSRR
	println("Plot Nr. $(n)")
	SRRplot = PyPlot.subplot(N,1,n, sharex = PTR3plot, label="SRR")
	cmap = get_cmap("tab20")
	colors = cmap(0:6)
	SRRplot.patch.set_alpha(0)
	SRRplot.stackplot(intpolDT, transpose(SRRs_06)[:,:], colors = colors)
	SRRplot.legend(SRR06names)
end

if (@isdefined SRRplot)
	n = n+1
end

# NO3cims Plot ###########################
if plotNO3cims
	println("Plot Nr. $(n)")
	NO3cimsPlot = PyPlot.subplot(N,1,n, sharex = PTR3plot, label="NO3cims")
	NO3cimsPlot.patch.set_alpha(0)
	NO3cimsPlot.plot(intpolDT, Nitrophenol,  color = "k" , linewidth = 2, label = "Nitrophenol")
	NO3cimsPlot.plot(intpolDT, MSAtrace, label = "MSA")
	NO3cimsPlot.plot(intpolDT, SO5trace, label = "SO5-")
	NO3cimsPlot.plot(intpolDT, HOMs, label = "HOMs")
	NO3cimsPlot.plot(intpolDT, SAtrace./DatasetFunctions.nanmean(SAtrace), label = "SA [ppq]")
	NO3cimsPlot.set_yscale("linear")
	NO3cimsPlot.set_ylim(1e-2, 1e2)
	NO3cimsPlot.legend()
end

if (@isdefined NO3cimsPlot)
	n = n+1
end

# SMPSplot ###############################
if plotSMPS
	println("Plot Nr. $(n)")
	SMPSplot = PyPlot.subplot(N,1,n, sharex = PTR3plot, label="SMPS")
	SMPSplot.patch.set_alpha(0)
	SMPSplot.contourf(intpolDT, log10.(SMPS_binsize.*1e9), transpose(SMPS_binconc), levels = 3 .^vcat(collect(0:0.4:8), collect(8.2:0.2:10.0), collect(10.1:0.1:10.8), collect(10.85:0.05:11.2)), cmap = "nipy_spectral")
end

if (@isdefined SMPSplot)
	n = n+1
end

# meteo plot #############################
if plotMeteo
	println("Plot Nr. $(n)")
	MeteoPlot = PyPlot.subplot(N,1,n, sharex = PTR3plot)
	MeteoPlot.patch.set_alpha(0)
	MeteoPlot.plot(intpolDT, SWd, color = "gold" , label = "sw rad [W/m²].")
	MeteoPlot.scatter(intpolDT, WD, s = 1, label = "WD [deg]")
	MeteoPlot.plot(intpolDT, WDstdv, label = "WDstdv")
	MeteoPlot.plot(intpolDT, TEMPstation, label = "Temp")
	MeteoPlot.plot(intpolDT, RHstation, label = "RH")
	MeteoPlot.plot(intpolDT, PRESSUREstation, label = "Pressure")

	#MeteoPlot.plot(intpolDT, BC_sum, label = "sum(BC)")
	#MeteoPlot.semilogy(intpolDT, RH, color = "slateblue", label = "RH [%]")
	if plotACSM == true && length(ACSM_org) > 0
		ACSMsum = (ACSM_org.+ACSM_sulf.+ACSM_nitr.+ACSM_ammon.+ACSM_chlor)
		MeteoPlot.plot(intpolDT, ACSMsum, color = "purple", label = "total ACSM signal")
	end
	MeteoPlot.legend()
	MeteoPlot.set_ylim(-5,370)
end

if (@isdefined MeteoPlot)
	n = n+1
end

# ACSM plot ##############################
if plotACSM == true && length(ACSM_org) > 0
	println("Plot Nr. $(n)")
	ACSMsum = (ACSM_org.+ACSM_sulf.+ACSM_nitr.+ACSM_ammon.+ACSM_chlor)
	ACSMplot = PyPlot.subplot(N,1,n, sharex = PTR3plot)
	ACSMplot.patch.set_alpha(0)
	ACSMplot.stackplot(intpolDT, ACSM_org./ACSMsum, ACSM_sulf./ACSMsum, ACSM_nitr./ACSMsum, ACSM_ammon./ACSMsum, ACSM_chlor./ACSMsum)
	ACSMlabels = ["org.", "sulf.", "nitr.", "ammon.", "chlor."]
	ACSMplot.legend(ACSMlabels)
	ACSMplot.set_ylabel("composition")
	ACSMplot.set_ylim([0,1])
	ACSMplot2 = ACSMplot.twinx()
	ACSMplot2.set_ylabel("ACSM sum [ug/m³]")
	ACSMplot2.semilogy(intpolDT, ACSMsum, color="black",linewidth = 2)
else println("no ACSM data available")
end

#winddir plot?
#WS		# wind strength
#WD		# wind direction

#write daily traces to file:
#global X = DatasetFunctions.nanmean(ptr3Traces2plot[(hour.(intpolDT).== 0),:], 1)
#for i in collect(1:23)
#	x = DatasetFunctions.nanmean(ptr3Traces2plot[(hour.(intpolDT).== i),:], 1)
#	global X = vcat(X, x)
#end
#
#open("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/DailyTrace/ptr3Dat.csv", "w") do io
#           writedlm(io, vcat(hcat(["hour"],permutedims(labels)), hcat(collect(0:23), X)))
#       end


text:
# Random Walk Metropolis for Gibbs.
"""
`name::Symbol`

`stepper::Function`  function that takes `(model::Model, state::T) where T` and returns updated value for parameter with name `name`.

`proposal::Union{Normal,MvNormal}`  proposal distribution
"""
struct RWM{S<:OneOrMoreSymbols, F<:Function, T<:Metropolis}
  name::S
  rwm::T
  stepper::F
end

"""
`bijector`: This should be the transformtion to the real space. e.g. if X has positive support, the Log transform should be supplied.
"""
function RWM(name::Symbol, logprob::Function, met::Metropolis; bijector=nothing)
  if bijector === nothing
    if met isa StaticRWM
      if met.proposal isa Normal
        bijector = Bijectors.Identity{0}()
      else
        bijector = Bijectors.Identity{1}()
      end
    elseif met isa MvAdaptiveRWM
      bijector = Bijectors.Identity{1}()
    else
      error("Not implemented!")
    end
  end
  invb = inv(bijector)

  function stepper(model::Model, state::S) where S
    function _logprob(real_x)
      x, labsdj = forward(invb, real_x)
      return logprob(model, state, x) + labsdj
    end
    return invb(update(met, bijector(state[name]), _logprob))
  end

  return RWM(name, met, stepper)
end

function RWM(name::Symbol, logprob::Function, proposal::Union{Normal, MvNormal};
             bijector=nothing)
  srwm = StaticRWM(proposal)
  return RWM(name, logprob, srwm, bijector=bijector)
end

# TODO: Implement this.
# function RWM(name::Tuple, logprob::Function, proposal::Union{Normal, MvNormal};
#              bijector=nothing)
# end