text:
#-------------------------------------------------------------------------------
# File: vns_lprc.jl
# Description: This files contains all functions that are used in VNS_LPRC.
# Date: October 31, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

"""
    perturbation_VNS_LPRC_insertion!(solution::Solution, k::Int,
                                     critical_cars::Array{Int, 1},
                                     instance::Instance)

Deletes `k` random critical cars from the sequence. Inserts back those cars
according to a greedy criterion.
"""
function perturbation_VNS_LPRC_insertion!(solution::Solution, k::Int,
                                          critical_cars::Array{Int, 1},
                                          instance::Instance)
    remove!(solution, instance, k, critical_cars)
    for i in 1:k
        greedy_add!(solution, instance, instance.nb_cars, 2)
    end
    return solution
end

"""
    perturbation_VNS_LPRC_exchange!(solution::Solution, k::Int, instance::Instance)

Randomly applies k exchange between cars that are involved in the same
constraints.
"""
function perturbation_VNS_LPRC_exchange!(solution::Solution, k::Int, instance::Instance)
    # Array that contain only key defining more than one car
    valid_key_HPRC = Array{Int,1}(undef, 0)
    # for key in keys(instance.same_HPRC)
    #   if length(instance.same_HPRC[key]) >= 2
    #       push!(valid_key_HPRC, key)
    # OR:
    for group in instance.same_HPRC
        if length(group.second) >= 2
            push!(valid_key_HPRC, group.first)
        end
    end

    for _ in 1:k
        same_HPRC_array = instance.same_HPRC[rand(valid_key_HPRC)]
        index_car_a = rand(same_HPRC_array)
        index_car_b = rand(same_HPRC_array)
        # Cannot be the same
        while index_car_a == index_car_b
            index_car_b = rand(same_HPRC_array)
        end
        move_exchange!(solution, index_car_a, index_car_b, instance)
    end

    return solution
end

"""
    perturbation_VNS_LPRC(solution_init::Solution, p::Int, k::Int, instance::Instance)

Calls both perturbations and return a new solution. This function does not
modify the solution given in parameters.
"""
function perturbation_VNS_LPRC(solution_init::Solution, p::Int, k::Int, instance::Instance)
    solution = deepcopy(solution_init)
    if p == 1
        perturbation_VNS_LPRC_exchange!(solution, k, instance)
    else
        critical_cars = find_critical_cars(solution, instance, 2)
        k = min(k, length(critical_cars))
        perturbation_VNS_LPRC_insertion!(solution, k, critical_cars, instance)
    end
    return solution
end

"""
    local_search_VNS_LPRC!(solution::Solution, perturbation_exchange::Bool, instance::Instance)

Optimizes the weighted sum of first and second objectives. When
`perturbation_exchange` is `true` this function does not degrade the first
objective.
"""
function local_search_VNS_LPRC!(solution::Solution, perturbation_exchange::Bool, instance::Instance, start_time::UInt)
    # useful variable
    b0 = instance.nb_late_prec_day+1

    improved = true
    while improved && TIME_PART_VNS_LPRC * TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        improved = false
        critical_cars = find_critical_cars(solution, instance, 2)
        for index_car_a in critical_cars
            best_delta = 0
            best_positions = Array{Int, 1}()
            hprc_value = instance.HPRC_keys[solution.sequence[index_car_a]]
            for index_car_b in b0:instance.nb_cars
                if !perturbation_exchange || (index_car_a != index_car_b && index_car_b in instance.same_HPRC[hprc_value])
                    delta = weighted_sum(cost_move_exchange(solution, index_car_a, index_car_b, instance, 2))
                    if delta < best_delta
                        best_positions = Array{Int, 1}([index_car_b])
                        best_delta = delta
                    elseif delta == best_delta
                        push!(best_positions, index_car_b)
                    end
                end
            end
            if !isempty(best_positions)
                index_car_b = rand(best_positions)
                move_exchange!(solution, index_car_a, index_car_b, instance)
                if best_delta < 0
                    improved = true
                end
            end
        end
    end

    return solution
end

"""
    local_search_intensification_VNS_LPRC_exchange!(solution::Solution, instance::Instance)

Optimizes the weighted sum of first and second objectives using `move_exchange!`.
"""
function local_search_intensification_VNS_LPRC_exchange!(solution::Solution, instance::Instance, start_time::UInt)
    # useful variable
    b0 = instance.nb_late_prec_day+1

    improved = true
    while improved && TIME_PART_VNS_LPRC * TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        improved = false
        critical_cars = find_critical_cars(solution, instance, 2)
        for index_car_a in critical_cars
            best_delta = 0
            best_positions = Array{Int, 1}()
            for index_car_b in b0:instance.nb_cars
                if index_car_a != index_car_b
                    delta = weighted_sum(cost_move_exchange(solution, index_car_a, index_car_b, instance, 2))
                    if delta < best_delta
                        best_positions = Array{Int, 1}([index_car_b])
                        best_delta = delta
                    elseif delta == best_delta
                        push!(best_positions, index_car_b)
                    end
                end
            end
            if !isempty(best_positions)
                index_car_b = rand(best_positions)
                move_exchange!(solution, index_car_a, index_car_b, instance)
                if best_delta < 0
                    improved = true
                end
            end
        end
    end

    return solution
end

"""
    local_search_intensification_VNS_LPRC_insertion!(solution::Solution, instance::Instance)

Optimizes the weighted sum of first and second objectives using `move_insertion!`.
"""
function local_search_intensification_VNS_LPRC_insertion!(solution::Solution, instance::Instance, start_time::UInt)
    # useful variable
    b0 = instance.nb_late_prec_day+1

    improved = true
    while improved && TIME_PART_VNS_LPRC * TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        improved = false
        critical_cars = find_critical_cars(solution, instance, 2)
        for index_car in critical_cars
            best_delta = 0
            best_positions = Array{Int, 1}()
            matrix_deltas = cost_move_insertion(solution, index_car, instance, 2)
            for position in b0:instance.nb_cars
                delta = weighted_sum(matrix_deltas[position, :])
                if delta < best_delta
                    best_positions = Array{Int, 1}([position])
                    best_delta = delta
                elseif delta == best_delta
                    push!(best_positions, position)
                end
            end

            if !isempty(best_positions)
                index_insert = rand(best_positions)
                move_insertion!(solution, index_car, index_insert, instance)
                if best_delta < 0
                    improved = true
                end
            end
        end
    end

    return solution
end

"""
    intensification_VNS_LPRC!(solution::Solution, instance::Instance)

Calls both intensification.
"""
function intensification_VNS_LPRC!(solution::Solution, instance::Instance, start_time::UInt)
    local_search_intensification_VNS_LPRC_insertion!(solution, instance, start_time)
    local_search_intensification_VNS_LPRC_exchange!(solution, instance, start_time)
    return solution
end

"""
    is_better_VNS_LPRC(solution_1::Solution, solution_2::Solution, instance::Instance)

Returns `true` if the weighted sum of `solution_1`'s objective value is better
than the weighted sum of `solution_2`'s objective value and if the HPRC value of
`solution_1` is better or equal than the HPRC value of `solution_2`.
Returns `false` otherwise.
"""
function is_better_VNS_LPRC(solution_1::Solution, solution_2::Solution, instance::Instance)
    solution_1_cost = cost(solution_1, instance, 2)
    solution_2_cost = cost(solution_2, instance, 2)

    cost_better = weighted_sum(solution_1_cost) <= weighted_sum(solution_2_cost)
    HPRC_not_worse = solution_1_cost[1] <= solution_2_cost[1]

    return cost_better && HPRC_not_worse
end

"""
    is_strictly_better_VNS_LPRC(solution_1::Solution, solution_2::Solution, instance::Instance)

Returns `true` if the weighted sum of `solution_1`'s objective value is better
than the weighted sum of `solution_2`'s objective value and if the HPRC value of
`solution_1` is better or equal than the HPRC value of `solution_2`.
Returns `false` otherwise.
"""
function is_strictly_better_VNS_LPRC(solution_1::Solution, solution_2::Solution, instance::Instance)
    solution_1_cost = cost(solution_1, instance, 2)
    solution_2_cost = cost(solution_2, instance, 2)

    cost_better = weighted_sum(solution_1_cost) < weighted_sum(solution_2_cost)
    HPRC_not_worse = solution_1_cost[1] <= solution_2_cost[1]

    return cost_better && HPRC_not_worse
end

"""
    VNS_LPRC(solution_init::Solution, instance::Instance, start_time::UInt)

Optimizes the weighted sum of first and second objectives.
"""
function VNS_LPRC(solution_init::Solution, instance::Instance, start_time::UInt)
    # p = 0 implies insertion move and p = 1 implies exchange move as stated
    # in section 6.1. Note that section 6.5 states the opposite.
    _bitarray = BitArray{1}([false, true, false])

    # variables of the algorithm
    solution_best = deepcopy(solution_init)
    k_min = (VNS_LPRC_MIN_INSERT, VNS_LPRC_MIN_EXCHANGE)
    k_max = (VNS_LPRC_MAX_INSERT, VNS_LPRC_MAX_EXCHANGE)
    p = 0
    k = k_min[p+1]
    nb_intens_not_better = 0
    vector_zero = zeros(Int, 3)
    while (nb_intens_not_better < VNS_LPRC_MAX_NON_IMPROVEMENT
          && TIME_PART_VNS_LPRC * TIME_LIMIT > (time_ns() - start_time) / 1.0e9
          && cost(solution_best, instance, _bitarray) != vector_zero)
        while (k <= k_max[p+1]
              && TIME_PART_VNS_LPRC * TIME_LIMIT > (time_ns() - start_time) / 1.0e9)
            neighbor = perturbation_VNS_LPRC(solution_best, p, k, instance)
            local_search_VNS_LPRC!(neighbor, p == 1, instance, start_time)
            if is_strictly_better_VNS_LPRC(neighbor, solution_best, instance)
                k = k_min[p+1]
                nb_intens_not_better = 0
            else
                k = k + 1
            end
            if is_better_VNS_LPRC(neighbor, solution_best, instance)
                solution_best = deepcopy(neighbor)
            end
        end
        intensification_VNS_LPRC!(solution_best, instance, start_time)
        nb_intens_not_better += 1

        p = 1 - p
        k = k_min[p+1]
    end

    return solution_best
end


text:
using Plots

function plot_results()
    x = map(x->10^x,1:6)
    elixir=[3.355,3.571,7.871,50.252,506.419,4764.373]
    julia=[0.025,0.190,1.954,22.580,254.876,5176]
    plot(x,elixir,xlabel="processes",ylabel="time [ms]",xscale=:log10,yscale=:log10,lab="elixir",lw=2)
    plot!(x,julia,lab="julia",lw=2)
    title!("Lightweight processes, Julia vs. Elixir", legend=:bottomright)
end

ext: jl


text:
using Test
using OrdinaryDiffEq
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

@testset "Test the different diffusion models" begin
    prob = prob_ode_fitzhughnagumo
    true_sol = solve(prob, Tsit5(), abstol=1e-12, reltol=1e-12)

    @testset "Time-Varying Diffusion" begin
        sol = solve(prob, EK0(diffusionmodel=:dynamic), adaptive=false, dt=1e-4)
        @test sol.u ≈ true_sol.(sol.t)
    end

    @testset "Time-Fixed Diffusion" begin
        sol = solve(prob, EK0(diffusionmodel=:fixed), adaptive=false, dt=1e-4)
        @test sol.u ≈ true_sol.(sol.t)
    end

    @testset "Time-Varying Diagonal Diffusion" begin
        sol = solve(prob, EK0(diffusionmodel=:dynamicMV), adaptive=false, dt=1e-4)
        @test sol.u ≈ true_sol.(sol.t)
    end

    @testset "Time-Fixed Diagonal Diffusion" begin
        sol = solve(prob, EK0(diffusionmodel=:fixedMV), adaptive=false, dt=1e-4)
        @test sol.u ≈ true_sol.(sol.t)
    end

    @testset "Time-Fixed Diffusion MAP" begin
        sol = solve(prob, EK0(diffusionmodel=:fixedMAP), adaptive=false, dt=1e-4)
        @test sol.u ≈ true_sol.(sol.t)
    end

end


text:
import Grid
using Base.Test

ctarget = ([1,1],[2,1],[1,2],[2,2],[1,3],[2,3])
k = 0
for c in Grid.Counter((2,3))
    k += 1
    @test c == ctarget[k]
end
@test k == 6
k = 0
for c in Grid.Counter((3,0))
    k += 1
    @test c == ctarget[k]
end
@test k == 0
k = 0
for c in Grid.Counter(())
    k += 1
    @test c == ctarget[k]
end
@test k == 0


text:
function solve_step!(
  x1::NTuple{3,AbstractVector},
  solver::Newmark,
  op::ConstantMatrixODEOperator,
  x0::NTuple{3,AbstractVector},
  t0::Real,
  cache) # -> (uF,tF)

  dt = solver.dt
  γ = solver.γ
  β = solver.β
  t1 = t0+dt
  u0, v0, a0 = x0
  u1, v1, a1 = x1

  if cache === nothing
    newmatrix = true

    # Allocate caches
    newmark_cache = allocate_cache(op,v0,a0)
    (v,a, ode_cache) = newmark_cache

    # Define Newmark operator
    newmark_affOp = NewmarkConstantMatrixOperator(op,t1,dt,γ,β,(u0,v0,a0),newmark_cache)

    # Allocate matrices and vectors
    A, b = _allocate_matrix_and_vector(op,x0,ode_cache)
    jacobian!(A,newmark_affOp,u1)

    # Create affine operator cache
    affOp_cache = (A,b,newmark_affOp,nothing)
  else
    newmatrix = false
    newmark_cache, affOp_cache = cache
  end

  # Unpack and update caches
  (v,a, ode_cache) = newmark_cache
  ode_cache = update_cache!(ode_cache,op,t1)
  A,b,newmark_affOp,l_cache = affOp_cache

  # Fill vector
  newmark_affOp = NewmarkConstantMatrixOperator(op,t1,dt,γ,β,x0,newmark_cache)
  residual!(b,newmark_affOp,u1)

  # Create affine operator with updated RHS
  affOp = AffineOperator(A,b)
  l_cache = solve!(u1,solver.nls,affOp,l_cache,newmatrix)

  # Update auxiliar variables
  u1 = u1 + u0
  v1 = γ/(β*dt)*(u1-u0) + (1-γ/β)*v0 + dt*(1-γ/(2*β))*a0
  a1 = 1.0/(β*dt^2)*(u1-u0) - 1.0/(β*dt)*v0 - (1-2*β)/(2*β)*a0

  # Pack caches
  affOp_cache = A,b,newmark_affOp,l_cache
  cache = (newmark_cache, affOp_cache)
  x1 = (u1,v1,a1)

  return (x1,t1,cache)

end

"""
Affine operator that represents the Newmark Affine operator with constant
matrix at a given time step, i.e., M(u_n+1-u_n)/dt + K u_n+1 + b(t)
"""
mutable struct NewmarkConstantMatrixOperator <: NonlinearOperator
  odeop::ConstantMatrixODEOperator
  t1::Float64
  dt::Float64
  γ::Float64
  β::Float64
  x0::NTuple{3,AbstractVector}
  ode_cache
end

function residual!(b::AbstractVector,op::NewmarkConstantMatrixOperator,x::AbstractVector)
  u1 = x
  u0, v0, a0 = op.x0
  v1, a1, cache = op.ode_cache
  a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
  v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
  residual!(b,op.odeop,op.t1,(u1,v1,a1),cache)
  b .*= -1.0
end

function jacobian!(A::AbstractMatrix,op::NewmarkConstantMatrixOperator,x::AbstractVector)
  u1 = x
  u0, v0, a0 = op.x0
  v1, a1, cache = op.ode_cache
  a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
  v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
  z = zero(eltype(A))
  fillstored!(A,z)
  jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
end

function _allocate_matrix(odeop::NewmarkConstantMatrixOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
  A = allocate_jacobian(odeop,x[1],ode_cache)
  return A
end

function _allocate_matrix_and_vector(odeop::NewmarkConstantMatrixOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
  b = allocate_residual(odeop,x[1],ode_cache)
  A = allocate_jacobian(odeop,x[1],ode_cache)
  return A, b
end


text:
using Flux
import KitML

# dataset
cd(@__DIR__)
dataset = KitML.load_data("data_mn.csv"; mode = :csv, dlm = ",")

baseDimension = 4 #length of u

for name in names(dataset)
    println(":", name, ", ")
end

L = size(dataset, 1)

# surrogate 1: moments -> entropy
begin
    X = zeros(baseDimension, size(dataset, 1))
    Y = zeros(1, size(dataset, 1))
    for j in axes(X, 2)
        X[1, j] = dataset.u_0[j]
        X[2, j] = dataset.u_1[j]
        X[3, j] = dataset.u_2[j]
        X[4, j] = dataset.u_3[j]
        Y[1, j] = dataset.h[j]
    end
    data = Flux.Data.DataLoader(X, Y, batchsize = 256, shuffle = true)

    #Base.show(nn)
    #nn = Chain(
    #    Dense(4, 16, relu),
    #    Dense(16, 16, relu),
    #    Dense(16, 16, relu),
    #    Dense(16, 1),
    #)
    nn = KitML.create_neural_closure(baseDimension, 1, 8, 3)

end
KitML.sci_train!(nn, data, ADAM(), 1000)
KitML.save_model(nn)
println("End ICNN Training")

# surrogate 2: ResNet
nn2 = Chain(
    Dense(4, 16, relu),
    KitML.Shortcut(Chain(Dense(16, 16, relu), Dense(16, 16, relu))),
    Dense(16, 1),
)
KitML.sci_train!(nn2, data, ADAM(), 2)

"""Notice the training of mapping u to Lagrangian multiplier λ seems not working well"""


text:
using NoiseJul
using Test

@testset "NoiseJul.jl" begin
    # Write your own tests here.
    @test isempty(NoiseJul.IO.SeisData())
    @test isempty(NoiseJul.Noise.FFTData())
end


text:
using BitFloats, Test

using BitFloats: BigFloat_mpfr, BuiltinInts, Int80, UInt80, decompose, explicit_bit,
                 exponent_half, exponent_mask, exponent_one, sign_mask, significand_mask,
                 uinttype, uniontypes

@testset "definitions" begin
    @test @isdefined Float80
    @test @isdefined UInt80
    @test @isdefined Float128
    @test sizeof(Float80)  == 10
    @test sizeof(UInt80)   == 10
    @test sizeof(Float128) == 16
    @test Float80  <: AbstractFloat
    @test UInt80   <: Unsigned
    @test Float128 <: AbstractFloat
    @test isprimitivetype(Float80)
    @test isprimitivetype(UInt80)
    @test isprimitivetype(Float128)
    @test uinttype(Float80)  == UInt80
    @test uinttype(Float128) == UInt128
end

# generate a random float in the whole range
function _rand(T::Type)
    r = rand(uinttype(T))
    s = r & significand_mask(T)
    e = r & exponent_mask(T)
    if T === Float80
        s >>= 1 # delete explicit bit
        if e != 0
            s |= explicit_bit(T)
        end
    end
    reinterpret(T, r & sign_mask(T) | e | s)
end

@testset "traits" begin
    @test explicit_bit() == explicit_bit(Float80) isa UInt80
    for T in (Float80, Float128)
        for ufun ∈ (sign_mask, exponent_mask, exponent_one, significand_mask)
            @test ufun(T) isa uinttype(T)
        end
        for ffun ∈ (eps, floatmin, floatmax, typemin, typemax)
            @test ffun(T) isa T
        end
        x = _rand(T)
        # @test floatmin(T) <= x <= floatmax(T)
        # @test nextfloat(T(1)) - T(1) == eps(T)
    end
    @test Inf80  isa Float80
    @test NaN80  isa Float80
    @test Inf128 isa Float128
    @test NaN128 isa Float128
    @test precision(Float80)  == 64
    @test precision(Float128) == 113
    @test Inf80  == typemax(Float80)  == -typemin(Float80)
    @test Inf128 == typemax(Float128) == -typemin(Float128)
    @test isnan(NaN80)
    @test isnan(NaN128)
    @test !isnan(_rand(Float80)) # very unlikely to fail
    @test !isnan(_rand(Float128))
    @test isinf(Inf80)
    @test isinf(Inf128)
    @test !isinf(_rand(Float80))
    @test !isinf(_rand(Float128))
end

@testset "exponent & significand" begin
    for F = (Float80, Float128)
        for k = -8.0:8.0
            n = 2.0^k
            x = rand(F) + 1
            @test exponent(n*x) == k
            @test significand(n*x) == significand(x)
        end
        @test significand(zero(F)) == zero(F)
        @test significand(F(NaN)) === F(NaN)
        @test significand(F(Inf)) === F(Inf)
        @test significand(F(-Inf)) === F(-Inf)
        @test_throws DomainError exponent(zero(F))
        @test_throws DomainError exponent(F(NaN))
        @test_throws DomainError exponent(F(Inf))
        @test_throws DomainError exponent(F(-Inf))
    end
end

@testset "nextfloat" begin
    for F = (Float80, Float128)
        while (x = _rand(F); isnan(x)) end
        setprecision(precision(x)) do
            b = BigFloat(x)
            @test x == b
            y = nextfloat(x)
            @test issubnormal(y) || issubnormal(x) || nextfloat(b) == y
            y = prevfloat(y)
            @test y == x
            y = prevfloat(y)
            @test issubnormal(y) || issubnormal(x) || prevfloat(b) == y
            @test nextfloat(y) == x
        end
    end
end

@testset "issubnormal" begin
    for F = (Float80, Float128)
        x = floatmin(F)
        @test !issubnormal(x)
        @test !issubnormal(-x)
        @test issubnormal(prevfloat(x))
        @test issubnormal(nextfloat(-x))
        @test !issubnormal(rand(F))
    end
end

@testset "ldexp & eps" begin
    for F = (Float80, Float128)
        x = _rand(F)
        b = big(x)
        d = 2^rand(1:16)
        n = rand(-d:d)
        r = ldexp(b, n)
        @test ldexp(x, n) == F(r) # conversion to F, in case of under or overflow

        # eps
        @test eps(F) == eps(F(1))
        @test isnan(eps(F(NaN)))
        @test isnan(eps(F(Inf)))
        @test isnan(eps(F(-Inf)))
    end
end

@testset "conversions" begin
    for F = (Float80, Float128)
        for T = (uniontypes(BuiltinInts)..., Float16, Float32, Float64, Float80, Float128)
            t = rand(T)
            @test F(t) isa F
            T == Bool && continue
            if T <: Integer
                @test unsafe_trunc(T, F(t)) isa T
                @test trunc(T, F(t)) isa T
                @test trunc(T, F(2.3)) == 2
                if T <: Signed
                    @test trunc(T, F(-2.3)) == -2
                else
                    @test_throws InexactError trunc(T, F(-2.3))
                end
                if F != Float128 # BROKEN
                    @test_throws InexactError T(F(2.3))
                end

                @test promote_type(F, T) == F
            else
                @test T(F(t)) isa T
                @test isnan(F(T(NaN)))
                @test isnan(T(F(NaN)))
                @test isinf(F(T(Inf)))
                @test isinf(T(F(Inf)))
                @test isinf(F(T(-Inf)))
                @test isinf(T(F(-Inf)))

                R = sizeof(T) < sizeof(F) ? F : T
                @test promote_type(F, T) == R
                @test F(zero(T)) == T(zero(F)) == 0
            end
            @test one(F) + one(T) === promote_type(F, T)(2)
        end
        @test one(F) isa F
        @test zero(F) isa F
        x = _rand(F)
        @test reinterpret(Unsigned, x) === reinterpret(uinttype(F), x)
        @test reinterpret(Signed, x) === reinterpret(F == Float80 ? Int80 : Int128, x)
        @test trunc(Signed, F(1.2)) === 1
        @test trunc(Integer, F(1.2)) === 1
        @test trunc(Unsigned, F(1.2)) === UInt(1)
        @test Signed(F(1)) === 1
        @test Unsigned(F(1)) === UInt(1)
        if F != Float128 #BROKEN
            @test_throws InexactError Signed(F(1.2))
            @test_throws InexactError Unsigned(F(1.2))
        end
        # BigFloat
        @test BigFloat(F(1)) == 1.0
        @test isnan(BigFloat(F(NaN)))
        @test BigFloat(F(Inf)) == Inf
        @test BigFloat(F(-Inf)) == -Inf
        x = rand(Int32)+rand()
        @test BigFloat(F(x)) == x
        @test F(BigFloat(F(x))) === F(x)
        y = _rand(F)
        @test isequal(F(BigFloat(y)), y)
    end
    x = _rand(Float128)
    @test isequal(BigFloat(x), BigFloat_mpfr(x))
end

@testset "round" begin
    for F = (Float80,) # broken for Float128
        @test round(F(1.2), RoundToZero) == F(1.0)
        @test round(F(1.2), RoundNearest) == F(1.0)
        @test round(F(1.2), RoundDown) == F(1.0)
        @test round(F(1.2), RoundUp) == F(2.0)
        @test round(F(1.8), RoundToZero) == F(1.0)
        @test round(F(1.8), RoundNearest) == F(2.0)
        @test round(F(1.8), RoundDown) == F(1.0)
        @test round(F(1.8), RoundUp) == F(2.0)
    end
end

@testset "comparisons" begin
    for F = (Float80, Float128)
        x, y = _rand(F), _rand(F)
        for op = (==, !=, <, <=, >, >=, isless, isequal)
            @test op(x, y) isa Bool
        end
        @test F(1) == F(1)
        @test F(1) != F(2)
        @test F(1) <  F(2)
        @test F(1) <= F(2)
        @test F(1) <= F(1)
        @test F(2) >  F(1)
        @test F(2) >= F(1)
        @test F(1) >= F(1)
        @test (x == y) == !(x != y)
        @test !(F(1) == F(2))
        @test !(F(1) != F(1))
        @test isequal(x, x)
        @test isequal(y, y)
        @test isequal(x, y) || isless(x, y) || isless(y, x)
        N = F(NaN)
        @test N != N
        @test isequal(N, N)
        @test !(N == N)
    end
end

@testset "arithmetic" begin
    for T = (Float80, Float128)
        n = rand(Int)
        for op = (*, /, +, -, rem, ^)
            for randfun = (_rand, rand)
                a, b = randfun(T), randfun(T)
                randfun == _rand && op == (^) && T == Float128 && continue # BUG crash otherwise
                r = op(a, b)
                @test r isa T
                op == (^) && T == Float128 && continue # BUG crash otherwise
                @test op(a, n) isa T
            end
        end
        x = _rand(T)
        @test -x  isa T
        @test -(-x) == x
        @test x >= 0 ? x == abs(x) : x == -abs(x)
        @test abs(T(-1)) == T(1)
        @test abs(T(1)) == T(1)
        @test T(2)^3 == 8
        @test T(2)^-1.0 == 0.5
        # LLVM intrinsics only called for Float80, otherwise conversion to BigFloat
        @test log2(T(16)) == 4
        @test log2(T(32)) == 5
        @test exp2(T(10)) == 1024
        @test sqrt(T(16)) == 4
        @test sqrt(T(25)) == 5
        # ≈ for sin/cos because this is not exact with Float128, due to intermediate conversion to BigFloat
        @test sin(T(big(pi)/2)) ≈ 1.0
        @test sin(T(-big(pi)/2)) ≈ -1.0
        @test cos(T(big(pi))) ≈ -1.0
        @test cos(T(0)) ≈ 1.0
        @test exp(T(0)) == 1
        @test exp(T(1)) ≈ exp(1)
        @test log(T(exp(1))) ≈ 1 atol=10^-16
        @test log10(T(10)) == 1.0
        T == Float128 && continue
        @test log(T(1)) == 0 # segfaults for Float128
    end
end

@testset "hashing" begin
    for T = (Float80, Float128)
        x = _rand(T)
        num, pow, den = BitFloats.decompose(x)
        # TODO: WAT? this intermediate value y needs to be computed for the test to not fail
        y = T(num) * T(2)^pow / den
        @test T(num) * T(2)^pow / den == x
    end
    n = rand(Int)
    @test hash(n) == hash(Float80(n)) == hash(Float128(n))
    f = _rand(Float64)
    @test hash(f) == hash(Float80(f)) == hash(Float128(f))
end

@testset "rand" begin
    for T = (Float80, Float128)
        x = rand(T)
        @test x isa T
        u = reinterpret(uinttype(T), x)
        @test u & sign_mask(T) == 0
        u2 = reinterpret(uinttype(T), (rand(T) + one(T)))
        u2 &= exponent_mask(T)
        if T == Float80
            u2 |= explicit_bit(T)
        end
        @test u2 == exponent_one(T)
    end
end

@testset "misc" begin
    for F = (Float80, Float128)
        x = _rand(F)
        @test bswap(x) !== x
        @test bswap(bswap(x)) === x

        @test "$(F(12))" == "12.0"
    end
end


text:
using BHAPtfem

ProjDir = dirname(@__FILE__)

#=
Compare formulas at:
http://www.awc.org/pdf/codes-standards/publications/design-aids/AWC-DA6-BeamFormulas-0710.pdf
=#

N = 20
Np1 = N + 1
Nhp1 = Int(N/2) + 1

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(N, Np1, 3, 1, 1, Line(2, 3)),
  :properties => [1.0e6 1.0e6 1.0e6 3.0e5;],
  :y_coords => zeros(Np1),
  :x_coords => zeros(Np1),
  :g_num => [
    collect(1:N)';
    collect(2:Np1)'],
  :support => [
    (1, [0 0 0 0 0 0]),
    #(Int(N/4),  [0 0 0 0 1 0]),
    (Int(3N/4),  [0 0 0 0 1 0]),
    (Np1, [1 0 0 0 1 0]),
    ],
  :loaded_nodes => [(i, [-20000.0/Np1 0.0 0.0 0.0 0.0 0.0]) for i in 1:Np1]
)

data[:z_coords] = VERSION.minor < 7 ? linspace(0, 4, Np1) :  range(0, stop=4, length=Np1)

@time m, dis_df, fm_df = PtFEM.p44(data)

data |> display
println()

m, dis_df, fm_df = PtFEM.p44(data)

println("Displacements:")
display(dis_df)
println()

println("Actions:")
display(fm_df)
println()

#=
using Plots
gr(size=(400,600))

p = Vector{Plots.Plot{Plots.GRBackend}}(3)
p[1] = plot(m.z_coords, m.displacements[1,:], ylim=(-0.0035, 0.003), lab="Displacement",
 xlabel="x [m]", ylabel="deflection [m]", color=:red)
p[2] = plot(m.actions[1,:], lab="Shear force", ylim=(-9000, 15000), xlabel="element",
  ylabel="shear force [N]", palette=:greens,fill=(0,:auto),α=0.6)
p[3] = plot(m.actions[11,:], lab="Moment", ylim=(-4200, 3000), xlabel="element",
  ylabel="moment [Nm]", palette=:grays,fill=(0,:auto),α=0.6)

plot(p..., layout=(3, 1))
savefig(ProjDir*"/figure-02.png")
=#
#=
plot!()
gui()
=#


text:
function Base.show(io::IO, data::Data)
    print(io, "CARP Data $(data.name)")
    print(io, " ($(length(data.vertices)) nodes,")
    print(io, " $(length(data.edges)) edges,")
    print(io, " $(length(data.requireds)) reqs)")
    print(io, " [$(data.lb), $(data.ub)]")
end

function Base.show(io::IO, vertex::Vertex)
    print(io, "v$(vertex.id)")
end

function Base.show(io::IO, edge::Edge)
    print(io, "e($(edge.from.id), $(edge.to.id))")
end


text:
using Base.Test
include("BasicCounting.jl")

## Start Test Linked List
# workspace()
include("LinkedList2.jl")
w = List()
@test isnull(w.head)

dat1 = NodeData(1, 3)
add!(w, dat1)
@test !isnull(w.head)

println(isnull(w.head.next))
# traverse(w)

# @test w.count == 1
#
# dat2 = ItemData(0, 2)
# addItem!(w, dat2)
# @test w.count == 2
## Start Test Linked List

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
"""
    fit_iht(y, x, z; k=10, d = Normal(), l=IdentityLink(),...)

Fits a model on design matrix (genotype data) `x`, response (phenotype) `y`,
and non-genetic covariates `z` on a specific sparsity parameter `k`. Only predictors in
`x` will be subject to sparsity constraint, unless `zkeep` keyword is specified.

# Arguments:
+ `y`: Phenotype vector or matrix. Should be an `Array{T, 1}` (single traits) or
    `Array{T, 2}` (multivariate Gaussian traits). For multivariate traits, each
    column of `y` should be a sample.
+ `x`: Genotype matrix (an `Array{T, 2}` or `SnpLinAlg`). For univariate
    analysis, samples are rows of `x`. For multivariate analysis, samples are
    columns of `x` (i.e. input `Transpose(x)` for `SnpLinAlg`)
+ `z`: Matrix of non-genetic covariates of type `Array{T, 2}` or `Array{T, 1}`.
    For univariate analysis, sample covariates are rows of `z`. For multivariate
    analysis, sample covariates are columns of `z`. If this is not specified, an
    intercept term will be included automatically. If `z` is specified, make sure
    the first column (row) is all 1s to represent the intercept.

# Optional Arguments:
+ `k`: Number of non-zero predictors. Can be a constant or a vector (for group IHT).
+ `J`: The number of maximum groups (set as 1 if no group infomation available)
+ `d`: Distribution of phenotypes. Specify `Normal()` for quantitative traits,
    `Bernoulli()` for binary traits, `Poisson()` or `NegativeBinomial()` for
    count traits, and `MvNormal()` for multiple quantitative traits.
+ `l`: A link function. The recommended link functions are `l=IdentityLink()` for
    quantitative traits, `l=LogitLink()` for binary traits, `l=LogLink()` for Poisson
    distribution, and `l=Loglink()` for NegativeBinomial distribution. For multivariate
    analysis, the choice of link does not matter.
+ `group`: vector storing (non-overlapping) group membership
+ `weight`: vector storing vector of weights containing prior knowledge on each SNP
+ `zkeep`: BitVector determining whether non-genetic covariates in `z` will be subject
    to sparsity constraint. `zkeep[i] = true` means covariate `i` will NOT be projected.
    Note covariates forced in the model are not subject to sparsity constraint `k`.
+ `est_r`: Symbol (`:MM`, `:Newton` or `:None`) to estimate nuisance parameters for negative binomial regression
+ `use_maf`: boolean indicating whether we want to scale projection with minor allele frequencies (see paper)
+ `debias`: boolean indicating whether we debias at each iteration
+ `verbose`: boolean indicating whether we want to print intermediate results
+ `tol`: used to track convergence
+ `max_iter`: is the maximum IHT iteration for a model to converge. Defaults to 200, or 100 for cross validation
+ `min_iter`: is the minimum IHT iteration before checking for convergence. Defaults to 5.
+ `max_step`: is the maximum number of backtracking per IHT iteration. Defaults 3
+ `io`: An `IO` object for displaying intermediate results. Default `stdout`.
+ `init_beta`: Whether to initialize beta values to univariate regression values.
    Currently only Gaussian traits can be initialized. Default `false`.
- `memory_efficient`: If `true,` it will cause ~1.1 times slow down but one only
    needs to store the genotype matrix (requiring 2np bits for PLINK binary files
    and `8np` bytes for other formats). If `memory_efficient=false`, one also need
    to store a `8nk` byte matrix where `k` is the sparsity levels.

# Output
+ An `IHTResult` (for single-trait analysis) or `mIHTResult` (for multivariate analysis).

# Group IHT
If `k` is a constant, then each group will have the same sparsity level. To run doubly
sparse IHT with varying group sparsities, construct `k` to be a vector where `k[i]`
indicates the max number of predictors for group `i`.
"""
function fit_iht(
    y         :: AbstractVecOrMat{T},
    x         :: AbstractMatrix{T},
    z         :: AbstractVecOrMat{T};
    k         :: Union{Int, Vector{Int}} = 10,
    J         :: Int = 1,
    d         :: Distribution = size(y, 2) > 1 ? MvNormal(T[]) : Normal(),
    l         :: Link = IdentityLink(),
    group     :: AbstractVector{Int} = Int[],
    weight    :: AbstractVector{T} = T[],
    zkeep     :: BitVector = trues(size(y, 2) > 1 ? size(z, 1) : size(z, 2)),
    est_r     :: Symbol = :None,
    use_maf   :: Bool = false,
    debias    :: Bool = false,
    verbose   :: Bool = true,          # print informative things to stdout
    tol       :: T = convert(T, 1e-4), # tolerance for tracking convergence
    max_iter  :: Int = 200,            # maximum IHT iterations
    min_iter  :: Int = 5,              # minimum IHT iterations
    max_step  :: Int = 3,              # maximum backtracking for each iteration
    io        :: IO = stdout,
    init_beta :: Bool = false,
    memory_efficient::Bool=true
    ) where T <: Float

    verbose && print_iht_signature(io)

    # first handle errors
    @assert J ≥ 0        "Value of J (max number of groups) must be nonnegative!\n"
    @assert max_iter ≥ 0 "Value of max_iter must be nonnegative!\n"
    @assert max_step ≥ 0 "Value of max_step must be nonnegative!\n"
    @assert tol > eps(T) "Value of global tol must exceed machine precision!\n"
    checky(y, d) # make sure response data y is in the form compatible with specified GLM
    check_group(k, group) # make sure sparsity parameter `k` is reasonable.
    !(typeof(d) <: NegativeBinomial) && est_r != :None &&
        error("Only negative binomial regression currently supports nuisance parameter estimation")
    typeof(x) <: AbstractSnpArray && error("x is a SnpArray! Please convert it to a SnpLinAlg first!")
    check_data_dim(y, x, z)
    if typeof(x) <: SnpLinAlg
        x.center || error("x is not centered! Please construct SnpLinAlg{Float64}(::SnpArray, center=true, scale=true)")
        x.scale || @warn("x is not scaled! We highly recommend `scale=true` in `SnpLinAlg` constructor")
        x.impute || @warn("x does not have impute flag! We highly recommend `impute=true` in `SnpLinAlg` constructor")
    end

    # initialize IHT variable
    v = initialize(x, z, y, J, k, d, l, group, weight, est_r, init_beta, zkeep,
        verbose=verbose, memory_efficient=memory_efficient)

    # print information
    verbose && print_parameters(io, k, d, l, use_maf, group, debias, tol, max_iter, min_iter)

    # fit IHT model
    tot_time, best_logl, mm_iter = fit_iht!(v, debias=debias, verbose=verbose,
        tol=tol, max_iter=max_iter, min_iter=min_iter, max_step=max_step, io=io)

    # compute phenotype's proportion of variation explained
    σ2 = pve(v)

    return IHTResult(tot_time, best_logl, mm_iter, σ2, v)
end

function fit_iht(
    y::AbstractVecOrMat{T},
    x::AbstractMatrix{T};
    kwargs...
    ) where T
    z = is_multivariate(y) ? ones(T, 1, size(y, 2)) : ones(T, length(y))
    return fit_iht(y, x, z; kwargs...)
end

"""
    fit_iht!(v; kwargs...)

Fits a IHT variable `v`.

# Arguments:
+ `v`: A properly initialized `mIHTVariable` or `IHTVariable`. Users should run [`fit_iht`](@ref)

# Optional Arguments:
+ `debias`: boolean indicating whether we debias at each iteration
+ `verbose`: boolean indicating whether we want to print results if model does not converge.
+ `tol`: used to track convergence
+ `max_iter`: is the maximum IHT iteration for a model to converge. Defaults to 200, or 100 for cross validation
+ `max_step`: is the maximum number of backtracking. Since l0 norm is not convex, we have no ascent guarantee
+ `io`: An `IO` object for displaying intermediate results. Default `stdout`.
"""
function fit_iht!(
    v         :: Union{mIHTVariable{T, M}, IHTVariable{T, M}};
    debias    :: Bool = false,
    verbose   :: Bool = true,          # print informative things
    tol       :: T = convert(T, 1e-4), # tolerance for tracking convergence
    max_iter  :: Int = 200,            # maximum IHT iterations
    min_iter  :: Int = 5,              # minimum IHT iterations
    max_step  :: Int = 3,              # maximum backtracking for each iteration
    io        :: IO = stdout
    ) where {T <: Float, M}

    #start timer
    start_time = time()

    # initialize constants
    mm_iter     = 0                 # number of iterations
    tot_time    = 0.0               # compute time *within* fit!
    next_logl   = typemin(T)        # loglikelihood
    best_logl   = typemin(T)        # best loglikelihood achieved
    η_step      = 0                 # counts number of backtracking steps for η

    # Begin 'iterative' hard thresholding algorithm
    for iter in 1:max_iter

        # notify and return current model if maximum iteration exceeded
        if iter ≥ max_iter
            best_logl = save_prev!(v, next_logl, best_logl)
            save_best_model!(v)
            mm_iter  = iter
            tot_time = time() - start_time
            verbose && printstyled(io, "Did not converge after $max_iter " *
                "iterations! IHT run time was " * string(tot_time) *
                " seconds\n", color=:red)
            break
        end

        # save values from previous iterate and update loglikelihood
        best_logl = save_prev!(v, next_logl, best_logl)

        # take one IHT step in positive score direction
        (η, η_step, next_logl) = iht_one_step!(v, next_logl, max_step)

        # perform debiasing if support didn't change
        debias && iter ≥ 5 && v.idx == v.idx0 && debias!(v)

        # track convergence
        # Note: estimated beta in first few iterations can be very small, so scaled_norm is very small
        # Thus we force IHT to iterate at least 5 times
        scaled_norm = check_convergence(v)
        progr = "Iteration $iter: loglikelihood = $next_logl, backtracks = $η_step, tol = $scaled_norm"
        verbose && println(io, progr)
        verbose && io != stdout && println(progr)
        if iter ≥ min_iter && scaled_norm < tol
            best_logl = save_prev!(v, next_logl, best_logl)
            save_best_model!(v)
            tot_time = time() - start_time
            mm_iter  = iter
            break
        end
    end

    return tot_time, best_logl, mm_iter
end

"""
Performs 1 iteration of the IHT algorithm, backtracking a maximum of `nstep` times.
We allow loglikelihood to potentially decrease to avoid bad boundary cases.
"""
function iht_one_step!(
    v::Union{IHTVariable{T, M}, mIHTVariable{T, M}},
    old_logl::T,
    nstep::Int
    ) where {T <: Float, M <: AbstractMatrix}

    # first calculate step size
    η = iht_stepsize!(v)

    # update b and c by taking gradient step v.b = P_k(β + ηv) where v is the score direction
    _iht_gradstep!(v, η)

    # update the linear predictors `xb`, `μ`, and residuals with the new proposed b
    update_xb!(v)
    update_μ!(v)

    # for multivariate IHT, also update precision matrix Γ = 1/n * (Y-BX)(Y-BX)'
    if typeof(v) <: mIHTVariable
        solve_Σ!(v)
    end

    # update r (nuisance parameter for negative binomial)
    if typeof(v) <: IHTVariable && v.est_r != :None
        v.d = mle_for_r(v)
    end

    # calculate current loglikelihood with the new computed xb and zc
    new_logl = loglikelihood(v)

    η_step = 0
    while _iht_backtrack_(new_logl, old_logl, η_step, nstep)

        # stephalving
        η /= 2

        # compute new loglikelihood after linesearch
        new_logl = backtrack!(v, η)

        # increment the counter
        η_step += 1
    end

    # compute score with the new mean
    score!(v)

    # check for finiteness before moving to the next iteration
    isnan(new_logl) && throw(error("Loglikelihood function is NaN, aborting..."))
    isinf(new_logl) && throw(error("Loglikelihood function is Inf, aborting..."))

    return η::T, η_step::Int, new_logl::T
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -12.548880784274985 <= q <= 15.06837792907793)

                     add_NL_constraint(m, :(softplus(0.4386422028907555 + 0.03049853723546203*softplus(-0.33806544253006265 + -0.5419050015026476*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.7963936693468932*softplus(0.35455126051631725 + -0.3141650319667155*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.19871150900740808*softplus(-0.9712553762218641 + -0.9788234903921222*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) + softplus(0.039843038969352484 + 0.2885956651060697*softplus(-0.33806544253006265 + -0.5419050015026476*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + 0.8477392827604406*softplus(0.35455126051631725 + -0.3141650319667155*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.38122333260621444*softplus(-0.9712553762218641 + -0.9788234903921222*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) + softplus(0.5234531332153054 + -0.26529753551452195*softplus(-0.33806544253006265 + -0.5419050015026476*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + 0.29155449146294243*softplus(0.35455126051631725 + -0.3141650319667155*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.357141816330512*softplus(-0.9712553762218641 + -0.9788234903921222*softplus(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*softplus(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*softplus(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Documenter, ImageReconstruction

makedocs(;
    modules = [ImageReconstruction],
    format = Documenter.HTML(),
    pages = ["Home" => "index.md"],
    repo = "https://github.com/JuliaImages/ImageReconstruction.jl/blob/{commit}{path}#L{line}",
    sitename = "ImageReconstruction.jl",
    assets = String[],
)

deploydocs(; repo = "github.com/JuliaImages/ImageReconstruction.jl")


text:
using GateDecompositions
using Test

@testset "GateDecompositions.jl" begin

      U4 = haar_unitary(4)

      @test test_KAK_decomp(U4) ≈ 0.0

end


text:
"""
    nlsolve!(nlsolver::NLSolver, nlcache::Union{NLFunctionalCache,NLAndersonCache,NLFunctionalConstantCache,NLAndersonConstantCache}, integrator)

Perform functional iteration that is used by implicit methods.

It solves

```math
G(z) = dt⋅f(tmp + γ⋅z, p, t + c⋅h)
z = G(z)
```

by iterating

```math
zᵏ⁺¹ = G(zᵏ),
```

where `dt` is the step size and `γ` is a constant.

It returns the tuple `(z, η, iter, fail_convergence)`, where `z` is the solution, `η` is
used to measure the iteration error (see [^HW96]), `iter` is the number of iterations, and
`fail_convergence` reports whether the algorithm succeeded.

[^HW96]: Ernst Hairer and Gerhard Wanner, "Solving Ordinary Differential
Equations II, Springer Series in Computational Mathematics. ISBN
978-3-642-05221-7. Section IV.8.
[doi:10.1007/978-3-642-05221-7](https://doi.org/10.1007/978-3-642-05221-7)
"""
@muladd function nlsolve!(nlsolver::NLSolver, nlcache::Union{NLFunctionalConstantCache,NLAndersonConstantCache}, integrator)
  @unpack t,dt,uprev,u,p = integrator
  @unpack z,tmp,κtol,c,γ,max_iter = nlsolver

  if nlcache isa NLAndersonConstantCache
    @unpack Δz₊s,Q,R,γs,aa_start,droptol = nlcache
  end

  # precalculations
  mass_matrix = integrator.f.mass_matrix
  f = nlsolve_f(integrator)
  tstep = t + c*dt
  η = nlsolver.ηold
  if nlcache isa NLAndersonConstantCache
    history = 0
    max_history = length(Δz₊s)
  end

  # fixed point iteration
  local ndz
  if nlcache isa NLAndersonConstantCache
    local Δdz, dzold, z₊old # cache variables
  end
  fail_convergence = true
  iter = 0
  while iter < max_iter
    iter += 1
    integrator.destats.nnonliniter += 1

    # evaluate function
    u = @. tmp + γ*z
    if mass_matrix == I
      z₊ = dt .* f(u, p, tstep)
      dz = z₊ .- z
    else
      mz = _reshape(mass_matrix * _vec(z), axes(z))
      dz = dt .* f(u, p, tstep) .- mz
      z₊ = z .+ dz
    end
    integrator.destats.nf += 1

    # compute norm of residuals
    iter > 1 && (ndzprev = ndz)
    ndz = integrator.opts.internalnorm(dz, tstep)

    # check divergence (not in initial step)
    if iter > 1
      θ = ndz / ndzprev
      if θ > 1 || ndz * θ^(max_iter - iter) > κtol * (1 - θ)
        # fixed-point iteration diverges
        break
      end
    end

    # update iterate
    z = z₊

    # check stopping criterion
    iter > 1 && (η = θ / (1 - θ))
    if η * ndz < κtol && (iter > 1 || iszero(ndz))
      # fixed-point iteration converges
      fail_convergence = false
      break
    end

    # perform Anderson acceleration
    if nlcache isa NLAndersonConstantCache && iter < max_iter
      if iter == aa_start
        # update cached values for next step of Anderson acceleration
        dzold = dz
        z₊old = z₊
      elseif iter > aa_start
        # increase size of history
        history += 1

        # remove oldest history if maximum size is exceeded
        if history > max_history
          # circularly shift differences of z₊
          for i in 1:(max_history-1)
            Δz₊s[i] = Δz₊s[i + 1]
          end

          # delete left-most column of QR decomposition
          qrdelete!(Q, R, max_history)

          # update size of history
          history = max_history
        end

        # update history of differences of z₊
        Δz₊s[history] = @. z₊ - z₊old

        # replace/add difference of residuals as right-most column to QR decomposition
        qradd!(Q, R, _vec(dz .- dzold), history)

        # update cached values
        dzold = dz
        z₊old = z₊

        # define current Q and R matrices
        Qcur, Rcur = view(Q, :, 1:history), UpperTriangular(view(R, 1:history, 1:history))

        # check condition (TODO: incremental estimation)
        if droptol !== nothing
          while cond(R) > droptol && history > 1
            qrdelete!(Q, R, history)
            history -= 1
            Qcur, Rcur = view(Q, :, 1:history), UpperTriangular(view(R, 1:history, 1:history))
          end
        end

        # solve least squares problem
        γscur = view(γs, 1:history)
        ldiv!(Rcur, mul!(γscur, Qcur', _vec(dz)))
        integrator.destats.nsolve += 1

        # update next iterate
        for i in 1:history
          z = @. z - γs[i] * Δz₊s[i]
        end

        # update norm of residuals
        ndz = integrator.opts.internalnorm(z .- z₊ .+ dz, tstep)
      end
    end
  end
  if fail_convergence
    integrator.destats.nnonlinconvfail += 1
  end
  integrator.force_stepfail = fail_convergence
  z, η, iter, fail_convergence
end

@muladd function nlsolve!(nlsolver::NLSolver, nlcache::Union{NLFunctionalCache,NLAndersonCache}, integrator)
  @unpack t,dt,uprev,u,p = integrator
  @unpack z,dz,tmp,ztmp,k,κtol,c,γ,max_iter = nlsolver

  if nlcache isa NLFunctionalCache
    @unpack z₊ = nlcache
  else
    @unpack z₊,dzold,z₊old,Δz₊s,Q,R,γs,aa_start,droptol = nlcache
  end

  # precalculations
  vecztmp = vec(ztmp); vecz = vec(z); vecz₊ = vec(z₊)
  mass_matrix = integrator.f.mass_matrix
  f = nlsolve_f(integrator)
  tstep = t + c*dt
  η = nlsolver.ηold
  if nlcache isa NLAndersonCache
    history = 0
    max_history = length(Δz₊s)
  end

  # fixed-point iteration without Newton
  local ndz
  fail_convergence = true
  iter = 0
  while iter < max_iter
    iter += 1
    integrator.destats.nnonliniter += 1

    # evaluate function
    @. u = tmp + γ*z
    f(k, u, p, tstep)
    integrator.destats.nf += 1
    if mass_matrix == I
      @. z₊ = dt*k
      @. dz = z₊ - z
    else
      mul!(vecztmp, mass_matrix, vecz)
      @. dz = dt*k - ztmp
      @. z₊ = z + dz
    end

    # compute norm of residuals
    iter > 1 && (ndzprev = ndz)
    ndz = integrator.opts.internalnorm(dz, tstep)

    # check divergence (not in initial step)
    if iter > 1
      θ = ndz / ndzprev
      if θ > 1 || ndz * θ^(max_iter - iter) > κtol * (1 - θ)
        # fixed-point iteration diverges
        break
      end
    end

    # update iterate
    @. z = z₊

    # check stopping criterion
    iter > 1 && (η = θ / (1 - θ))
    if η * ndz < κtol && (iter > 1 || iszero(ndz))
      # fixed-point iteration converges
      fail_convergence = false
      break
    end

    # perform Anderson acceleration
    if nlcache isa NLAndersonCache && iter < max_iter
      if iter == aa_start
        # update cached values for next step of Anderson acceleration
        @. dzold = dz
        @. z₊old = z₊
      elseif iter > aa_start
        # increase size of history
        history += 1

        # remove oldest history if maximum size is exceeded
        if history > max_history
          # circularly shift differences of z₊
          ptr = Δgs[1]
          for i in 1:(max_history-1)
            Δz₊s[i] = Δz₊s[i + 1]
          end
          Δz₊s[max_history] = ptr

          # delete left-most column of QR decomposition
          qrdelete!(Q, R, max_history)

          # update size of history
          history = max_history
        end

        # update history of differences of z₊
        @. Δz₊s[history] = z₊ - z₊old

        # replace/add difference of residuals as right-most column to QR decomposition
        @. dzold = dz - dzold
        qradd!(Q, R, vec(dzold), history)

        # update cached values
        @. dzold = dz
        @. z₊old = z₊

        # define current Q and R matrices
        Qcur, Rcur = view(Q, :, 1:history), UpperTriangular(view(R, 1:history, 1:history))

        # check condition (TODO: incremental estimation)
        if droptol !== nothing
          while cond(R) > droptol && history > 1
            qrdelete!(Q, R, history)
            history -= 1
            Qcur, Rcur = view(Q, :, 1:history), UpperTriangular(view(R, 1:history, 1:history))
          end
        end

        # solve least squares problem
        γscur = view(γs, 1:history)
        ldiv!(Rcur, mul!(γscur, Qcur', vec(dz)))
        integrator.destats.nsolve += 1

        # update next iterate
        for i in 1:history
          @. z = z - γs[i] * Δz₊s[i]
        end

        # update norm of residuals
        @. dz = z - z₊ + dz
        ndz = integrator.opts.internalnorm(dz, tstep)
      end
    end
  end
  if fail_convergence
    integrator.destats.nnonlinconvfail += 1
  end
  integrator.force_stepfail = fail_convergence
  z, η, iter, fail_convergence
end


text:
############################################################
# IterativeSolvers support #################################
############################################################

Adivtype(A::joAbstractLinearOperator{DDT,RDT}, b::AbstractVector) where {DDT,RDT} = DDT


text:
export fetchval

const SHOULD_SEND = Dict(:PATCH => true, :POST => true, :PUT => true)

# Lists of events fired by sending a request to each endpoint.
const EVENTS_FIRED = Dict(
    :DELETE => [
        (r"^/channels/\d+$", ChannelDelete),
        (r"^/channels/\d+$", ChannelUpdate),
        (r"^/channels/\d+/messages/\d+/reactions/.+/@me$", MessageReactionRemove),
        (r"^/channels/\d+/messages/\d+/reactions/.+/\d+$", MessageReactionRemove),
        (r"^/channels/\d+/messages/\d+/reactions$", MessageReactionRemoveAll),
        (r"^/channels/\d+/messages/\d+$", MessageDelete),
        (r"^/channels/\d+/messages/bulk-delete$", MessageDeleteBulk),
        (r"^/channels/\d+/permissions/\d+$", ChannelUpdate),
        (r"^/channels/\d+/pins/\d+$", ChannelPinsUpdate),
        (r"^/channels/\d+/pins/\d+$", MessageUpdate),
        (r"^/channels/\d+/recipients/\d+$", ChannelUpdate),  # Untested.
        (r"^/guilds/\d+/emojis/\d+$", GuildEmojisUpdate),
        (r"^/guilds/\d+$", GuildDelete),
        (r"^/guilds/\d+/members/\d+/roles/\d+$", GuildMemberUpdate),
        (r"^/guilds/\d+/members/\d+/$", GuildMemberRemove),
        (r"^/guilds/\d+/bans/\d+$", GuildBanRemove),
        (r"^/guilds/\d+/roles/\d+$", GuildRoleDelete),
        (r"^/guilds/\d+/integrations$", GuildIntegrationsUpdate),
        (r"^/users/@me/guilds/\d+$", GuildMemberRemove),
        (r"^/users/@me/guilds/\d+$", PresenceUpdate),
        (r"^/webhooks/\d+$", WebhooksUpdate),
        (r"^/webhooks/\d+/.+$", WebhooksUpdate),
    ],
    :POST => [
        (r"^/channels/\d+/messages$", MessageCreate),
        (r"^/channels/\d+/typing$", TypingStart),
        (r"^/guilds/\d+/emojis$", GuildEmojisUpdate),
        (r"^/guilds$", GuildCreate),
        (r"^/guilds/\d+/channels$", ChannelCreate),
        (r"^/guilds/\d+/roles$", GuildRoleCreate),
        (r"^/guilds/\d+/prune$", GuildMemberRemove),
        (r"^/guilds/\d+/integrations$", GuildIntegrationsUpdate),
        (r"^/users/@me/channels$", ChannelCreate),
        (r"^/channels/\d+/webhooks$", WebhooksUpdate),
    ],
    :PUT => [
        (r"^/channels/\d+$", ChannelUpdate),
        (r"^/channels/\d+/messages/\d+/reactions/.+/@me$", MessageReactionAdd),
        (r"^/channels/\d+/permissions/\d+$", ChannelUpdate),
        (r"^/channels/\d+/pins/\d+$", ChannelPinsUpdate),
        (r"^/channels/\d+/pins/\d+$", MessageUpdate),
        (r"^/channels/\d+/recipients/\d+$", ChannelUpdate),  # Untested.
        (r"^/guilds/\d+/members/\d+$", GuildMemberAdd),
        (r"^/guilds/\d+/members/\d+/roles/\d+$", GuildMemberUpdate),
        (r"^/guilds/\d+/bans/\d+$", GuildBanAdd),
    ],
    :PATCH => [
        (r"^/channels/\d+$", ChannelUpdate),
        (r"^/channels/\d+/messages/\d+$", MessageUpdate),
        (r"^/guilds/\d+/emojis/\d+$", GuildEmojisUpdate),
        (r"^/guilds/\d+$", GuildUpdate),
        (r"^/guilds/\d+/channels$", ChannelUpdate),
        (r"^/guilds/\d+/members/\d+$", GuildMemberUpdate),
        (r"^/guilds/\d+/members/@me/nick$", GuildMemberUpdate),
        (r"^/guilds/\d+/roles$", GuildRoleUpdate),
        (r"^/guilds/\d+/roles/\d+$", GuildRoleUpdate),
        (r"^/guilds/\d+/integrations$", GuildIntegrationsUpdate),
        (r"^/users/@me$", UserUpdate),
        (r"^/webhooks/\d+$", WebhooksUpdate),
        (r"^/webhooks/\d+/.+$", WebhooksUpdate),
    ],
)

"""
A wrapper around a response from the REST API. Every function which wraps a Discord REST
API endpoint returns a `Future` which will contain a value of this type. To retrieve the
`Response` from the `Future`, use `fetch`. See also: [`fetchval`](@ref).

## Fields
- `val::Nullable{T}`: The object contained in the HTTP response. For example, for a
  call to [`get_channel_message`](@ref), this value will be a [`Message`](@ref).
- `ok::Bool`: The state of the request. If `true`, then it is safe to access `val`.
- `http_response::Nullable{HTTP.Messages.Response}`: The underlying HTTP response, if
  a request was made.
- `exception::Nullable{Exception}`: The caught exception, if one is thrown.

## Examples
Multiple API calls which immediately return `Future`s and can be awaited:
```julia
futures = map(i -> create_message(c, channel_id; content=string(i)), 1:10);
other_work_here()
resps = fetch.(futures)
```
Skipping error checks and returning the value directly:
```julia
guild = fetchval(create_guild(c; name="foo"))
```
"""
struct Response{T}
    val::Nullable{T}
    ok::Bool
    http_response::Nullable{HTTP.Messages.Response}
    exception::Nullable{Exception}
end

# HTTP response with no body.
function Response{Nothing}(r::HTTP.Messages.Response)
    return Response{Nothing}(nothing, r.status < 300, r, nothing)
end

# HTTP response with body (maybe).
function Response{T}(c::Client, r::HTTP.Messages.Response) where T
    r.status == 204 && return Response{T}(nothing, true, r, nothing)
    r.status >= 300 && return Response{T}(nothing, false, r, nothing)

    data = if HTTP.header(r, "Content-Type") == "application/json"
        JSON.parse(String(copy(r.body)); dicttype=Dict{Symbol, Any})
    else
        copy(r.body)
    end

    val, e = tryparse(c, T, data)
    return Response{T}(val, e === nothing, r, e)
end

# HTTP request with no expected response body.
function Response(
    c::Client,
    method::Symbol,
    endpoint::AbstractString;
    body="",
    kwargs...
)
    return Response{Nothing}(c, method, endpoint; body=body, kwargs...)
end

# HTTP request.
function Response{T}(
    c::Client,
    method::Symbol,
    endpoint::AbstractString,
    cachetest::Function=alwaystrue;
    headers=Dict(),
    body=HTTP.nobody,
    kwargs...
) where T
    f = Future()

    @async begin
        kws = (
            channel=cap("channels", endpoint),
            guild=cap("guilds", endpoint),
            message=cap("messages", endpoint),
            user=cap("users", endpoint),
        )

        # Retrieve the value from the cache, maybe.
        if c.use_cache && method === :GET
            val = get(c.state, T; kws...)
            if val !== nothing
                try
                    # Only use the value if it satisfies the cache test.
                    if cachetest(val) === true
                        put!(f, Response{T}(val, true, nothing, nothing))
                        return
                    end
                catch
                    # TODO: Log this?
                end
            end
        end

        # Prepare the request.
        url = "$DISCORD_API/v$(c.version)$endpoint"
        isempty(kwargs) || (url *= "?" * HTTP.escapeuri(kwargs))
        headers = merge(Dict(
            "User-Agent" => "Discord.jl $DISCORD_JL_VERSION",
            "Content-Type" => "application/json",
            "Authorization" => c.token,
        ), Dict(headers))
        args = [method, url, headers]
        if get(SHOULD_SEND, method, false)
            push!(args, headers["Content-Type"] == "application/json" ? json(body) : body)
        end

        # Queue a job to be run within the rate-limiting constraints.
        enqueue!(c.limiter, method, endpoint) do
            http_r = nothing

            try
                # Make an HTTP request, and generate a Response.
                # If we hit an upstream rate limit, return the response immediately.
                http_r = HTTP.request(args...; status_exception=false)
                http_r.status == 429 && return http_r
                r = Response{T}(c, http_r)

                if r.ok && r.val !== nothing && should_put(c, method, endpoint)
                    # Cache the value, and also replace any missing data
                    # in the response with data from the cache.
                    r = @set r.val = put!(c.state, r.val; kws...)
                end

                # Store the successful Response to the Future.
                put!(f, r)
            catch e
                kws = logkws(c; conn=undef, exception=(e, catch_backtrace()))
                @error "Creating response failed" kws...
                put!(f, Response{T}(nothing, false, http_r, e))
            end

            return http_r  # Return the HTTP response to update the rate limits.
        end
    end

    return f
end

Base.eltype(::Response{T}) where T = T
Base.eltype(::Type{Response{T}}) where T = T

"""
    fetchval(f::Future{Response{T}}) -> Nullable{T}

Shortcut for `fetch(f).val`: Fetch a [`Response`](@ref) and return its value. Note that
there are no guarantees about the response's success and the value being returned, and it
discards context that can be useful for debugging, such as HTTP responses and caught
exceptions.
"""
fetchval(f::Future) = fetch(f).val

# Capture an ID from a string.
function cap(path::AbstractString, s::AbstractString)
    m = match(Regex("/$path/(\\d+)"), s)
    return m === nothing ? nothing : parse(Snowflake, first(m.captures))
end

# Determine whether we should cache something for the request.
function should_put(c::Client, method::Symbol, endpoint::AbstractString)
    # TODO: Maybe we could use the type of the Response value to be more granular.
    # e.g. if we got a Message, then only check handlers for types with a Message field.
    c.use_cache || return false
    isopen(c) || return true
    method === :GET && return true

    Ts = map(
        last,
        Iterators.filter(t -> match(t[1], endpoint) !== nothing, EVENTS_FIRED[method]),
    )
    return !all(T -> hasdefault(c, T), Ts)
end

include(joinpath("endpoints", "endpoints.jl"))
include(joinpath("crud", "crud.jl"))


text:
using ArgParse, Core
include("../src/debruijnnets.jl");
include("../src/sampling.jl");
include("../src/motifs.jl")
include("../src/randomwalks.jl")
include("../src/hypa.jl")

"""
Run a hypa simulation
"""
function hypa_sample(fo, ko, ko_map, fo_map, k, num_samples, empirical_motifs)
    ko_rev_map = Dict(val=>key for (key,val) in ko_map)
    # Get the adjacency and Xi matrices
    adjacency = adjacency_matrix(ko)
    pvals, Xi = hypa(k, ko, ko_map, fo, fo_map)
    sampled_counts = Dict(m=>Dict("frequency"=>empirical_motifs[m], "samples"=>zeros(num_samples)) for m in keys(empirical_motifs))
    # TODO: For some reason this is not thread safe and hangs if I use Threads.@threads
    for run in range(1, num_samples)
        # Draw a sample
        Core.println("Sampling adjacency.")
        sampled_adj = draw_sample(adjacency, Xi)
        Core.println("Sampled.")
        # Transform into walks
        cidxs = findall(>(0), sampled_adj)
        walks = Vector{Tuple}(undef, length(cidxs))
        weights = Vector{Int64}(undef, length(cidxs))
        widx = 1
        for cidx in cidxs
            u, v = Tuple(cidx)
            walk = (ko_rev_map[u]...,ko_rev_map[v][end])
            walks[widx] = walk
            weight = sampled_adj[cidx]
            weights[widx] = weight
            widx += 1
        end
        # Count motifs
        motifs = count_motifs(walks, weights)
        for (motif, val) in motifs
            sampled_counts[motif]["samples"][run] = val
        end
    end
    return sampled_counts
end

"""
Run a Biased DeBruijn graph simulation.
"""
function bdg(fo, ko, ko_map, fo_map, k, M, num_samples, empirical_motifs)
    # Get all of the walks from the korder graph
    all_walks, weights = walks_from_edges(ko, ko_map, fo_map)
    sampled_counts = Dict(m=>Dict("frequency"=>empirical_motifs[m], "samples"=>zeros(num_samples)) for m in keys(empirical_motifs))
    Threads.@threads for run in range(1, num_samples)
        # Sample M walks at a time
        walks = sample(all_walks, M)
        # Count motifs
        motifs = count_motifs(walks)
        for (motif, val) in motifs
            sampled_counts[motif]["samples"][run] = val
        end
    end
    return sampled_counts
end

"""
Run a line-graph sampling simulation.
"""
function lgs(fo, k, M, num_samples,  walks_per_node, empirical_motifs)
    println("wpn: $walks_per_node.")
    if walks_per_node < 0
        println("Computing all walks for every node.")
        all_kedge_walks = get_all_walks(fo, k)
        println("Done with computation.")
    end

    # Sample a bunch of times
    sampled_counts = Dict(m=>Dict("frequency"=>empirical_motifs[m], "samples"=>zeros(num_samples)) for m in keys(empirical_motifs))
    Threads.@threads for run in range(1, num_samples)
        if walks_per_node > 0
            walks = uniform_walk_sample(fo, k, M, walks_per_node, true, true)
        else
            walks = sample(all_kedge_walks, M)
        end
        # Count motifs
        motifs = count_motifs(walks)
        for (motif, val) in motifs
            sampled_counts[motif]["samples"][run] = val
        end
    end
    return sampled_counts
end

"""
    Run a random-walk sampling simulation.
"""
function rw(fo, k::Integer, M::Integer , num_samples::Integer, empirical_motifs,
        bias_nodes::Bool=false, weighted::Bool=false)
    nodes, A, A_bin, path_counts = filter_nodes(fo, k)
    node_probabilities = compute_node_probs(nodes, path_counts)
    # Sample a bunch of times
    sampled_counts = Dict(m=>Dict("frequency"=>empirical_motifs[m], "samples"=>zeros(num_samples)) for m in keys(empirical_motifs))
    Threads.@threads for run in range(1, num_samples)
        walks = random_walks(fo, k, M, nodes, node_probabilities, bias_nodes, weighted)
        # Count motifs
        motifs = count_motifs(walks)
        for (motif, val) in motifs
            sampled_counts[motif]["samples"][run] = val
        end
    end
    return sampled_counts
end

function write_counts(sampled_counts, output_file)
    open(output_file, "w") do file
        for (m, mdict) in sampled_counts
            out_str = join(m, ',') * '|' * string(mdict["frequency"]) * '|' * join(mdict["samples"], ",") * "\n"
            write(file, out_str)
        end
    end
    return nothing
end

function parse_commandline(args)
    s = ArgParseSettings()
    @add_arg_table s begin
	"--input_file", "-i"
        help = "Full path to input file"
        #arg_type = String
	"--frequency", "-f"
	    help = "Frequencies in ngram?"
	    action = :store_true
        "-k"
	    help = "Order of the motifs."
	    arg_type = Int
	"-w", "--walks_per_node"
	    help = "Walks per node. -1 for all."
            arg_type = Int
	"-n", "--num_samples"
	    help = "Number of samples to draw."
	    arg_type = Int
	    default = 10
	"-e", "--ensemble"
	    help = "Ensemble to use."
	    arg_type = String
    end
return parse_args(args, s)
end

arguments = parse_commandline(ARGS)

println(arguments)
input_filepath = arguments["input_file"]
println(input_filepath)
frequency = arguments["frequency"]
println(frequency)
k = arguments["k"]
walks_per_node = arguments["walks_per_node"]
num_samples = arguments["num_samples"]
ensemble = arguments["ensemble"]
num_cpus = Threads.nthreads()
splitpath = split(input_filepath, '/')
ngram_filename = splitpath[end]
weighted = false
if ensemble == "rw-w"
    weighted = true
end
println("k: $k")
fo, fo_map, ko, ko_map = from_ngram(input_filepath, frequency, k);
if has_self_loops(fo)
   println("fo has selfloops..")
end
# Count motifs
walks, walk_freqs = walks_from_edges(ko, ko_map, fo_map)
empirical_motifs = count_motifs(walks, walk_freqs)
println(empirical_motifs)
M = sum(values(empirical_motifs))
if ensemble == "lg-s"
    sampled_counts = lgs(fo, k, M, num_samples, walks_per_node, empirical_motifs)
    output_filename = "../../debruijn-nets/results/motifs/$(ngram_filename)_k-$(k)_e-$(ensemble)_bn_ww_wpn-$(walks_per_node).csv"
elseif ensemble == "rw-uw" || ensemble == "rw-w"
    sampled_counts = rw(fo, k, M, num_samples, empirical_motifs, false, weighted)
    output_filename = "../../debruijn-nets/results/motifs/$(ngram_filename)_k-$(k)_e-$(ensemble).csv"
elseif ensemble == "bdg"
    sampled_counts = bdg(fo, ko, ko_map, fo_map, k, M, num_samples, empirical_motifs)
    output_filename = "../../debruijn-nets/results/motifs/$(ngram_filename)_k-$(k)_e-$(ensemble).csv"
elseif ensemble == "hypa"
    println("WARNING: hypa implementation has an over/underflow issue that has not been resolved. Use at your own risk.")
    sampled_counts = hypa_sample(fo, ko, ko_map, fo_map, k, num_samples, empirical_motifs)
    output_filename = "../../debruijn-nets/results/motifs/$(ngram_filename)_k-$(k)_e-h.csv"
end

println(sampled_counts)
write_counts(sampled_counts, output_filename)


text:
include("run_sim.jl")

T = 100000
p_E_range = 0.5
p_H_range = 0.02
λ_E_range = 0.5:0.5:(10-0.5)
λ_H_range = 2
transitions_list = [tr_chains_d1, tr_chains_d10,
                    tr_chains_d50, tr_chains_d100]
table1 = run_multiple_sim(λ_E_range, λ_H_range, p_E_range, p_H_range, transitions_list);
writetable("csv/fig_4.1.csv", table1)

T = 100000
p_E_range = 0.5
p_H_range = 0.02
λ_H_range = 0.5:0.5:(10-0.5)
λ_E_range = 2
transitions_list = [tr_chains_d1, tr_chains_d10,
                    tr_chains_d50, tr_chains_d100]
table2 = run_multiple_sim(λ_E_range, λ_H_range, p_E_range, p_H_range, transitions_list);
writetable("csv/fig_4.2.csv", table2)

T = 1000000
p_E_range = 0.5
p_H_range = [0.1, 0.05, 0.01, 0.005, 0.001, 0.0005, 0.0001, 0.00005, 0.00001]
λ_E_range = 1
λ_H_range = 2
transitions_list = [tr_chains_d1, tr_chains_d10,
                    tr_chains_d50, tr_chains_d100]
table3 = run_multiple_sim(λ_E_range, λ_H_range, p_E_range, p_H_range, transitions_list);

table3[:norm_d1] = table3[:tr_chains_d1] .* table3[:p_H]
table3[:norm_d10] = table3[:tr_chains_d10] .* table3[:p_H]
table3[:norm_d50] = table3[:tr_chains_d50] .* table3[:p_H]
table3[:norm_d100] = table3[:tr_chains_d100] .* table3[:p_H]
table3[:log_p_H] = - log(table3[:p_H])

writetable("csv/fig_4.3.csv", table3)

T = 1000000
p_E_range = 0.1:0.1:1
p_H_range = 0.01
λ_E_range = 1
λ_H_range = 2
transitions_list = [tr_chains_d1, tr_chains_d10,
                    tr_chains_d50, tr_chains_d100]
table4 = run_multiple_sim(λ_E_range, λ_H_range, p_E_range, p_H_range, transitions_list);

writetable("csv/fig_4.4.csv", table4)


text:
function soc1_test()
    #=
        max 0x + 1y + 1z
    s.t.
        x == 1
        x >= ||(y,z)||
    =#
    model = TestModel{Float64}()

    x, y, z = MOI.add_variables(model, 3)

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [y, z]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    ceq = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
            [-1.0],
        ),
        MOI.Zeros(1),
    )

    vov = MOI.VectorOfVariables([x, y, z])

    csoc = MOI.add_constraint(model, vov, MOI.SecondOrderCone(3))

    return model
end

function soc2_test()
    #=
        max 0x + 1y + 1z
      s.t.
        x == 1
        x >= ||(y,z)||
    =#
    model = TestModel{Float64}()

    x, y, z = MOI.add_variables(model, 3)

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [y, z]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    ceq = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
            [-1.0],
        ),
        MOI.Zeros(1),
    )

    vov = MOI.VectorOfVariables([x, y, z])

    csoc = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{Float64}(vov),
        MOI.SecondOrderCone(3),
    )

    return model
end

function soc3_test()
    #=
        min  x
    s.t.
        y ≥ 1/√2
        x² + y² ≤ 1

    in conic form:
        min  x
    s.t.
        -1/√2 + y ∈ R₊
        1 - t ∈ {0}
        (t,x,y) ∈ SOC₃
    =#
    model = TestModel{Float64}()

    x, y, t = MOI.add_variables(model, 3)

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    cnon = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y))],
            [-1 / √2],
        ),
        MOI.Nonnegatives(1),
    )
    ceq = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, t))],
            [1.0],
        ),
        MOI.Zeros(1),
    )
    csoc = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 2, 3],
                MOI.ScalarAffineTerm.(1.0, [t, x, y]),
            ),
            zeros(3),
        ),
        MOI.SecondOrderCone(3),
    )

    return model
end

function soc4_test()
    #=
        min  x
    s.t.
        y ≥ 1/√2
        x² + y² ≤ 1

    in conic form:
        min  x
    s.t.
        1/√2 - y ∈ R₋
        1 - t ∈ {0}
        (t,x,y) ∈ SOC₃
    =#
    model = TestModel{Float64}()

    x, y, t = MOI.add_variables(model, 3)

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    cnon = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, y))],
            [1 / √2],
        ),
        MOI.Nonpositives(1),
    )

    ceq = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, t))],
            [1.0],
        ),
        MOI.Zeros(1),
    )
    csoc = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 2, 3],
                MOI.ScalarAffineTerm.(1.0, [t, x, y]),
            ),
            zeros(3),
        ),
        MOI.SecondOrderCone(3),
    )

    return model
end

function soc5_test()
    #=
        min x
    s.t.
        y ≥ 2
        x ≤ 1
        |y| ≤ x

    in conic form:
        min x
    s.t.
        -2 + y ∈ R₊
        -1 + x ∈ R₋
        (x,y) ∈ SOC₂
    =#
    model = TestModel{Float64}()

    x, y = MOI.add_variables(model, 2)

    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y))],
            [-2.0],
        ),
        MOI.Nonnegatives(1),
    )
    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
            [-1.0],
        ),
        MOI.Nonpositives(1),
    )
    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.(1.0, [x, y])),
            zeros(2),
        ),
        MOI.SecondOrderCone(2),
    )

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function soc6_test()
    #=
        min 0x[1] - 2x[2] - 1x[3]
    s.t.
        x[1]                    == 1
            x[2]   - x[4]       == 0
               x[3]      - x[5] == 0
        x[1] >= ||(x[4],x[5])||

    in conic form:
        min  c^Tx
    s.t.
        Ax + b ∈ {0}₃
        (x[1],x[4],x[5]) ∈ SOC₃
    =#
    model = TestModel{Float64}()

    x = MOI.add_variables(model, 5)

    c1 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 2, 3, 2, 3],
                MOI.ScalarAffineTerm.([1.0, 1.0, 1.0, -1.0, -1.0], x),
            ),
            [-1.0, 0.0, 0.0],
        ),
        MOI.Zeros(3),
    )
    c2 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x[1], x[4], x[5]]),
        MOI.SecondOrderCone(3),
    )

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([0.0, -2.0, -1.0, 0.0, 0.0], x),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end


text:
# # "Gaussian advection-diffusion" Spatial resolution convergence test

using PyPlot

# Define a few utilities for running tests and unpacking and plotting results

include("ConvergenceTests/ConvergenceTests.jl") # we use the GaussianAdvectionDiffusion module here.

run_test = ConvergenceTests.OneDimensionalGaussianAdvectionDiffusion.run_test

""" Run advection-diffusion test for all Nx in resolutions. """
function run_convergence_test(κ, U, resolutions...)

    # Determine save time-step
           Lx = 2.5
    stop_time = 0.25
            h = Lx / maximum(resolutions)
           Δt = min(0.1 * h / U, 0.01 * h^2 / κ)

    # Adjust time-step
    stop_iteration = round(Int, stop_time / Δt)
                Δt = stop_time / stop_iteration

    # Run the tests
    results = [run_test(Nx=Nx, Δt=Δt, stop_iteration=stop_iteration, U=U, κ=κ, width=0.1) for Nx in resolutions]

    return results
end

#####
##### Run test
#####

Nx = 2 .^ (6:8) # N = 64 through N = 256
diffusion_results = run_convergence_test(1e-1, 0, Nx...)
advection_results = run_convergence_test(1e-6, 3, Nx...)
advection_diffusion_results = run_convergence_test(1e-2, 1, Nx...)

#####
##### Plot solution and error profile
#####

all_results = (diffusion_results, advection_results, advection_diffusion_results)
names = ("diffusion only",  "advection only",  "advection-diffusion")
linestyles = ("-", "--", ":")
specialcolors = ("xkcd:black", "xkcd:indigo", "xkcd:wine red")

# Solutions
close("all")

fig, axs = subplots(nrows=2, figsize=(12, 6), sharex=true)

ConvergenceTests.OneDimensionalUtils.plot_solutions!(axs, all_results, names, linestyles, specialcolors)

savefig("figs/gaussian_advection_diffusion_solutions.png", dpi=480)

# Error profile
fig, axs = subplots()

ConvergenceTests.OneDimensionalUtils.plot_error_convergence!(axs, Nx, all_results, names)

savefig("figs/gaussian_advection_diffusion_error_convergence.png", dpi=480)


text:
# Usage instruction:
# julia quantum_dot.jl groundstate # calculates groundstate at t0 = 0
# julia quantum_dot.jl ref [tol] # tol::Float64, calc ref solution with specified tol for AdaptiveAdamsLawson
# julia quantum_dot.jl {conv | stab} [method] [num_steps] # run experiment for equidistant method as listed below
# conv compares the result with a reference solution and prints the norm distance function
# stab calculates until t = 70 without external potential
#
# Parameters are changed as global variables otherwise
# Experiment results are being saved to subdirectories, please review those lines and mkdir before running!
# Saved data includes the wave_function as .h5 for ref and conv,
# the time, norm, energy (, stepsize) for every timestep in a .jld file.
# For conv, the error and the number of B_Calls is additionally saved.
using TSSM
using JLD
include("mctdhf1d.jl")
include("check.jl")
include("../examples/time_propagators.jl")

V1(x) = 1. / 32 * (x^2)
V2(x,y) = 1/sqrt((x-y)^2+1/16)
V(x,y) = V1(x) + V1(y) + V2(x,y)

const N = 4
const interval_bound = 20
const nx = Int(512*interval_bound/10)
const t0 = 0
const groundstate_tol=1e-9

println("N = ", N, ", nx = ", nx, " on [-", interval_bound, ",", interval_bound, "]")

function method_from_string(variant::String)
    if variant == "yoshida"
        g = [1.351207191959657634, -1.702414383919315268, 1.351207191959657634] # Yoshida
        a, b = get_coeffs_composition(g)
        method = SplittingRK4BMethod(a,b)
    end
    if variant == "suzuki"
        g = [1/(4-4^(1/3)),1/(4-4^(1/3)),-4^(1/3)/(4-4^(1/3)), 1/(4-4^(1/3)), 1/(4-4^(1/3))] # Suzuki
        a, b = get_coeffs_composition(g)
        method = SplittingRK4BMethod(a,b)
    end
    if variant == "strang"
        a, b = [0.5, 0.5], [1.0, 0.0]
        method = SplittingRK4BMethod(a, b, secondorder=true)
    end
    if variant == "krogstad"
        method = ExponentialRungeKutta(:krogstad)
    end
    if variant == "rk4"
        method = ExponentialRungeKutta(:rk4)
    end
    if variant == "rklawson"
        method = ExponentialRungeKutta(:lawson)
    end
    if variant == "strehmel_weiner"
        method = ExponentialRungeKutta(:strehmel_weiner)
    end
    if variant == "adamslawson3"
        method = ExponentialMultistep(3, version=2, iters=1, starting_method=ExponentialRungeKutta(:lawson), starting_subdivision=50)
    end
    if variant == "expmultistep3"
        method = ExponentialMultistep(3, version=1, iters=1, starting_method=ExponentialRungeKutta(:lawson), starting_subdivision=50)
    end
    if variant == "adamslawson5"
        method = ExponentialMultistep(5, version=2, iters=1, starting_method=ExponentialRungeKutta(:lawson), starting_subdivision=50)
        println("adamslawson5")
    end
    if variant == "expmultistep5"
        method = ExponentialMultistep(5, version=1, iters=1, starting_method=ExponentialRungeKutta(:lawson), starting_subdivision=50)
        println("exponential multistep5")
    end
    if variant == "adamslawson3_nocorr"
        method = ExponentialMultistep(3, version=2, iters=0, starting_method=ExponentialRungeKutta(:lawson), starting_subdivision=50)
        println("adamslawson3 without correction step")
    end
    if variant == "expmultistep3_nocorr"
        method = ExponentialMultistep(3, version=1, iters=0, starting_method=ExponentialRungeKutta(:lawson), starting_subdivision=50)
        println("exponential multistep 3 without correction step")
    end
    return method
end

function ref(tol::Float64=1e-5)
    E(t) = sin(2*t)
    tend = 12

    V1_t(x, t) = E(t)*x
    V_t(x,y,t) = E(t)*(x+y)

    m = MCTDHF1D(2, N, nx, -interval_bound, interval_bound, potential1=V1, potential1_t=V1_t, potential2=V2, spin_restricted=true);
    psi = wave_function(m);
    TSSM.load!(psi, "wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_t0_1e-9.h5")
    set_time!(psi, t0)

    times, energies, norms, stepsizes = propagate_adaptive!(psi, t0, tend, tol)

    TSSM.save(psi, "results_tolerance/wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_t"*string(tend)*"_"*string(tol)*".h5")
    JLD.save("results_tolerance/wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_"*string(tol)*".jld", "times", times, "energies", energies, "norms", norms, "stepsizes", stepsizes)
end

function conv(variant::String, steps::Int)
    #laser field
    E(t) = sin(2*t)
    tend = 12
    dt=tend/steps

    V1_t(x, t) = E(t)*x
    V_t(x,y,t) = E(t)*(x+y)

    m = MCTDHF1D(2, N, nx, -interval_bound, interval_bound, potential1=V1, potential1_t=V1_t, potential2=V2, spin_restricted=true)
    psi = wave_function(m)

    TSSM.load!(psi, "wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_t0_1e-9.h5")
    set_time!(psi, t0)

    psi_ref = wave_function(m);
    TSSM.load!(psi_ref, "results_tolerance/wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_t"*string(tend)*"_1.0e-9.h5")

    method = method_from_string(variant)

    times, energies, norms = propagate_equidistant!(method, psi, t0, dt, steps)

    err = distance(psi, psi_ref)
    println("err:", err)
    TSSM.save(psi, "results_conv/wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_t"*string(tend)*"_"*variant*"_"*string(steps)*".h5")
    JLD.save("results_conv/wfquantumdot_conv_f2_n"*string(N)*"_nx"*string(nx)*"_t"*string(tend)*"_"*variant*"_"*string(steps)*".jld", "times", times, "energies", energies, "dt", dt, "err", err, "B_calls", COUNT_B)

end

function stab(variant::String, steps::Int)
    #laser field

    E(t) = 0
    tend = 70
    dt=tend/steps

    V1_t(x, t) = E(t)*x
    V_t(x,y,t) = E(t)*(x+y)

    m = MCTDHF1D(2, N, nx, -interval_bound, interval_bound, potential1=V1, potential1_t=V1_t, potential2=V2, spin_restricted=true);
    psi = wave_function(m);

    TSSM.load!(psi, "wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_t0_1e-9.h5")
    set_time!(psi, t0)

    method = method_from_string(variant)

    times, energies, norms = propagate_equidistant!(method, psi, t0, dt, steps)

    JLD.save("results_stab/wfquantumdot_stab_f2_n"*string(N)*"_nx"*string(nx)*"_t"*string(tend)*"_"*variant*"_"*string(steps)*".jld", "times", times, "energies", energies, "dt", dt, "norms", norms)
end

function propagate_equidistant!(method, psi, t0, dt, steps)
    times = Vector{Float64}()
    energies = Vector{Float64}()
    norms = Vector{Float64}()
    step_counter = 1
    time0 = time()
    for (step, tsi) in EquidistantTimeStepper(method, psi, t0, dt, steps)
        n = norm(psi)
        E_pot = potential_energy(psi)
        E_kin = kinetic_energy(psi)
        E_tot = E_pot+E_kin
        if step_counter % 50 == 0 || step == steps
            @printf("%5i  %14.10f  %14.10f  %14.10f  %14.10f  %14.10f  %14.10f %6i %10.2f\n",
                step_counter, get_time(psi), E_pot, E_kin, E_tot, dt, n, COUNT_B, time()-time0)
        end
        append!(times, get_time(psi))
        append!(energies, E_tot)
        append!(norms, n)
        step_counter += 1
    end
    return times, energies, norms
end

function propagate_adaptive!(psi, t0, tend, tol)
    method = AdaptiveAdamsLawson(6)
    times = Vector{Float64}()
    energies = Vector{Float64}()
    stepsizes = Vector{Float64}()
    norms = Vector{Float64}()
    step_counter = 1
    time0 = time()
    told = t0
    for (t, tsi) in AdaptiveTimeStepper(method, psi, t0, tend, tol, 0.00001)
        n = norm(psi)
        E_pot = potential_energy(psi)
        E_kin = kinetic_energy(psi)
        E_tot = E_pot+E_kin
        stepsize = t - told
        if step_counter % 50 == 0
            @printf("%5i  %14.10f  %14.10f  %14.10f  %14.10f  %14.10f  %14.10f  %10.2f\n",
                step_counter, get_time(psi), n, E_pot, E_kin, E_tot, stepsize, time()-time0)
        end
        append!(times, get_time(psi))
        append!(energies, E_tot)
        append!(stepsizes, stepsize)
        append!(norms, n)
        told = t
        step_counter += 1
    end
    return times, energies, norms, stepsizes
end

if ARGS[1] == "groundstate"
    println("calculate groundstate with tolerance ", groundstate_tol)
    include("propagators.jl")
    m = MCTDHF1D(2, N, nx, -interval_bound, interval_bound, potential1=V1, potential2=V2, spin_restricted=true)
    psi = wave_function(m)
    groundstate!(psi, dt=0.1, max_iter=500000, output_step=100, tol=groundstate_tol)
    TSSM.save(psi, "wfquantumdot_f2_n"*string(N)*"_nx"*string(nx)*"_t0_1e-9.h5")
elseif ARGS[1] == "ref"
    println("calculate reference solution with tolerance ", ARGS[2])
    ref(parse(Float64, ARGS[2]))
elseif ARGS[1] == "conv"
    println("calculate solution using method ", ARGS[2], " with ", ARGS[3], " steps")
    conv(ARGS[2], parse(Int, ARGS[3]))
elseif ARGS[1] == "stab"
    stab(ARGS[2], parse(Int, ARGS[3]))
else
    throw(ErrorException("first cmdline argument is not 'ref', 'conv' or 'stab'"))
end