text:
include("../src/RationalFunctionFields.jl")
using .RationalFunctionFields: change_parent_ring, singular2aa, aa2singular,
                        double_singular2aa, double_aa2singular

@testset "Utils tests" begin

    R, xs = Sing.PolynomialRing(Sing.QQ, ["x", "y"])
    x, y = xs
    f = x + y

    faa = singular2aa(f)

    revfaa = aa2singular(faa, new_ring=R)
    @test revfaa == f

    R2, xs2 = Sing.PolynomialRing(R, ["x2", "y2"])
    x2, y2 = xs2
    f2 = x2 + y2

    f2aa = double_singular2aa(f2)

    revf2s = double_aa2singular(f2aa, base=R, new_ring=R2)
    @test revf2s == f2

    AA_R, (x, y) = AA.PolynomialRing(Sing.QQ, ["x", "y"])
    f = x + y
    newparent, (xx, yy, tt) = AA.PolynomialRing(Sing.QQ, ["x", "y", "t"])

    ft, t = change_parent_ring(f, newparent)

    @test parent(ft) == newparent
    @test t == tt

    frev = change_parent_ring(ft, AA_R)

    @test f == frev

end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   BondHoppingHamiltonianSimple <: AbstractBondHamiltonian{L,1}
#   --> L is the label type of bonds
#   --> N is the dimension the bond term matrix (NxN matrix)
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - generator function
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
mutable struct BondHoppingHamiltonianSimple{L} <: AbstractBondHamiltonian{L,1}

    # coupling
    coupling :: Float64

end

# export the type
export BondHoppingHamiltonianSimple

################################################################################
#
#	INTERFACING / ACCESSING BOND HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# get bond term of Hamiltonian
function bondterm(
            h :: BondHoppingHamiltonianSimple{L},
            b :: AbstractBond{L,NB}
        ) :: Float64 where {N,L,NB,NS}

    # every bond has the same coupling
    return h.coupling
end

# export the function
export bondterm

function getHoppingHamiltonianSimple(
            unitcell :: U
        ) :: BondHoppingHamiltonianSimple{L} where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B}}

    # just create a new hamiltonian
    return BondHoppingHamiltonianSimple{L}(1.0)
end

export getHoppingHamiltonianSimple

function getBondHamiltonianType(
        ::Val{:BondHoppingHamiltonianSimple}
    )

    # return the type
    return BondHoppingHamiltonianSimple
end

function saveBondHamiltonian(
        hb :: HB,
        fn :: AbstractString,
        group :: AbstractString = "bond_hamiltonian"
        ;
        append :: Bool = false
    ) where {L,HB<:BondHoppingHamiltonianSimple{L}}

    # determine the mode based on if one wants to append stuff
    if append
        mode = "r+"
    else
        mode = "w"
    end

    # open the file in mode
    h5open(fn, mode) do file
        # create the group in which the bonds are saved
        group_hb = g_create(file, group)
        # save the type identifier
        attrs(group_hb)["type"] = "BondHoppingHamiltonianSimple"
        # save the parameters
        attrs(group_hb)["N"] = 1
        attrs(group_hb)["L"] = string(L)
        # save the Float64 coupling
        group_hb["coupling"] = hb.coupling
    end

    # return nothing
    return nothing
end

function loadBondHamiltonian(
        ::Type{HB},
        fn :: AbstractString,
        group :: AbstractString = "bond_hamiltonian"
    ) where {LI,HB<:Union{BondHoppingHamiltonianSimple{LI},BondHoppingHamiltonianSimple}}

    # read attribute data
    attr_data = h5readattr(fn, group)
    # determine D based on this
    L = Meta.eval(Meta.parse(attr_data["L"]))
    N = attr_data["N"]

    # load all remaining data
    coupling = h5read(fn, group*"/coupling")

    # return the new bond hamiltonian
    return BondHoppingHamiltonianSimple{L}(coupling)
end

ext: jl


text:
using UniqueVectors: UniqueVector

export Blackbox
export impedance_matrices, admittance_matrices, scattering_matrices

struct Blackbox{T, W<:Number, Y<:AbstractMatrix{<:Number},
        P<:AbstractMatrix{<:Number}} <: AdmittanceModel{T}
    ω::Vector{W}
    Y::Vector{Y}
    P::P
    ports::UniqueVector{T}
    function Blackbox{T}(ω::AbstractVector, y::AbstractVector, p::AbstractMatrix,
            ports::UniqueVector{T}) where {T}
        @assert all(isreal(m) for m in (ω, p))
        @assert length(ω) == length(y)
        if !isempty(y)
            l = size(first(y), 1)
            @assert size(p) == (l, length(ports))
            @assert all(size(yy) == (l, l) for yy in y)
        end
        ww, pp = float.((ω, p))
        yy = float.(y)
        W = eltype(ww)
        Y = eltype(yy)
        P = typeof(pp)
        return new{T, W, Y, P}(ww, yy, pp, ports)
    end
end

"""
    Blackbox(ω::AbstractVector{<:Number}, Y::AbstractVector,
        P::AbstractMatrix, ports::AbstractVector{T}) where {T}
Create a `Blackbox`, which has an admittance matrix `Y(ω′)` sampled over a
discrete set of frequencies `ω` such that `Y[i] == Y(ω[i])`, for valid `i`.
The port names are of type `T`.
"""
function Blackbox(ω::AbstractVector{<:Number}, Y::AbstractVector,
        P::AbstractMatrix, ports::AbstractVector{T}) where {T}
    return Blackbox{T}(ω, Y, P, UniqueVector(ports))
end

get_Y(bbox::Blackbox) = bbox.Y

get_P(bbox::Blackbox) = bbox.P

get_ports(bbox::Blackbox) = bbox.ports

function partial_copy(bbox::Blackbox{T, U};
        Y::Union{Vector{V}, Nothing} = nothing,
        P::Union{V, Nothing} = nothing,
        ports::Union{AbstractVector{W}, Nothing} = nothing) where {T, U, V, W}
    Y = isnothing(Y) ? get_Y(bbox) : Y
    P = isnothing(P) ? get_P(bbox) : P
    ports = isnothing(ports) ? get_ports(bbox) : ports
    return Blackbox(bbox.ω, Y, P, ports)
end

function compatible(bboxes::AbstractVector{<:Blackbox})
    if isempty(bboxes)
        return true
    end
    ω = bboxes[1].ω
    return all(ω == bbox.ω for bbox in bboxes[2:end])
end

# TODO for sparse matrices use sparse linear solver instead of \ and /

function impedance_to_scattering(Z::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    U = diagm(0 => Z0)
    return (Z + U) \ (Z - U)
end

function scattering_to_impedance(S::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    U = diagm(0 => Z0)
    return U * (I + S) / (I - S)
end

function admittance_to_scattering(Y::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    U = diagm(0 => Z0)
    return (I + Y * U) \ (I - Y * U)
end

function scattering_to_admittance(S::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    V = diagm(0 => 1 ./ Z0)
    return (I - S) * ((I + S) \ V)
end

"""
    impedance_matrices(bbox::Blackbox)

Find the impedance matrices `Z` with `y = Zx` for each `ω`.
"""
impedance_matrices(bbox::Blackbox) = [transpose(bbox.P) * (Y \ collect(bbox.P))
    for Y in bbox.Y]

"""
    admittance_matrices(bbox::Blackbox)

Find the admittance matrices `Y` with `Yy = x` for each `ω`.
"""
function admittance_matrices(bbox::Blackbox)
    if bbox.P == I # avoid unnecessary matrix inversion
        return map(collect, bbox.Y)
    else
        return inv.(impedance_matrices(bbox))
    end
end

"""
    scattering_matrices(bbox::Blackbox, Z0::AbstractVector{<:Real})

Find the scattering matrices that map incoming waves to outgoing waves on transmission
lines with characteristic impedance `Z0` for each `ω`.
"""
function scattering_matrices(bbox::Blackbox, Z0::AbstractVector{<:Number})
    @assert isreal(Z0)
    if bbox.P == I # avoid unnecessary matrix inversion
        return [admittance_to_scattering(Y, Z0) for Y in admittance_matrices(bbox)]
    else
        return [impedance_to_scattering(Z, Z0) for Z in impedance_matrices(bbox)]
    end
end

function canonical_gauge(bbox::Blackbox)
    n = length(bbox.ports)
    Y = admittance_matrices(bbox)
    P = Matrix{eltype(eltype(Y))}(I, n, n)
    return Blackbox(bbox.ω, Y, P, bbox.ports)
end

"""
    Blackbox(ω::AbstractVector{<:Number}, pso::PSOModel)

Create a Blackbox for a PSOModel where `ω` represents angular frequency.
"""
function Blackbox(ω::AbstractVector{<:Number}, pso::PSOModel)
    @assert isreal(ω)
    K, G, C = get_Y(pso)
    Y = [K ./ s + G + C .* s for s in im * ω]
    P = complex(float(get_P(pso)))
    return Blackbox(ω, Y, P, get_ports(pso))
end


text:
@info "Hello! It's sim time!"
include("helperFunctions.jl")
outputFileBaseName, io = setupDirectories(directoryname)
@info "Imported helperFunctions.jl"
flush(io)

# Setup multi threading
# threadsAvailable = length(Sys.cpu_info())
# if numThreads > threadsAvailable
    # Threads.nthreads() = threadsAvailable
    # @warn "Assigned more threads than available."
# elseif numThreads < threadsAvailable
    Threads.nthreads() = numThreads
# end
@info "using $(Threads.nthreads()) CPU threads"# of $threadsAvailable available threads."
flush(io)

## Initialize particles
h0, f0, η, ε, resolution = generateFlatParticleDistribution(numParticles, ICrange, z0, λ0);
numParticles = length(h0[:,1]);      # new total number of particles may have changed.
@info "All prepped, running model on $numParticles particles."
flush(io)
# hf0 = [h0  f0][shuffle(1:end), :];   # shuffle so batches all take same-ish time

## Setup model
params = @SVector [η, ε, Omegape, omegam, a, dPhi, dλ1, dλ2, B_w_normalizer];
@info "Using these parameters: $params"
flush(io)
prob = ODEProblem(eom!, ~, tspan, params);
# _prob, _alg = auto_optimize(prob) TODO: get auto-optimize to work

@everywhere probGeneratorList, nPerBatch, percentage = generateModifiableFunction(batches); # consider making all of these constant
flush(io)

## Run simulation
tick()
flush(io)
ensemble()
@info "Total sim time:"
tock()
flush(io)

max_stars_repo_name: ethantsai/nlwhistlers

max_issues_repo_name: ethantsai/nlwhistlers

max_forks_repo_name: ethantsai/nlwhistlers


text:
# -----------------------------------------------
# Minimal Distributions routines
# -----------------------------------------------

# TEST:
# using Rmath

#=
macro QD2(A, B) quote
    esc(:(sum((A - B).^2)))
end
=#

abstract type AbstractDistributionType end

# #########   NORMAL DISTRIBUTION   #########

abstract type Normal<:AbstractDistributionType end

function PDF{T<:Number}(x::T, dtype::Type{Normal})
    return dnorm(x, 0, 1)
end

function PDF{Tx<:Number, Tmu<:Number, Tsig<:Number}(x::Tx,  mu::Tmu, sigma²::Tsig, dtype::Type{Normal})
    return dnorm(x, mu, sqrt(sigma²))
end

function PDF{Tx<:Number, Tmu<:Number, Tsig<:Number}(X::Array{Tx,1}, mu::Tmu, sigma²::Tsig, dtype::Type{Normal})
    # return exp.(-(X - mu).^2 / (2*sigma²)) / sqrt(2*pi*sigma²);
    return [ dnorm(x, mu, sqrt(sigma²)) for x in X ]
end

function PDF{T<:Number}(X::Array{T}, mu::Array{T,1}, Σ::Array{T},  dtype::Type{Normal})
    m, n = size(X)
    if length(dtype.mu) == n && size(dtype.Σ)[1] == n
        return sqrt(2*pi*det(Σ)) * exp.(-0.5 * (X' - row(mu)') * inv(Σ) * (X - row(mu)));
    else
        error("X, mu and Σ have not valid dimensions")
    end
end

function CDF{T<:Number}(x::T, dtype::Type{Normal})
    return pnorm(x, 0, 1)
end

function CDF{Tx<:Number, Tmu<:Number, Tsig<:Number}(x::Tx,  mu::Tmu, sigma²::Tsig, dtype::Type{Normal})
    # return CDF((sqrt(sigma²)*x + mu), Normal)
    return pnorm(x, mu, sqrt(sigma²))
end

function CDF{Tx<:Number, Tmu<:Number, Tsig<:Number}(X::Array{Tx,1}, mu::Tmu, sigma²::Tsig, dtype::Type{Normal})
    # return exp.(-(X - mu).^2 / (2*sigma²)) / sqrt(2*pi*sigma²);
    return [ pnorm(x, mu, sqrt(sigma²)) for x in X ]
end

#=
function CDF{T<:Number}(X::Array{T}, mu::Array{T,1}, Σ::Array{T},  dtype::Normal)
    return 0;
end
=#

function logLike{T<:Number}(X::Array{T,1}, mu::T, sigma²::T, dtype::Type{Normal})
    N = length(X)
    return - N/2*log(sigma²) - (0.5/sigma²)*sum((X - mu).^2)
end

function logLike{T<:Number}(X::Array{T,1}, mu::Array{T,1}, Σ::Array{T}, dtype::Type{Normal})
    N = length(X)
    return - N/2*log(det(Σ)) - 0.5 * sum(X' - row(mu)') * inv(Σ) * (X - row(mu));
end

# println(@time CDF([-1.0, 1], 1.0, 4.0, Normal))

# #########   STUDENT DISTRIBUTION   #########

abstract type Student<:AbstractDistributionType end

function PDF{Tx<:Number, Tnu<:Number}(x::Tx, ν::Tnu, dtype::Type{Student})
    # dt(x, ν) == gamma((ν+1)/2) / (sqrt(ν*pi) * gamma(ν/2)) * (1 + x.^2/ν).^(-(ν+1)/2)
    return dt(x, ν)
end

function PDF{Tx<:Number, Tnu<:Number}(X::Array{Tx,1}, ν::Tnu, dtype::Type{Student})
    return [ dt(x, ν) for x in X ]
end

function CDF{Tx<:Number, Tnu<:Number}(x::Tx, ν::Tnu, dtype::Type{Student})
    return pt(x, ν)
end

function CDF{Tx<:Number, Tnu<:Number}(X::Array{Tx,1}, ν::Tnu, dtype::Type{Student})
    return [ pt(x, ν) for x in X ]
end

# println(@time CDF(1.0, 1.0, Student))

# #########   LOGISTIC   #########

abstract type Logistic<:AbstractDistributionType end

logistic{T<:Number}(X::Array{T,1}) = 1 ./ (1 + exp.(X));

function logLike{T<:Number}(X::Array{T}, Y::Array{T}, theta::Array{T}, dtype::Type{Logistic})
    N = length(Y)
    return sum((logistic(X * theta)).^Y .* (1 - logistic(X * theta)).^(1 - Y)) / N
end


text:
#=
$title Economic Load Dispatch

$onText
For more details please refer to Chapter 3 (Gcode3.1), of the following book:
Soroudi, Alireza. Power System Optimization Modeling in GAMS. Springer, 2017.
--------------------------------------------------------------------------------
Model type: QCP
--------------------------------------------------------------------------------
Contributed by
Dr. Alireza Soroudi
IEEE Senior Member
Email: alireza.soroudi@gmail.com
We do request that publications derived from the use of the developed GAMS code
explicitly acknowledge that fact by citing
Soroudi, Alireza. Power System Optimization Modeling in GAMS. Springer, 2017.
DOI: doi.org/10.1007/978-3-319-62350-4
$offText

Converstion from GAMS to JuMP by Jon Duan
=#

N_Gen = 5
Gen_set = 1:N_Gen
Gen = [Symbol("G$i") for i in Gen_set]

Other_set = ["a"     ,"b",      "c",       "Pmin",  "Pmax"]
Other = [Symbol("$i") for i in Other_set]

Load = 400

# https://github.com/davidavdav/NamedArrays.jl
# using DataStructures
using NamedArrays
data = Float64[  [3     20     100     28    206];
                 [4.05  18.07  98.87   90    284];
                 [4.05  15.55  104.26  68    189];
                 [3.99  19.21  107.21  76    266];
                 [3.88  26.18  95.31   19    53]]

data =  NamedArray(data, ( Gen, Other ), ("Rows", "Cols"))
# df[:y2] = parse.(Int,df[:yearsAsString])

using JuMP
using Gurobi

ED = Model(solver=GurobiSolver())

@variable(ED,  data[i, :Pmin]<= P[i in Gen]<=data[i, :Pmax])

@variable(ED, OF>=0)

@objective(ED, Min, sum( data[i, :a]*P[i]^2 + data[i,:b]*P[i] for i in Gen))

@constraint(ED, Balance, sum( P[i]for i in Gen ) >= Load)

# Solve the problem
solve(ED)

# Print our final solution
println("Objective value: ", getobjectivevalue(ED))
println("P = ", getvalue(P))
println("Marginal Value of Balance = ", getdual(Balance))
# check balance
sum(getvalue(P))
# check margian value for
getdual(P)

# Pkg.add("ExcelFiles")
# Info: Using the Python distribution in the Conda package by default.
# To use a different Python version, set ENV["PYTHON"]="pythoncommand" and re-run Pkg.build("PyCall").

using ExcelFiles, DataFrames

years = [Symbol("$i") for i in 2006:2015];

demand = DataFrame(load("Data.xlsx", "Load!A1:K8761"));

head(demand)

demand = demand[2:end];

demand.colindex

size(demand)

windSouthEast = DataFrame(load("Data.xlsx", "SouthEast!A1:K8761"));

head(windSouthEast)

windNorth = DataFrame(load("Data.xlsx", "North!A1:K8761"));

windNorth =windNorth[2:end];

windSouthWest = DataFrame(load("Data.xlsx", "SouthWest!A1:K8761"));

windSouthWest = windSouthWest[2:end];

#SouthsolarInput rng=PINCHER-CREEK
solarSouth = DataFrame(load("Data.xlsx", "PINCHER-CREEK"));

solarSouth = solarSouth[2:end-1];

# par= MiddleSolarInput rng=EDMONTON-STONY-PLAIN-CS!A1:K8761 Rdim=1 Cdim=1
# par= NorthSolarInput rng=GRANDE-PRAIRIE-A!A1:K8761 Rdim=1 Cdim=1
solarMiddle = DataFrame(load("Data.xlsx", "EDMONTON-STONY-PLAIN-CS!A1:K8761"));

#solarMiddle

solarMiddle = solarMiddle[2:end];

solarNorth = DataFrame(load("Data.xlsx", "GRANDE-PRAIRIE-A!A1:K8761"));

solarNorth = solarNorth[2:end];

#solarNorth[Symbol("2006")]

size(solarSouth)

#solarSouth

#solarSouth

typeof(solarSouth)

#solarSouth1 = similar(solarSouth)

#* the unit of output in data table is w. transfor to MW. A solar farm with 50000 PV panel? may cause floating point problem, number is too small

for (i,j) in eachcol(solarSouth) # a tuple with symbol and array
  solarSouth[i] = j*50000/1000000;
end

for (i,j) in eachcol(solarMiddle) # a tuple with symbol and array
  solarMiddle[i] = j*50000/1000000;
end

for (i,j) in eachcol(solarNorth) # a tuple with symbol and array
  solarNorth[i] = j*50000/1000000;
end

#Para.csv
using CSV
para =  CSV.read("Para.csv")

para

para[1].== "wind"

turbine = parse(para[ para[1].== "wind",:UnitCap][1]);

turbine

panel =  parse(para[ para[1].== "solar",:UnitCap][1]);
panel

generation_assets = [Symbol(i) for i in para[:,1][2:end]]

fossil_fuel = [Symbol(i) for i in para[:,1][end-1:end]]

renewables = [Symbol(i) for i in para[:,1][2:3]]

windLoc = [:se, :sw, :n]
solarLoc = [:s, :m, :n]

size(para)

#para.columns[1]

para.colindex

FixOMCosts = Dict( zip( generation_assets    , ([parse(i) for i in para[:FixOMCost][2:end]]) )   )

VarOMCosts = Dict( zip( generation_assets    , ([parse(i) for i in para[:VarOMCost][2:end]]) )   )

ONights= Dict( zip( generation_assets    , ([parse(i) for i in para[:ONight][2:end]]) )   )

UnitCaps  = Dict( zip( generation_assets    , ([parse(i) for i in para[:UnitCap][2:end]]) )   )

Maxis  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Max][2:end]]) )   )

Emits  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Emit][2:end]]) )   )

Ramps  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Ramp][2:end]]) )   )

LTimes= Dict( zip( generation_assets    , ([parse(i) for i in para[:LTime][2:end]]) )   )

# # parameters
# ''''''
# round = Para("battery", "Round") ;
# duration = Para("battery", "duration") ;
# taxEm = Para("CT","CTax");
# turbine = Para("wind","UnitCap");
# maxturbine = Para("wind","Max");
# panel = Para("solar","UnitCap");
# maxPanel = Para("solar","Max");
# rate = Para("CT","Discount");
# ''''''

maximums = Dict( zip( generation_assets    , ([parse(i) for i in para[:Max][2:end]]) )   )

roundEfficiency = Dict( zip( generation_assets    , ([parse(i) for i in para[:Round][2:end]]) )   )

duration= Dict( zip( generation_assets    , ([i for i in para[:duration][2:end]]) )   )

taxEm = Dict( zip( generation_assets    , ([parse(i) for i in para[:Ctax][2:end]]) )   )

rate = Dict( zip( generation_assets    , ([parse(i) for i in para[:Discount][2:end]]) )   )

Vcost = Dict( zip( generation_assets    , ( VarOMCosts[i] + taxEm[i]* Emits[i]  for i in generation_assets) )   ) ;

Vcost

Fcost =   Dict( zip(generation_assets, (ONights[i]*rate[i]*((1+rate[i])^LTimes[i])/(((1+rate[i])^LTimes[i])-1) +
                                            FixOMCosts[i] for i in generation_assets)
                    )
                );

Fcost

using JuMP
using Gurobi

m = Model(solver=GurobiSolver())

T = 8760

@variable(m, batteryCap>=0)

# energy in battery
@variable(m, v[1:T]>=0)

@variable(m, charge[1:T]>=0)
@variable(m, discharge[1:T]>=0)

@variable(m, fuelCap[fossil_fuel] >=0)

@variable(m, gen[ 1:T,fossil_fuel] >=0)

@variable(m, nSolarPlant[solarLoc] >=0)
@variable(m, nWindTurbine[windLoc] >=0)

wind = DataFrame(s = 1.0:T,m = 1.0:T, n = 1.0:T);

solar = DataFrame(s = 1.0:T,m = 1.0:T, n = 1.0:T);

year = 2006

Symbol("$year")

wind[:se] = windSouthEast[Symbol("$year")];
wind[:sw] = windSouthWest[Symbol("$year")];
wind[:n] = windNorth[ Symbol("$year")];

#* solar unit output per hour  (MW), number of panel is in 50000, so total solar ouput is MW
solar[:s] = solarSouth[Symbol("$year")];
solar[:m] = solarMiddle[Symbol("$year")];
solar[:n] = solarNorth[Symbol("$year")];

demand0 = demand[Symbol("$year")];

@objective(m, Min, sum(sum(wind[t,L]*nWindTurbine[L]*Vcost[:wind] for L in windLoc) +
                            (sum(solar[t,Ls]*nSolarPlant[Ls]*Vcost[:solar] for Ls in solarLoc))  +
                            (sum(gen[t,f] *Vcost[f] for f in fossil_fuel )) +
                            discharge[t]*Vcost[:battery]
                            for t in 1:T
                            ) +
                    sum(Fcost[:wind]*nWindTurbine[L]*turbine for L in windLoc) +
                    sum(Fcost[:solar]*nSolarPlant[Ls]*turbine for Ls in solarLoc) +
                    sum(Fcost[f]*fuelCap[f] for f in fossil_fuel) +
                    Fcost[:battery]*batteryCap
            );

@constraint(m,  sum(wind[1,L]*nWindTurbine[L] for L in windLoc) +
                                sum(solar[1,Ls]*nSolarPlant[Ls] for Ls in solarLoc) +
                                  sum(gen[1,f] for f in fossil_fuel) +
                                   discharge[1] -charge[1] >= demand0[1]);

# for t in 1:T
#     @constraint(m, Balance[t], sum(wind[t,L].*nWindTurbine[L] for L in windLoc) +
#                                 sum(solar[t,Ls].*nSolarPlant[Ls] for Ls in solarLoc) +
#                                   sum(gen[t,f] for f in fossil_fuel) +
#                                    discharge[t] -charge[t] >= demand0[t]);
# end

@constraint(m, Balance0[t=1:T], sum(wind[t,L].*nWindTurbine[L] for L in windLoc) +
                            sum(solar[t,Ls].*nSolarPlant[Ls] for Ls in solarLoc) +
                              sum(gen[t,f] for f in fossil_fuel) +
                               discharge[t] -charge[t] >= demand0[t]);

 @constraint(m, capacityFuel[t=1:T, f = fossil_fuel], gen[t,f] <= fuelCap[f] );

 @constraint(m, capacityWind[L=windLoc], nWindTurbine[L] <= maximums[:wind] );

 @constraint(m, capacitySolar[Ls=solarLoc], nSolarPlant[Ls] <= maximums[:solar] );

@constraint(m, RampUp[t=2:T, f = fossil_fuel], gen[t, f] <= gen[t-1,f]+Ramps[f]*fuelCap[f] );

@constraint(m, RampDn[t=2:T, f = fossil_fuel], gen[t, f] >= gen[t-1,f]-Ramps[f]*fuelCap[f] );

@constraint(m, Storage[t=2:T], v[t] ==   v[t-1]+(roundEfficiency[:battery]*charge[t-1]-discharge[t-1]));

@constraint(m, MaxStorage[t=1:T], v[t] <= duration[:battery] * batteryCap);

@constraint(m, MaxCharge[t=1:T], charge[t] <= duration[:battery] * batteryCap);

@constraint(m, MaxDischarge[t=1:T], discharge[t] <= batteryCap);

solve(m)

# # include("examples_julia/2_uc.jl");

# n_machines = 3;
# time_steps = 12;

# fixed_costs = [5, 10, 1000];
# variable_costs = [800, 700, 20];
# max_power = [500, 700, 1500];
# startup_costs = [10, 5, 850];
# demands = [1500, 1200, 1800, 2000, 2200, 1750, 1600, 1500, 1200, 1500, 1600, 1800];

# using JuMP
# using Gurobi

# m = Model(solver=GurobiSolver())

# @variable(m, on[1:time_steps, 1:n_machines], Bin)
# @variable(m, power[1:time_steps, 1:n_machines] >= 0)
# @variable(m, start[1:time_steps, 1:n_machines], Bin)

# @objective(m, Min, sum(dot(fixed_costs, vec(on[t, :])) +
#         dot(variable_costs, vec(power[t, :])) + dot(startup_costs, vec(start[t, :])) for t in 1:time_steps))

# for t in 1:time_steps
#     for i in 1:n_machines
#         @constraint(m, power[t, i] <= max_power[i] * on[t, i])
#     end
#     @constraint(m, sum(power[t, :]) == demands[t])
# end

# for i in 1:n_machines
#     @constraint(m, start[1, i] == on[1, i])
# end

# for t in 2:time_steps
#     for i in 1:n_machines
#         @constraint(m, start[t, i] >= on[t, i] - on[t - 1, i])
#     end
# end

# for t in 1:time_steps
#     for i in 1:n_machines
#         @constraint(m, on[t, i] >= start[t, i])
#         if t <= time_steps - 1
#             @constraint(m, on[t + 1, i] >= start[t, i])
#         end
#         if t <= time_steps - 2
#             @constraint(m, on[t + 2, i] >= start[t, i])
#         end
#     end
# end

# print(m)
# solve(m)
# getvalue(on)

Gurobi.Env()


text:
module RandUtils

using Compat

if VERSION >= v"0.7.0-DEV.3406"
    using Random
end
if VERSION >= v"0.7.0-DEV.3389"
    using SparseArrays
end

const SEED = MersenneTwister(1)
const DEFAULT_ELTYPE = typeof(rand())

samerand(args...) = rand(deepcopy(SEED), args...)
samesprand(args...) = sprand(deepcopy(SEED), args...)

samesprandbool(args...) = sprand(deepcopy(SEED), Bool, args...)

randvec(T, n) = samerand(T, n)
randvec(n) = samerand(n)

randmat(T, n) = samerand(T, n, n)
randmat(n) = samerand(n, n)

export samerand, samesprand, samesprandbool, randvec, randmat

end


text:
using Documenter, SpeedDate, Pkg

makedocs(
    modules = [SpeedDate],
    format = Documenter.HTML(),
    sitename = "SpeedDate.jl",
    authors = replace(join(Pkg.TOML.parsefile("Project.toml")["authors"], ", "), r" <.*?>" => "" ),
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Introduction" => "man/intro.md",
            "Loading Data" => "man/loading.md",
            "Dating Sequences" => "man/dating.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/vanOosterhoutLab/SpeedDate.jl.git",
    push_preview = true,
    deps = nothing,
    make = nothing
)


text:
# This file includes code that was formerly a part of Julia. License is MIT: http://julialang.org/license

## basic UTF-8 decoding & iteration ##

const utf8_offset = [
    0x00000000, 0x00003080,
    0x000e2080, 0x03c82080,
    0xfa082080, 0x82082080,
]

const utf8_trailing = [
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,
]

## required core functionality ##

function lastindex(s::UTF8String)
    d = s.data
    i = length(d)
    i == 0 && return i
    while is_valid_continuation(d[i])
        i -= 1
    end
    i
end

codeunit(s::UTF8String) = UInt8
ncodeunits(s::UTF8String) = length(s.data)

function length(s::UTF8String)
    d = s.data
    cnum = 0
    for i = 1:length(d)
        @inbounds cnum += !is_valid_continuation(d[i])
    end
    cnum
end

function next(s::UTF8String, i::Int)
    # potentially faster version
    # d = s.data
    # a::UInt32 = d[i]
    # if a < 0x80; return Char(a); end
    # #if a&0xc0==0x80; return '\ufffd'; end
    # b::UInt32 = a<<6 + d[i+1]
    # if a < 0xe0; return Char(b - 0x00003080); end
    # c::UInt32 = b<<6 + d[i+2]
    # if a < 0xf0; return Char(c - 0x000e2080); end
    # return Char(c<<6 + d[i+3] - 0x03c82080)

    d = s.data
    b = d[i]
    if is_valid_continuation(b)
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, d[i]))
    end
    trailing = utf8_trailing[b+1]
    if length(d) < i + trailing
        return '\ufffd', i+1
    end
    c::UInt32 = 0
    for j = 1:trailing+1
        c <<= 6
        c += d[i]
        i += 1
    end
    c -= utf8_offset[trailing+1]
    Char(c), i
end

if isdefined(Base, :iterate)
    function iterate(s::UTF8String, i::Int = firstindex(s))
        i > ncodeunits(s) && return nothing
        return next(s, i)
    end
end

function first_utf8_byte(ch::Char)
    c = UInt32(ch)
    c < 0x80    ? c%UInt8 :
    c < 0x800   ? ((c>>6)  | 0xc0)%UInt8 :
    c < 0x10000 ? ((c>>12) | 0xe0)%UInt8 :
                  ((c>>18) | 0xf0)%UInt8
end

function reverseind(s::UTF8String, i::Integer)
    j = lastidx(s) + 1 - i
    d = s.data
    while is_valid_continuation(d[j])
        j -= 1
    end
    return j
end

## overload methods for efficiency ##

bytestring(s::UTF8String) = s

sizeof(s::UTF8String) = sizeof(s.data)

lastidx(s::UTF8String) = length(s.data)

isvalid(s::UTF8String, i::Integer) =
    (1 <= i <= lastindex(s.data)) && !is_valid_continuation(s.data[i])

const empty_utf8 = UTF8String(UInt8[])

function getindex(s::UTF8String, r::UnitRange{Int})
    isempty(r) && return empty_utf8
    i, j = first(r), last(r)
    d = s.data
    if i < 1 || i > length(s.data)
        throw(BoundsError(s, i))
    end
    if is_valid_continuation(d[i])
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, d[i]))
    end
    if j > length(d)
        throw(BoundsError())
    end
    j = nextind(s,j)-1
    UTF8String(d[i:j])
end

function search(s::UTF8String, c::Char, i::Integer)
    if i < 1 || i > sizeof(s)
        i == sizeof(s) + 1 && return 0
        throw(BoundsError(s, i))
    end
    d = s.data
    if is_valid_continuation(d[i])
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, d[i]))
    end
    c < Char(0x80) && return search(d, c%UInt8, i)
    while true
        i = search(d, first_utf8_byte(c), i)
        (i==0 || s[i] == c) && return i
        i = next(s,i)[2]
    end
end

function rsearch(s::UTF8String, c::Char, i::Integer)
    c < Char(0x80) && return rsearch(s.data, c%UInt8, i)
    b = first_utf8_byte(c)
    while true
        i = rsearch(s.data, b, i)
        (i==0 || s[i] == c) && return i
        i = prevind(s,i)
    end
end

function string(a::ByteString...)
    if length(a) == 1
        return a[1]::UTF8String
    end
    # ^^ at least one must be UTF-8 or the ASCII-only method would get called
    data = Vector{UInt8}(undef, 0)
    for d in a
        append!(data,d.data)
    end
    UTF8String(data)
end

function reverse(s::UTF8String)
    dat = s.data
    n = length(dat)
    n <= 1 && return s
    buf = Vector{UInt8}(undef, n)
    out = n
    pos = 1
    @inbounds while out > 0
        ch = dat[pos]
        if ch > 0xdf
            if ch < 0xf0
                (out -= 3) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                buf[out + 1], buf[out + 2], buf[out + 3] = ch, dat[pos + 1], dat[pos + 2]
                pos += 3
            else
                (out -= 4) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                buf[out+1], buf[out+2], buf[out+3], buf[out+4] = ch, dat[pos+1], dat[pos+2], dat[pos+3]
                pos += 4
            end
        elseif ch > 0x7f
            (out -= 2) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
            buf[out + 1], buf[out + 2] = ch, dat[pos + 1]
            pos += 2
        else
            buf[out] = ch
            out -= 1
            pos += 1
        end
    end
    UTF8String(buf)
end

## outputting UTF-8 strings ##

write(io::IO, s::UTF8String) = write(io, s.data)

## transcoding to UTF-8 ##

utf8(x) = convert(UTF8String, x)
convert(::Type{UTF8String}, s::UTF8String) = s
convert(::Type{UTF8String}, s::ASCIIString) = UTF8String(s.data)
convert(::Type{SubString{UTF8String}}, s::SubString{ASCIIString}) =
    SubString(utf8(s.string), s.offset+1, ncodeunits(s)+s.offset)

function convert(::Type{UTF8String}, dat::Vector{UInt8})
    # handle zero length string quickly
    isempty(dat) && return empty_utf8
    # get number of bytes to allocate
    len, flags, num4byte, num3byte, num2byte = unsafe_checkstring(dat)
    if (flags & (UTF_LONG | UTF_SURROGATE)) == 0
        len = sizeof(dat)
        @inbounds return UTF8String(copyto!(Vector{UInt8}(undef, len), 1, dat, 1, len))
    end
    # Copy, but eliminate over-long encodings and surrogate pairs
    len += num2byte + num3byte*2 + num4byte*3
    buf = Vector{UInt8}(undef, len)
    out = 0
    pos = 0
    @inbounds while out < len
        ch::UInt32 = dat[pos += 1]
        # Handle ASCII characters
        if ch <= 0x7f
            buf[out += 1] = ch
        # Handle overlong < 0x100
        elseif ch < 0xc2
            buf[out += 1] = ((ch & 3) << 6) | (dat[pos += 1] & 0x3f)
        # Handle 0x100-0x7ff
        elseif ch < 0xe0
            buf[out += 1] = ch
            buf[out += 1] = dat[pos += 1]
        elseif ch != 0xed
            buf[out += 1] = ch
            buf[out += 1] = dat[pos += 1]
            buf[out += 1] = dat[pos += 1]
            # Copy 4-byte encoded value
            ch >= 0xf0 && (buf[out += 1] = dat[pos += 1])
        # Handle surrogate pairs
        else
            ch = dat[pos += 1]
            if ch < 0xa0 # not surrogate pairs
                buf[out += 1] = 0xed
                buf[out += 1] = ch
                buf[out += 1] = dat[pos += 1]
            else
                # Pick up surrogate pairs (CESU-8 format)
                ch = ((((((ch & 0x3f) << 6) | (dat[pos + 1] & 0x3f)) << 10)
                       + ((((dat[pos + 3] & 0x3f)%UInt32) << 6) | (dat[pos + 4] & 0x3f)))
                      - 0x01f0c00)
                pos += 4
                output_utf8_4byte!(buf, out, ch)
                out += 4
            end
        end
    end
    UTF8String(buf)
end

function convert(::Type{UTF8String}, a::Vector{UInt8}, invalids_as::AbstractString)
    l = length(a)
    idx = 1
    iscopy = false
    while idx <= l
        if !is_valid_continuation(a[idx])
            nextidx = idx+1+utf8_trailing[a[idx]+1]
            (nextidx <= (l+1)) && (idx = nextidx; continue)
        end
        !iscopy && (a = copy(a); iscopy = true)
        endn = idx
        while endn <= l
            !is_valid_continuation(a[endn]) && break
            endn += 1
        end
        (endn > idx) && (endn -= 1)
        splice!(a, idx:endn, Vector{UInt8}(invalids_as))
        l = length(a)
    end
    UTF8String(a)
end
convert(::Type{UTF8String}, s::AbstractString) = utf8(bytestring(s))

if isdefined(Base, :CodeUnits)
    convert(::Type{UTF8String}, s::Base.CodeUnits{UInt8,String}) = convert(UTF8String, Vector{UInt8}(s))
end

"""
Converts an already validated vector of `UInt16` or `UInt32` to a `UTF8String`

Input Arguments:

* `dat` Vector of code units (`UInt16` or `UInt32`), explicit `\0` is not converted
* `len` length of output in bytes

Returns:

* `UTF8String`
"""
function encode_to_utf8(::Type{T}, dat, len) where {T<:Union{UInt16, UInt32}}
    buf = Vector{UInt8}(undef, len)
    out = 0
    pos = 0
    @inbounds while out < len
        ch::UInt32 = dat[pos += 1]
        # Handle ASCII characters
        if ch <= 0x7f
            buf[out += 1] = ch
        # Handle 0x80-0x7ff
        elseif ch < 0x800
            buf[out += 1] = 0xc0 | (ch >>> 6)
            buf[out += 1] = 0x80 | (ch & 0x3f)
        # Handle 0x10000-0x10ffff (if input is UInt32)
        elseif ch > 0xffff # this is only for T == UInt32, should not be generated for UInt16
            output_utf8_4byte!(buf, out, ch)
            out += 4
        # Handle surrogate pairs
        elseif is_surrogate_codeunit(ch)
            output_utf8_4byte!(buf, out, get_supplementary(ch, dat[pos += 1]))
            out += 4
        # Handle 0x800-0xd7ff, 0xe000-0xffff UCS-2 characters
        else
            buf[out += 1] = 0xe0 | ((ch >>> 12) & 0x3f)
            buf[out += 1] = 0x80 | ((ch >>> 6) & 0x3f)
            buf[out += 1] = 0x80 | (ch & 0x3f)
        end
    end
    UTF8String(buf)
end

utf8(p::Ptr{UInt8}) =
    utf8(p, p == C_NULL ? Csize_t(0) : ccall(:strlen, Csize_t, (Ptr{UInt8},), p))
function utf8(p::Ptr{UInt8}, len::Integer)
    p == C_NULL && throw(ArgumentError("cannot convert NULL to string"))
    UTF8String(ccall(:jl_pchar_to_array, Vector{UInt8},
                     (Ptr{UInt8}, Csize_t), p, len))
end


text:
using Statistics
using LinearAlgebra
using Flux
using Flux: @functor
using LaTeXStrings
using KrylovMethods
using Distributions: Normal
using BSON: @load
using Plots
using Dates
using CSV, DataFrames
using Random

use_gpu = true
if use_gpu == true
    using CUDA
    #CUDA.allowscalar(false)
    cgpu = gpu
else
    cgpu = cpu
    pyplot()
end

pu = cpu # gpu
r_type = Float64
c_type = ComplexF64
u_type = Float32
gmres_type = ComplexF64
a_type = CuArray{gmres_type}
a_type = Array{gmres_type}
run_title = "64bit_cpu"

include("../src/multigrid/helmholtz_methods.jl")
include("../src/unet/model.jl")
include("../src/unet/data.jl")
include("../src/unet/train.jl")
include("../src/kappa_models.jl")
include("../src/gpu_krylov.jl")
include("test_utils.jl")

fgmres_func = KrylovMethods.fgmres # gpu_flexible_gmres #

# Test Parameters

level = 3
v2_iter = 10
gamma_val = 0.00001
pad_cells = [10;10]
point_sorce_results = true
check_unet_as_preconditioner = true
dataset_size = 1
blocks = 10

# Model Parameters

model_type = FFSDNUnet
k_type = FFKappa
resnet_type = SResidualBlock
arch = 2
k_chs = 10
indexes = 3
σ = elu
kernel = (3,3)
e_vcycle_input = false
kappa_type = 4 # 0 - uniform, 1 - CIFAR10, 2 - STL10
kappa_threshold = 25 # kappa ∈ [0.01*threshold, 1]
kappa_input = true
gamma_input = true
axb = false
norm_input = false
smooth = true
k_kernel = 5
before_jacobi = false
unet_in_vcycle = false

n = m = 128
h = r_type(2.0 / (n+m))
f = 10.0
restrt = 10
max_iter = 30
kappa = r_type.(generate_kappa!(n, m; type=kappa_type, smooth=smooth, threshold=kappa_threshold, kernel=k_kernel)|>pu) #ones(r_type,n-1,n-1)
omega = r_type(2*pi*f);
gamma = gamma_val*2*pi * ones(r_type,size(kappa))
gamma = r_type.(absorbing_layer!(gamma, pad_cells, omega))|>pu
bs = 10

retrain_size = 300
iter = 30

test_name = "09_40_06 SDNUnet1 g=-1 t=Float32 g=t e=f k=0 50 n=128 f=10_0 m=20000 bs=20 lr=0_0001 each=70 i=155"
test_name = "23_48_23 RADAM FFSDNUnet FFKappa SResidualBlock 10 elu 3 5 g=-1 t=Float32 g=t e=f r=f k=1 25 n=128 f=10_0 m=20000 bs=20 lr=0_0001 each=48 i=100"

sm_test_name = "23_48_23_$(run_title)_b$(blocks)_m$(kappa_type)_f$(Int32(f))_$(retrain_size)_$(iter)"
sm_test_name_r = "$(sm_test_name)_retrain"

model = load_model!("../models/$(test_name).bson", e_vcycle_input, kappa_input, gamma_input;kernel = kernel,model_type=model_type, k_type=k_type, resnet_type=resnet_type, k_chs=k_chs, indexes=indexes, σ=σ, arch=arch)

model_r128 = load_model!("../models/23_48_23 128 100 10 30 5 f f -1 r=f.bson", e_vcycle_input, kappa_input, gamma_input;kernel = kernel,model_type=model_type, k_type=k_type, resnet_type=resnet_type, k_chs=k_chs, indexes=indexes, σ=σ, arch=arch)
model_r256 = load_model!("../models/23_48_23 10 blocks 256 300 10 20 3 f f -1 r=f.bson", e_vcycle_input, kappa_input, gamma_input;kernel = kernel,model_type=model_type, k_type=k_type, resnet_type=resnet_type, k_chs=k_chs, indexes=indexes, σ=σ, arch=arch)
model_r512 = load_model!("../models/23_48_23 10 blocks 512 500 20 30 3 f f -1 r=f.bson", e_vcycle_input, kappa_input, gamma_input;kernel = kernel,model_type=model_type, k_type=k_type, resnet_type=resnet_type, k_chs=k_chs, indexes=indexes, σ=σ, arch=arch)

# model_r128, _ = model_tuning!(model1, sm_test_name, kappa, omega, gamma, n, m, f, retrain_size, bs, iter, 0.001, kappa_threshold, false, false, k_kernel, -1,kappa_type;residual_loss=false)

if point_sorce_results == false
    check_point_source_problem!("$(Dates.format(now(), "HH_MM_SS")) $(sm_test_name)", model, n, m, kappa, omega, gamma, e_vcycle_input, kappa_input, gamma_input; v2_iter=v2_iter, level=level)
end
if check_unet_as_preconditioner == true #gmres_alternatively_
    # check_model_times!(sm_test_name, model, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
    check_model!(sm_test_name, model, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
    check_model!(sm_test_name_r, model_r128, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
end

n = m = 256
f = 20.0
restrt = 15
max_iter = 40
k_kernel = 3
kappa = r_type.((generate_kappa!(n, m; type=kappa_type, smooth=true, threshold=kappa_threshold, kernel=k_kernel))|>pu) #ones(r_type,n-1,n-1)
omega = r_type(2*pi*f)
gamma = gamma_val*2*pi * ones(r_type,size(kappa))
gamma = r_type.((absorbing_layer!(gamma, pad_cells, omega))|>pu)
bs = 10
iter = 40
# model, _ = model_tuning!(model, sm_test_name, kappa, omega, gamma, n, f, 50, bs, iter, 0.001, kappa_threshold, false, false, k_kernel, -1)
# model, sm_test_name = model_tuning!(model, sm_test_name, kappa, omega, gamma, n, f,500, bs, iter, 0.001, kappa_threshold, false, false, k_kernel, -1;residual_loss=false)

if point_sorce_results == false
    check_point_source_problem!("$(Dates.format(now(), "HH_MM_SS")) $(sm_test_name)", model, n, m, kappa, omega, gamma, e_vcycle_input, kappa_input, gamma_input; v2_iter=v2_iter, level=level)
end
if check_unet_as_preconditioner == true
    # check_model_times!("$(sm_test_name)", model, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
    check_model!(sm_test_name, model, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
    check_model!(sm_test_name_r, model_r256, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
end

n = m = 512
f = 40.0
restrt = 20
max_iter = 50
k_kernel = 3
kappa = r_type.(generate_kappa!(n, m; type=kappa_type, smooth=false, threshold=kappa_threshold, kernel=k_kernel))|>pu
omega = r_type(2*pi*f)
gamma = gamma_val*2*pi * ones(r_type,size(kappa))
gamma = r_type.(absorbing_layer!(gamma, pad_cells, omega))|>pu
bs = 20
iter = 40
# model, _ = model_tuning!(model, sm_test_name, kappa, omega, gamma, n, f, 50, bs, iter, 0.001, kappa_threshold, false, false, k_kernel, -1)
# model, sm_test_name = model_tuning!(model, sm_test_name, kappa, omega, gamma, n, f, 500, bs, iter, 0.001, kappa_threshold, false, false, k_kernel, -1;residual_loss=false)

if point_sorce_results == false
    check_point_source_problem!(test_name, model, n, m, kappa, omega, gamma, e_vcycle_input, kappa_input, gamma_input; v2_iter=v2_iter, level=level)
end
if check_unet_as_preconditioner == true
    # check_model_times!("$(sm_test_name)", model, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
    check_model!(sm_test_name, model, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
    check_model!(sm_test_name_r, model_r512, n, m, f, kappa, omega, gamma, e_vcycle_input, kappa_type, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=smooth, k_kernel=k_kernel, threshold=kappa_threshold, axb=axb, norm_input=norm_input, before_jacobi=false,log_error=false,unet_in_vcycle=unet_in_vcycle,indexes=indexes, arch=arch)
end


text:
"""
    SteadyStateCallback(; abstol=1.0e-8, reltol=1.0e-6)

Terminates the integration when the [`residual_steady_state(du, equations)`](@ref)
falls below the threshold specified by `abstol, reltol`.
"""
mutable struct SteadyStateCallback{RealT<:Real}
  abstol::RealT
  reltol::RealT
end

function SteadyStateCallback(; abstol=1.0e-8, reltol=1.0e-6)
  abstol, reltol = promote(abstol, reltol)
  steady_state_callback = SteadyStateCallback(abstol, reltol)

  DiscreteCallback(steady_state_callback, steady_state_callback,
                   save_positions=(false,false))
end

function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SteadyStateCallback})
  @nospecialize cb # reduce precompilation time

  steady_state_callback = cb.affect!
  print(io, "SteadyStateCallback(abstol=", steady_state_callback.abstol, ", ",
                                "reltol=", steady_state_callback.reltol, ")")
end

function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:SteadyStateCallback})
  @nospecialize cb # reduce precompilation time

  if get(io, :compact, false)
    show(io, cb)
  else
    steady_state_callback = cb.affect!

    setup = [
             "absolute tolerance" => steady_state_callback.abstol,
             "relative tolerance" => steady_state_callback.reltol,
            ]
    summary_box(io, "SteadyStateCallback", setup)
  end
end

# affect!
(::SteadyStateCallback)(integrator) = terminate!(integrator)

# the condition
function (steady_state_callback::SteadyStateCallback)(u_ode::AbstractVector, t, integrator)
  semi = integrator.p

  u  = wrap_array(u_ode, semi)
  du = wrap_array(get_du(integrator), semi)
  terminate = steady_state_callback(du, u, semi)
  if mpi_isparallel()
    # MPI.jl doesn't seem to have MPI_C_BOOL
    terminate_integer = Int(terminate)
    terminate = !iszero(MPI.Allreduce!(Ref(terminate_integer), +, mpi_comm())[])
  end
  if mpi_isroot() && terminate
    @info "  Steady state tolerance reached" steady_state_callback t
  end
  return terminate
end

function (steady_state_callback::SteadyStateCallback)(du, u, semi::AbstractSemidiscretization)
  steady_state_callback(du, u, mesh_equations_solver_cache(semi)...)
end

include("steady_state_dg1d.jl")
include("steady_state_dg2d.jl")
include("steady_state_dg3d.jl")


text:
function testNHBaroWaveCart()
close all
clear all

% Grid
nz=20;
Param.Lx=30000*1.e3;
Param.Ly=6000*1.e3;
Param.H=30*1.e3;
dx=400*1.e3;
dy=400*1.e3;
nx=Param.Lx/dx;
ny=2*Param.Ly/dy;
x0=0;
y0=-Param.Ly;

Param.Boundary.WE='Period';
Param.Boundary.BT='Period';
Param.hS='';
Param.Grid=CartGrid(nx,ny,Param.Lx,2*Param.Ly,x0,y0,OrientFaceCart,Param.Boundary,Param);

Param.Grid.nz=nz;
Param.Grid.zP=zeros(nz,1);
Param.Grid.z=zeros(nz+1,1);
Param.Grid.dz=Param.H/nz;
Param.Grid.zP(1)=Param.Grid.dz/2;
for i=2:nz
  Param.Grid.zP(i)=Param.Grid.zP(i-1)+Param.Grid.dz;
end
for i=2:nz+1
  Param.Grid.z(i)=Param.Grid.z(i-1)+Param.Grid.dz;
end

% Model
LRef=11*1.e5;
Param.HyperD=0; %1.e14*(dx/LRef)^3.2;
Param.StrideDamp=6000;
Param.Relax=5.e-1;
Param.Damping=true;
Param.Flat=false;
Param.Coriolis=true;
Param.CoriolisType='Beta-Plane';
Param.Buoyancy=true;
Param.Th0=300;
Param.uMax=10;
Param.vMax=0;
Param.NBr=1.e-2;
Param.lat0=pi/2;
Param.Equation='Compressible';
Param.TopoS='';
Param.Omega=2*pi/(24*3600);
Param.RadEarth=6.37122d+6;
Param.f0=2*sin(Param.lat0)*Param.Omega;
Param.beta0=2*cos(Param.lat0)*Param.Omega/Param.RadEarth;
Param.y0=Param.Ly/2;
Param.xC=2000*1.e3;
Param.yC=2500*1.e3;
Param.Lp=600*1.e3;
Param.uP=1;
Param.b=2;
Param.Lapse=0.005;
Param.T0=288;
Param.u0=35;
Param.ProfVel='BaroWaveCart';
Param.ProfRho='BaroWaveCart';
Param.ProfTheta='BaroWaveCart';
Param.NumV=5;
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
Param.Thermo='Energy';

% Physical parameters
Param.Grav=9.81;
Param.cS=360;
Param.Grav=9.81d0;
Param.Cpd=1004.0d0;
Param.Cvd=717.0d0;
Param.Rd=Param.Cpd-Param.Cvd;
Param.p0=1.0d5;
Param.Cpv=1885.0d0;
Param.Gamma=Param.Cpd/Param.Cvd;
Param.kappa=Param.Rd/Param.Cpd;

% Output
Param.Flat=false;
Param.level=1;
Param.fig=1;
Param.SliceYZ.Type='YZ';
Param.SliceYZ.x=x0+Param.Lx/nx/2;
Param.SliceYZ.ix=10;
Param.SliceXY.Type='XY';
Param.SliceXY.z=1500;
Param.SliceXY.iz=1;

% Discretization

OrdPoly=3;
[CG,Param]=Discretization(OrdPoly,@JacobiDG3,Param);

% Initial conditions

U=zeros(CG.NumG,nz,Param.NumV);
U(:,:,Param.RhoPos)=Project(@fRho,CG,Param);
[U(:,:,Param.uPos),U(:,:,Param.vPos)]=ProjectVec(@fVel,CG,Param);
U(:,:,Param.ThPos)=Project(@fTheta,CG,Param).*U(:,:,Param.RhoPos);
if strcmp(Param.Thermo,'Energy')
  U(:,:,Param.ThPos)=PotToEnergy(U,CG,Param);
end
ThB=Project(@fThetaBGrd,CG,Param);
Param.fig=PlotCG(U(:,:,Param.uPos)...
  ,CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXY);

% Integration
CFL=0.125;
dtau=240;
time=0;

IntMethod='Rosenbrock'; %'RungeKutta';
%IntMethod='RungeKutta';
nIter=20000;
%v = VideoWriter ('Galewsky.avi');
%open (v);
Param.RK=RungeKuttaMethod('RK4');
Param.ROS=RosenbrockMethod('ROSRK3');
switch IntMethod
  case 'Rosenbrock'
    for i=1:nIter
      U=Rosenbrock(U,dtau,@FcnNHCurlVec,@Jac,CG,Param);
      time=time+dtau;
      if mod(i,360)==0
        if strcmp(Param.Thermo,'Energy')
          Th=EnergyToPot(U,CG,Param);
        else
          Th=U(:,:,Param.ThPos);
        end
        Param.fig=PlotCG(Th./U(:,:,Param.RhoPos)-ThB...
          ,CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXY);
%         Param.fig=PlotCG(U(:,:,Param.uPos),CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXY);
%         Param.fig=PlotCG(U(:,:,Param.vPos),CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXY);
        W=zeros(size(U,1),nz+1);
        W(:,2:nz+1)=U(:,:,Param.wPos);
        Param.fig=PlotCG(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,@TransCart...
          ,@Topo,Param,Param.fig,Param.SliceXY);
      end
    end
  case 'RungeKutta'
    for i=1:nIter
      U=RungeKuttaExplicit(U,dtau,@FcnNHCurlVec,CG,Param);

      time=time+dtau;
      if mod(i,1)==0
        Param.fig=PlotCG(U(:,:,Param.ThPos)./U(:,:,Param.RhoPos)...
          ,CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXY);
        Param.fig=PlotCG(U(:,:,Param.uPos),CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);
        Param.fig=PlotCG(U(:,:,Param.vPos),CG,@TransCart,@Topo,Param,Param.fig,Param.SliceXZ);
        W=zeros(size(U,1),nz+1);
        W(:,2:nz+1)=U(:,:,Param.wPos);
        Param.fig=PlotCG(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,@TransCart...
          ,@Topo,Param,Param.fig,Param.SliceXZ);
        frame=getframe(gcf);
        writeVideo (v, frame);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
%close(v)
end


text:
@enum OpCode::UInt8 begin
    # ============================
    # PROTOCOL 0/1

    # push special markobeject on stack
    OP_MARK = 0x28

    # every pickle ends with STOP
    OP_STOP = 0x2E

    # discard topmost stack item
    OP_POP = 0x30

    # discard stack top through topmost markobject
    OP_POP_MARK = 0x31

    # duplicate top stack item
    OP_DUP = 0x32

    # push float object; decimal string argument
    OP_FLOAT = 0x36

    # push integer or bool; decimal string argument
    OP_INT = 0x49

    # push four-byte signed int
    OP_BININT = 0x4A

    # push 1-byte unsigned int
    OP_BININT1 = 0x4B

    # push long; decimal string argument
    OP_LONG = 0x4C

    # push 2-byte unsigned int
    OP_BININT2 = 0x4D

    # push None
    OP_NONE = 0x4E

    # push persistent object; id is taken from string arg
    OP_PERSID = 0x50

    # push persistent object; id is taken from stack
    OP_BINPERSID = 0x51

    # apply callable to argtuple, both on stack
    OP_REDUCE = 0x52

    # push string; NL-terminated string argument
    OP_STRING = 0x53

    # push string; counted binary string argument
    OP_BINSTRING = 0x54

    # push string; counted binary string argument 256 bytes
    OP_SHORT_BINSTRING = 0x55

    # push Unicode string; raw-unicode-escaped argument
    OP_UNICODE = 0x56

    # push Unicode string; counted UTF-8 string argument
    OP_BINUNICODE = 0x58

    # append stack top to list below it
    OP_APPEND = 0x61

    # call __setstate__ or __dict__.update()
    OP_BUILD = 0x62

    # push self.find_class(modname, name); 2 string args
    OP_GLOBAL = 0x63

    # build a dict from stack items
    OP_DICT = 0x64

    # push empty dict
    OP_EMPTY_DICT = 0x7d

    # extend list on stack by topmost stack slice
    OP_APPENDS = 0x65

    # push item from memo on stack; index is string arg
    OP_GET = 0x67

    # push item from memo on stack; index is 1-byte arg
    OP_BINGET = 0x68

    # build & push class instance
    OP_INST = 0x69

    # push item from memo on stack; index is 4-byte arg
    OP_LONG_BINGET = 0x6a

    # build list from topmost stack items
    OP_LIST = 0x6c

    # push empty list
    OP_EMPTY_LIST = 0x5d

    # build & push class instance
    OP_OBJ = 0x6f

    # store stack top in memo; index is string arg
    OP_PUT = 0x70

    # store stack top in memo; index is 1-byte arg
    OP_BINPUT = 0x71

    # store stack top in memo; index is 4-byte arg
    OP_LONG_BINPUT = 0x72

    # add key+value pair to dict
    OP_SETITEM = 0x73

    # build tuple from topmost stack items
    OP_TUPLE = 0x74

    # push empty tuple
    OP_EMPTY_TUPLE = 0x29

    # modify dict by adding topmost key+value pairs
    OP_SETITEMS = 0x75

    # push float; arg is 8-byte float encoding
    OP_BINFLOAT = 0x47

    # ============================
    # PROTOCOL 2

    # identify pickle protocol
    OP_PROTO = 0x80

    # build object by applying cls.__new__ to argtuple
    OP_NEWOBJ = 0x81

    # push object from extension registry; 1-byte index
    OP_EXT1 = 0x82

    # ditto, but 2-byte index
    OP_EXT2 = 0x83

    # ditto, but 4-byte index
    OP_EXT4 = 0x84

    # build 1-tuple from stack top
    OP_TUPLE1 = 0x85

    # build 2-tuple from two topmost stack items
    OP_TUPLE2 = 0x86

    # build 3-tuple from three topmost stack items
    OP_TUPLE3 = 0x87

    # push True
    OP_NEWTRUE = 0x88

    # push False
    OP_NEWFALSE = 0x89

    # push long from < 256 bytes
    OP_LONG1 = 0x8a

    # push really big long
    OP_LONG4 = 0x8b

    # ============================
    # PROTOCOL 3

    # push bytes; counted binary string argument
    OP_BINBYTES = 0x42

    # push bytes; counted binary string argument < 256 bytes
    OP_SHORT_BINBYTES = 0x43

    # ============================
    # PROTOCOL 4

    # push short string; UTF-8 length < 256 bytes
    OP_SHORT_BINUNICODE = 0x8c

    # push very long string
    OP_BINUNICODE8 = 0x8d

    # push very long bytes string
    OP_BINBYTES8 = 0x8e

    # push empty set on the stack
    OP_EMPTY_SET = 0x8f

    # modify set by adding topmost stack items
    OP_ADDITEMS = 0x90

    # build frozenset from topmost stack items
    OP_FROZENSET = 0x91

    # like NEWOBJ but work with keyword only arguments
    OP_NEWOBJ_EX = 0x92

    # same as GLOBAL but using names on the stacks
    OP_STACK_GLOBAL = 0x93

    # store top of the stack in memo
    OP_MEMOIZE = 0x94

    # indicate the beginning of a new frame
    OP_FRAME = 0x95

    # ============================
    # PROTOCOL 5

    # push bytearray
    OP_BYTEARRAY8 = 0x96

    # push next out-of-band buffer
    OP_NEXT_BUFFER = 0x97

    # make top of stack readonly
    OP_READONLY_BUFFER = 0x98

end


text:
module Covid
include("Stats.jl") #incluimos los archivos con las funciones
include("Utils.jl")
end

max_stars_repo_name: MarleenMs/Covid

max_issues_repo_name: MarleenMs/Covid

max_forks_repo_name: MarleenMs/Covid


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct RemindersApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_remindersAdd(_api::RemindersApi, text::String, time::String; user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/reminders.add", ["slackAuth"])
    Swagger.set_param(_ctx.form, "text", text)  # type String
    Swagger.set_param(_ctx.form, "time", time)  # type String
    Swagger.set_param(_ctx.form, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Creates a reminder.
Param: text::String (required)
Param: time::String (required)
Param: user::String
Return: Dict{String, Any}
"""
function remindersAdd(_api::RemindersApi, text::String, time::String; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersAdd(_api, text, time; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function remindersAdd(_api::RemindersApi, response_stream::Channel, text::String, time::String; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersAdd(_api, text, time; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_remindersComplete(_api::RemindersApi; reminder=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/reminders.complete", ["slackAuth"])
    Swagger.set_param(_ctx.form, "reminder", reminder)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Marks a reminder as complete.
Param: reminder::String
Return: Dict{String, Any}
"""
function remindersComplete(_api::RemindersApi; reminder=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersComplete(_api; reminder=reminder, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function remindersComplete(_api::RemindersApi, response_stream::Channel; reminder=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersComplete(_api; reminder=reminder, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_remindersDelete(_api::RemindersApi; reminder=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/reminders.delete", ["slackAuth"])
    Swagger.set_param(_ctx.form, "reminder", reminder)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Deletes a reminder.
Param: reminder::String
Return: Dict{String, Any}
"""
function remindersDelete(_api::RemindersApi; reminder=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersDelete(_api; reminder=reminder, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function remindersDelete(_api::RemindersApi, response_stream::Channel; reminder=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersDelete(_api; reminder=reminder, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_remindersInfo(_api::RemindersApi; reminder=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/reminders.info", ["slackAuth"])
    Swagger.set_param(_ctx.query, "reminder", reminder)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets information about a reminder.
Param: reminder::String
Return: Dict{String, Any}
"""
function remindersInfo(_api::RemindersApi; reminder=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersInfo(_api; reminder=reminder, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function remindersInfo(_api::RemindersApi, response_stream::Channel; reminder=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_remindersInfo(_api; reminder=reminder, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_remindersList(_api::RemindersApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/reminders.list", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Lists all reminders created by or for a given user.
Return: Dict{String, Any}
"""
function remindersList(_api::RemindersApi; _mediaType=nothing)
    _ctx = _swaggerinternal_remindersList(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function remindersList(_api::RemindersApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_remindersList(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export remindersAdd, remindersComplete, remindersDelete, remindersInfo, remindersList


text:
using Makie
using MakieLayout

begin
    scene = Scene(camera=campixel!)
    screen = display(scene)

    outergrid = GridLayout(scene, alignmode=Outside(30))

    innergrid = outergrid[1, 1] = GridLayout(2, 2)

    las = innergrid[:, :] = [
        LAxis(scene, xautolimitmargin=(0, 0), yautolimitmargin=(0, 0))
        for i in 1:2, j in 1:2]

    xs = LinRange(0, 3, 100)
    ys = LinRange(0, 3, 100)
    vals = [sin(x * y) for x in xs, y in ys]

    gl1 = gridnest!(innergrid, 1, 1)

    hm1 = heatmap!(las[1, 1], vals, colormap=:viridis)
    cb1 = gl1[1, 2] = LColorbar(scene, hm1, width=30f0,
        height=Relative(0.66), alignment=(:center, :center), label = "amplitude",
        ticklabelspace = 60f0)

    gl2 = gridnest!(innergrid, 1, 2)

    hm2 = heatmap!(las[1, 2], vals, colormap=:heat)

    gl2[2, 1] = LSlider(scene, height = 30, range = 0:100)

    gl3 = gridnest!(innergrid, 2, 1)

    las[2, 1].titlevisible = false
    hm3 = heatmap!(las[2, 1], vals, colormap=:inferno)
    cb3 = gl3[0, 1] = LColorbar(scene, hm3, vertical = false,
        flipaxisposition=true, ticklabelalign = (:center, :bottom), height=30f0,
        label = "amplitude")

    gl4 = gridnest!(innergrid, 2, 2)

    las[2, 2].yaxisposition = :right
    las[2, 2].yticklabelalign = (:left, :center)

    hm4 = heatmap!(las[2, 2], vals, colormap=:rainbow)
    cb4 = gl4[1, 0] = LColorbar(scene, hm4,
        flipaxisposition = false, ticklabelalign = (:right, :center), width=30f0,
        label = "amplitude", ticklabelspace=60)

    innergrid[0, :] = LText(scene, text="Colorbars", textsize=50)

    cb2 = innergrid[0, :] = LColorbar(scene, hm2, vertical = false,
        flipaxisposition=false, ticklabelalign = (:center, :top), height=30f0,
        alignment=(:center, :center), label = "amplitude")
end


text:
struct Bioequivalence
    data
    design
    dropout
    subject
    period
    formulation
    sequence
    reference
    subjects
    periods
    formulations
    sequences
end

function bioquivalence(data;
    variable = nothing,
    subject = :subject,
    period = :period, formulation = :formulation, sequence = :sequence, reference = nothing, design = nothing, io::IO = stdout, seqcheck = true, dropcheck = true)

    dfnames = Symbol.(names(data))

    fac = [subject, formulation]

    fac ⊆ dfnames || error("Subject or formulation column not found in dataframe!")

    disallowmissing!(data, subject)

    disallowmissing!(data, formulation)

    subjects = unique(data[!, subject])

    formulations = sort!(unique(data[!, formulation]))

    if isnothing(reference)
        @info "Reference formulation not specified. First used: \"$(first(formulations))\"."
        reference = first(formulations)
    else
        reference ∈ formulations || error("Reference formulation \"$(reference)\" not found in dataframe.")
    end

    local dropout = nothing
    local periods = nothing
    local sequences = nothing

    if isnothing(period) && isnothing(sequence) && isnothing(design)
        length(subjects) == length(data[!, subject]) || error("Trial design seems parallel, but subjects not unique!")
        design = :parallel
        println(io, "Parallel desigh used.")
    end

    if isnothing(design) || design != :parallel

        !isnothing(period) || error("Trial design seems NOT parallel, but period is nothing")

        !isnothing(sequence) || error("Trial design seems NOT parallel, but sequence is nothing")

        period ∈ dfnames || error("Period not found in dataframe!")

        sequence ∈ dfnames || error("Sequence not found in dataframe!")

        periods = unique(data[!, period])

        sequences = unique(data[!, sequence])

        push!(fac, period, sequence)

        disallowmissing!(data, period)

        disallowmissing!(data, sequence)

        unstdata = unstack(data, [subject, sequence], period, formulation)

        pnames = Symbol.(periods)
        if dropcheck
            if !all(completecases(unstdata, pnames))
                dropout = true
                dropmissing!(unstdata)
                 @info "Dropuot(s) found in dataframe!"
            else
                @info "No dropuot(s) found in dataframe!"
                dropout = false
            end
        end

        if seqcheck
            seqdict = Dict()
            for i  in 1:size(unstdata, 1)
                if ht_keyindex(seqdict, unstdata[i, pnames]) > 0
                    if seqdict[unstdata[i, pnames]] != unstdata[i, sequence] error("Sequence error!") end
                else
                    seqdict[unstdata[i, pnames]] = unstdata[i, sequence]
                end
            end
            @info "Sequences looks correct..."
        end

        if isnothing(design)
            @info "Trying to find out the design..."
            design = Symbol("$(length(formulations))X$(length(sequences))X$(length(periods))")
            @info  "Seems design type is: $design"
        else
            spldes = split(uppercase(string(design)), "X")
            if length(spldes) != 3 &&  uppercase(string(design)) != "2X2" error("Unknown design type. Use fXsXp format or \"2Х2\".") end
            if length(formulations) != parse(Int, spldes[1]) error("Design error: formulations count wrong!") end
            if length(sequences) != parse(Int, spldes[2]) error("Design error: sequences count wrong!") end
            if length(periods) != parse(Int, spldes[3]) error("Design error: periods count wrong!") end
            @info "Design type seems fine..."
        end
    end

    Bioequivalence(
        data,
        design,
        dropout,
        subject,
        period,
        formulation,
        sequence,
        reference,
        subjects,
        periods,
        formulations,
        sequences)
end


text:
using QuantLib

tests = ["cash_flows",
        "indexes",
        "bonds",
        "forwards",
        "options",
        "swaps",
        "swaptions",
        "fd_mesher",
        "curves"]

println("Running tests:")

for t in tests
    tfile = string(t, ".jl")
    println(" * $(tfile) ...")
    include(tfile)
end


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# camera/format.jl: Camera object image format features

"""
  adcbits(::Camera) -> String

  Return ADC bit depth.
"""
adcbits(cam::Camera) = get(SpinEnumNode(cam, "AdcBitDepth"))

"""
  adcbits!(::Camera, ::AbstractString) -> String

  Set ADC bit depth, returns set depth.
"""
adcbits!(cam::Camera, bits) = set!(SpinEnumNode(cam, "AdcBitDepth"), bits)

"""
  pixelformat(::Camera) -> String

  Return camera pixel format.
"""
pixelformat(cam::Camera) = get(SpinEnumNode(cam, "PixelFormat"))

"""
  pixelformat!(::Camera, ::AbstractString) -> String

  Set camera pixel format, returns set format.
"""
pixelformat!(cam::Camera, fmt) = set!(SpinEnumNode(cam, "PixelFormat"), fmt)

"""
  sensordims(::Camera) -> (Int, Int)

  Return the width and height of the sensor.
"""
sensordims(cam::Camera) = (get(SpinIntegerNode(cam, "SensorWidth")), get(SpinIntegerNode(cam, "SensorHeight")))

"""
  imagedims(::Camera) -> (Int, Int)

  Return the width and height of the image.
"""
imagedims(cam::Camera) = (get(SpinIntegerNode(cam, "Width")), get(SpinIntegerNode(cam, "Height")))

"""
  imagedims!(::Camera, (width, height)) -> (Int, Int)

  Set the width and height of the image, return the set values.
"""
function imagedims!(cam::Camera, dims)
  set!(SpinIntegerNode(cam, "Width"), dims[1])
  set!(SpinIntegerNode(cam, "Height"), dims[2])
  imagedims(cam)
end

"""
  imagedims_limits(::Camera) -> ((Int, Int), (Int, Int))

  Image dimension limits in width, and height.
"""
imagedims_limits(cam::Camera) = (range(SpinIntegerNode(cam, "Width")), range(SpinIntegerNode(cam, "Height")))

"""
  offsetdims(::Camera) -> (Int, Int)

  Return the offset in x and y of the image.
"""
offsetdims(cam::Camera) = (get(SpinIntegerNode(cam, "OffsetX")), get(SpinIntegerNode(cam, "OffsetY")))

"""
  offsetdims!(::Camera, (X,Y)) -> (Int, Int)

  Set the image offset in x and y, return the set values.
"""
function offsetdims!(cam::Camera, dims)
  set!(SpinIntegerNode(cam, "OffsetX"), dims[1])
  set!(SpinIntegerNode(cam, "OffsetY"), dims[2])
  offsetdims(cam)
end

"""
  offsetdims_limits(::Camera) -> ((Int, Int), (Int, Int))

  Image offsets limits in width, and height.
"""
offsetdims_limits(cam::Camera) = (range(SpinIntegerNode(cam, "OffsetX")), range(SpinIntegerNode(cam, "OffsetY")))


text:
using JuDGE

leafnodesA = [ [1,1,1,1],
               [1,1,2,1],
               [1,2,1,1,1],
               [1,2,1,1,2],
               [1,2,1,2],
               [1,2,1,3] ]

probsA = [0.2,0.3,0.1,0.2,0.05,0.15]

treeA, probabilityA = tree_from_leaves(leafnodesA,probsA)

nodesA = collect(treeA)

probabilityA[nodesA[10]]
probabilityA[get_node(treeA,leafnodesA[3])]
probabilityA[get_node(treeA,[1,2,1,1,1])]
JuDGE.print_tree(treeA)

nodeprobsB = [1.0,
               [0.5,
                 [0.4,
                   [1.0]],
                 [0.6,
                   [1.0]]],
               [0.5,
                 [1.0,
                   [0.6,
                     [1/3],
                     [2/3]],
                   [0.1],
                   [0.3]]]]

treeB,probabilityB = tree_from_nodes(nodeprobsB)

nodesB = collect(treeB)

probabilityB[treeB.children[2].children[1].children[1].children[1]]
probabilityB[nodesB[10]]
JuDGE.print_tree(treeB)

treeC,data = tree_from_file(joinpath(@__DIR__,"treeC.tree"))

nodesC = collect(treeC)
JuDGE.print_tree(treeC,data[:pr])
probabilityC=JuDGE.convert_probabilities(treeC,data[:pr])
JuDGE.print_tree(treeC,probabilityC)
probabilityC[get_node(treeC,[1,2,1,1,1])]

data[:pr_absolute]=probabilityC
JuDGE.visualize_tree(treeC,data)


text:
using AeroMDAO

## Fuselage definition
rads = range(0, 1, length = 10)
lens = 0.05:0.05:1.0
fuse = Fuselage(lens[1:end-1], [ rads; reverse(rads) ] )

## Property checks
projected_area(fuse)
length(fuse)

## Cosine interpolation
lens_rads = cosine_spacing(fuse, 40)

# Circles for plotting
n_pts          = 20
circle3D(r, n) = let arcs = 0:2π/n:2π; [ zeros(length(arcs)) r * cos.(arcs) r * sin.(arcs) ] end

##
xs = lens_rads[:,1]
circs = [ reduce(hcat, eachrow(circ) .+ Ref([x; 0; 0]))' for (x, circ) in zip(xs, circle3D.(lens_rads[:,2], n_pts)) ]

##
using Plots
pyplot()
plt = plot()
[ plot!(circ[:,1], circ[:,2], circ[:,3], color = :blue, label = :none) for circ in circs ]
plot!()


text:
using Base.Test

using Celeste: Model, SensitiveFloats
using Celeste.DeterministicVI.ConstraintTransforms: ParameterConstraint,
                        ConstraintBatch, BoxConstraint, SimplexConstraint
using Celeste.DeterministicVI: ElboArgs
using Celeste.DeterministicVI.ElboMaximize: ElboConfig, maximize!, elbo_optim_options
using Optim

function verify_sample_galaxy(vs, pos)
    @test vs[ids.is_star[2]] >= 0.99

    @test isapprox(vs[ids.pos[1]], pos[1], atol=0.1)
    @test isapprox(vs[ids.pos[2]], pos[2], atol=0.1)

    @test isapprox(vs[ids.gal_axis_ratio] , 0.7, atol=0.05)
    @test isapprox(vs[ids.gal_frac_dev]  , 0.1, atol=0.08)
    @test isapprox(vs[ids.gal_radius_px], 4.0, atol=0.2)

    phi_hat = vs[ids.gal_angle]
    phi_hat -= floor(phi_hat / pi) * pi
    five_deg = 5 * pi/180
    @test isapprox(phi_hat, pi/4, atol=five_deg)

    brightness_hat = exp(vs[ids.flux_loc[2]] + 0.5 * vs[ids.flux_scale[2]])
    @test isapprox(brightness_hat / sample_galaxy_fluxes[3], 1.0, atol=0.05)

    true_colors = log.(sample_galaxy_fluxes[2:5] ./ sample_galaxy_fluxes[1:4])
    for b in 1:4
        @test isapprox(vs[ids.color_mean[b, 2]], true_colors[b], atol=0.2)
    end
end

#########################################################

function test_single_source_optimization()
    ea, vp, catalog = gen_three_body_dataset();

    s = 2
    ea = make_elbo_args(ea.images, catalog, active_source=s, include_kl=false);
    vp_original = deepcopy(vp);

    cfg = ElboConfig(ea, vp; loc_width=1.0)
    maximize!(ea, vp, cfg)

    # Test that it only optimized source s
    @test vp[s] != vp_original[s]
    for other_s in setdiff(1:ea.S, s)
        @test vp[other_s] ≈ vp_original[other_s]
    end
end

function test_galaxy_optimization()
    ea, vp, catalog = gen_sample_galaxy_dataset(; include_kl = false);
    cfg = ElboConfig(ea, vp; loc_width=3.0)
    maximize!(ea, vp, cfg)
    verify_sample_galaxy(vp[1], [8.5, 9.6])
end

function test_full_elbo_optimization()
    ea, vp, catalog = gen_sample_galaxy_dataset(perturb=true);
    cfg = ElboConfig(ea, vp; loc_width=1.0,
                 optim_options=elbo_optim_options(xtol_abs=0.0))
    maximize!(ea, vp, cfg)
    verify_sample_galaxy(vp[1], [8.5, 9.6]);
end

@testset "optimization" begin
    test_single_source_optimization()
    test_full_elbo_optimization()
    test_galaxy_optimization()
end


text:
using Base.Test
using CUJulia

path = joinpath(dirname(@__FILE__), "cudnn.jl")
include(path)


text:
in_channel_v = 6
out_channel_v = 30
in_channel_e = 1
out_channel_e = 4
N = 4
T = Float32
adj = T[0. 1. 0. 1.;
        1. 0. 1. 0.;
        0. 1. 0. 1.;
        1. 0. 1. 0.]

adj_single_vertex =   T[0. 0. 0. 1.;
                        0. 0. 0. 0.;
                        0. 0. 0. 1.;
                        1. 0. 1. 0.]

@testset "EdgeFtLayer" begin
    @testset "Constructor" begin
        ft_layer = SeaPearl.EdgeFtLayer(;v_dim = in_channel_v=>out_channel_v, e_dim = in_channel_e=>out_channel_e)

        @test size(ft_layer.W_a) == (out_channel_v, 2 * in_channel_v + in_channel_e)
        @test size(ft_layer.W_T) == (out_channel_v, 2 * in_channel_v + in_channel_e)
        @test size(ft_layer.b_T) == (out_channel_v,)
        @test size(ft_layer.W_e) == (out_channel_e, in_channel_v)
        @test size(ft_layer.W_ee) == (out_channel_e, in_channel_e)
        @test isa(ft_layer.prelu_α, T)
    end

    @testset "EdgeFtLayer(fg)" begin
        ft_layer = SeaPearl.EdgeFtLayer(;v_dim = in_channel_v=>out_channel_v, e_dim = in_channel_e=>out_channel_e)

        nf = rand(Float32, in_channel_v, N)
        ef = rand(Float32, in_channel_e, 8)
        fg = GraphSignals.FeaturedGraph(adj; nf=nf, ef=ef)

        fg_ = ft_layer(fg)

        @test size(GraphSignals.node_feature(fg_)) == (out_channel_v, N)
        @test size(GraphSignals.edge_feature(fg_)) == (out_channel_e, 8)

        g = Zygote.gradient(x -> sum(node_feature(ft_layer(x))), fg)[1]
        @test size(g[].nf) == size(nf)

        ft_layer2 = SeaPearl.EdgeFtLayer(;v_dim = out_channel_v=>10, e_dim = out_channel_e=>4)
        fg2_ = ft_layer2(fg_)

        @test size(GraphSignals.node_feature(fg2_)) == (10, N)
        @test size(GraphSignals.edge_feature(fg2_)) == (4, 8)
    end

    @testset "message()" begin
        ft_layer = SeaPearl.EdgeFtLayer(;v_dim = in_channel_v=>out_channel_v, e_dim = in_channel_e=>out_channel_e)

        nf = rand(Float32, in_channel_v, N)
        ef = rand(Float32, in_channel_e, 8)
        message_vec = GeometricFlux.message(ft_layer, nf[:, 1], nf[:, 2], ef[:, 1])

        @test size(message_vec) == (out_channel_v*2 + out_channel_e,)
    end

    @testset "apply_batch_message()" begin
        ft_layer = SeaPearl.EdgeFtLayer(;v_dim = in_channel_v=>out_channel_v, e_dim = in_channel_e=>out_channel_e)
        nf = rand(Float32, in_channel_v, N)
        ef = rand(Float32, in_channel_e, 8)
        fg = GraphSignals.FeaturedGraph(adj; nf=nf, ef=ef)

        adj_list = GeometricFlux.adjacency_list(fg)

        batch_message = GeometricFlux.apply_batch_message(ft_layer, 2, adj_list[2], GeometricFlux.edge_index_table(adj_list), ef, nf)

        @test size(batch_message) == (out_channel_v + out_channel_e, 2)
    end

    @testset "update_batch_edge()" begin
        ft_layer = SeaPearl.EdgeFtLayer(;v_dim = in_channel_v=>out_channel_v, e_dim = in_channel_e=>out_channel_e)
        nf = rand(Float32, in_channel_v, N)
        ef = rand(Float32, in_channel_e, 8)
        fg = GraphSignals.FeaturedGraph(adj; nf=nf, ef=ef)

        adj_list = GeometricFlux.adjacency_list(fg)

        every_messages = GeometricFlux.update_batch_edge(ft_layer, adj_list, ef, nf)

        @test size(every_messages) == (out_channel_v + out_channel_e, 8)
    end
end


text:
using Test
using FunctionIndices
using OffsetArrays

struct YANI{T} <: FunctionIndices.AbstractNotIndex{T}
    parent::T
end
Base.parent(I::YANI) = I.parent
YANI(x::T, xs::T...) where {T} = YANI((x, xs...))

@testset "FunctionIndices" begin
    A = collect(reshape(1:12, 3, 4))
    OA = OffsetArray(A, 2:4, 0:3)

    @testset "FI" begin
        @test A[FI(iseven), FI(iseven)] == A[[2], [2, 4]]
        @test OA[FI(iseven), FI(iseven)] == OA[[2, 4], [0, 2]] == A[[1, 3], [1, 3]]
        AdjA = A'
        @test AdjA[FI(iseven), FI(iseven)] ==
              AdjA[[2, 4], [2]] ==
              A[FI(iseven), FI(iseven)]'
    end

    @testset "$_not" for _not in (not, YANI)
        @testset "not base" begin
            @test ones()[_not(CartesianIndex())] == 1
            @test A[_not(1), _not(2)] ==
                A[axes(A, 1).!=1, axes(A, 2).!=2] ==
                A[_not(CartesianIndex(1, 2))] ==
                A[_not(1), _not(2)] ==
                A[_not(CartesianIndex(1)), _not(2)] ==
                A[_not(CartesianIndex()), _not(1), _not(2)]
            @test isempty(A[_not(1), _not(2), _not(1)])
            @test !isempty(A[_not(1), _not(2), _not(2)])
            I = Bool[0, 1, 0]
            J = Bool[1, 0, 1, 0]
            @test A[_not(I), _not(J)] ==
                A[.!I, .!J] ==
                A[_not(Base.to_index(I)), _not(Base.to_index(J))]
            @test OA[_not(1), _not(2)] == OA[axes(OA, 1).!=1, axes(OA, 2).!=2] == A[:, _not(3)]
            @test A[_not(2)] ==
                A[_not(CartesianIndex(2))] ==
                A[1:length(A).!=2] ==
                OA[_not(2)] ==
                OA[_not(CartesianIndex(2))] ==
                OA[1:length(A).!=2]
            @test A[_not(2, 3)] == A[_not(2:3)] == OA[_not(2, 3)] == OA[_not(2:3)]
            @test A[_not(2, 3), _not(1, 3)] ==
                A[_not(2:3), _not(1:2:3)] ==
                OA[_not(3, 4), _not(0, 2)] ==
                OA[_not(3:4), _not(0:2:2)]
            @test A[_not(CartesianIndices((1:2, 2:3)))] ==
                A[_not(1:2), _not(2:3)] ==
                OA[_not(CartesianIndices((2:3, 1:2)))] ==
                OA[_not(2:3), _not(1:2)]
            # adjoint
            AdjA = A'
            @test AdjA[_not(1), _not(2)] ==
                AdjA[axes(AdjA, 1).!=1, axes(AdjA, 2).!=2] ==
                A[_not(2), _not(1)]'
            # slice and colon
            slc1, slc2 = Base.to_indices(A, (:, :)) # slices
            @test A[_not(:), :] == A[_not(slc1), :] == A[[], :] == similar(A, 0, size(A, 2))
            @test A[:, _not(:)] == A[:, _not(slc2)] == A[:, []] == similar(A, size(A, 1), 0)
            @test A[_not(:), _not(:)] == A[_not(slc1), _not(slc2)] == A[[], []] == similar(A, 0, 0)
            # slice is _not equal to axe
            @test OA[_not(slc1), :] == OA[[4], :]
            @test OA[:, _not(slc2)] == OA[:, [0]]
            @test OA[_not(slc1), _not(slc2)] == OA[[4], [0]]
            # converted index
            I, J = to_indices(A, (_not([2, 3]), FI(iseven)))
            K, L = to_indices(A, (_not(1, 2), FI(isodd)))
            @test OA[_not(I), _not(J)] ==
                OA[[2, 3], FI(isodd)] ==
                A[_not(K), _not(L)] ==
                A[[1, 2], FI(iseven)]
        end

        @testset "_not for outbound index" begin
            @test A[_not(0)] == A[_not(20)] == OA[_not(0)] == OA[_not(20)] == vec(A) == vec(OA)
            @test A[_not(0:13)] == A[_not(1:12)] == OA[_not(0:13)] == OA[_not(1:12)]
            @test A[_not(-1:1)] == A[_not(1:1)] == OA[_not(-1:1)] == OA[_not(1:1)]
            @test A[_not(0:2), _not(0)] == A[_not(1:2), _not(5)] == OA[_not(1:3), _not(-1)] == OA[_not(1:3), _not(4)] == A[3:3, 1:4] == OA[4:4, 0:3]
        end

        @testset "not for StepRange" begin
            @test A[not(1:0)] == A[not(0:-1:1)] == OA[not(1:0)] == OA[not(0:-1:1)] == vec(A) == vec(OA)
            @test A[not(2:10)] == OA[not(2:10)] == A[not(2:1:10)] == OA[not(2:1:10)] == A[[1, 11, 12]] == OA[[1, 11, 12]]
            @test A[not(10:-1:2)] == OA[not(10:-1:2)] == A[[1, 11, 12]] == OA[[1, 11, 12]]
            @test A[not(4:2:10)] == OA[not(4:2:10)] == A[[1:3; 5:2:9; 11:12]] == OA[[1:3; 5:2:9; 11:12]]
            @test A[not(10:-2:4)] == OA[not(10:-2:4)] == A[[1:3; 5:2:9; 11:12]] == OA[[1:3; 5:2:9; 11:12]]
            @test A[not(1:2:3), not(1:3:4)] == OA[not(2:2:4), not(0:3:3)] == A[[2], 2:3] == OA[[3], 1:2]
            @test A[not(1:2:3), not(4:-3:1)] == OA[not(2:2:4), not(3:-3:0)] == A[[2], 2:3] == OA[[3], 1:2]
        end
    end

    @testset "enable/disable optimization" begin
        @test to_indices(A, (not(0:2:2),))[1] isa Vector{Int}
        # disable optimization for not(::StepRange)
        FunctionIndices.indextype(
            ::Type{<:AbstractArray},
            ::Type{<:NotIndex{<:StepRange}},
        ) = Base.LogicalIndex
        @test to_indices(A, (not(0:2:2),))[1] isa Base.LogicalIndex
    end
end


text:
# Faster version of permga.jl

module Spmap

using SparseArrays
import SparseArrays: indtype, getcolptr
using Random

@inline colstartind(A::SparseMatrixCSC, j) = getcolptr(A)[j]
@inline colboundind(A::SparseMatrixCSC, j) = getcolptr(A)[j + 1]
@inline setcolptr!(A::SparseMatrixCSC, j, val) = getcolptr(A)[j] = val

function trimstorage!(A::SparseMatrixCSC, maxstored)
    resize!(rowvals(A), maxstored)
    resize!(nonzeros(A), maxstored)
    return maxstored
end
function expandstorage!(A::SparseMatrixCSC, maxstored)
    length(rowvals(A)) < maxstored && resize!(rowvals(A), maxstored)
    length(nonzeros(A)) < maxstored && resize!(nonzeros(A), maxstored)
    return maxstored
end

# map!(+, C, A, B[h,h][1:size(A,1),1:size(A,1)]), where size(C) == size(A)
# This version assumes that all values in A and B are the value 1
function map_index_plus!(C::SparseMatrixCSC, A::SparseMatrixCSC, B::SparseMatrixCSC,
    h::AbstractVector, hinv::AbstractVector, hinvrowvals::AbstractVector)

    spaceC::Int = min(length(rowvals(C)), length(nonzeros(C)))
    rowsentinelA = convert(indtype(A), size(C,1) + 1)
    rowsentinelB = convert(indtype(B), size(C,1) + 1)

    (length(hinvrowvals) < length(rowvals(B))) && error("length hinvrowvals")
    @inbounds for hi in 1:length(rowvals(B))
        hinvrowvals[hi] = min(rowsentinelB, hinv[rowvals(B)[hi]])
    end
    @inbounds for j in 1:size(B,2)
        Bk, stopBk = colstartind(B, j), colboundind(B, j)
        sort!(hinvrowvals, Bk, stopBk-1, QuickSort, Base.Order.Forward)
    end

    Ck = 1
    @inbounds for j in 1:size(C,2)
        bj = h[j]
        setcolptr!(C, j, Ck)
        Ak, stopAk = colstartind(A, j), colboundind(A, j)
        Bk, stopBk = colstartind(B, bj), colboundind(B, bj)
        Ai = Ak < stopAk ? rowvals(A)[Ak] : rowsentinelA

        #Bi = Bk < stopBk ? hinv[rowvals(B)[Bk]] : rowsentinelB
        Bi = Bk < stopBk ? hinvrowvals[Bk] : rowsentinelB
        while true
            #println("j $j bj $bj Ak $Ak stopAk $stopAk Ai $Ai Bk $Bk stopBk $stopBk Bi $Bi Ck $Ck")
            if Ai == Bi
                Ai == rowsentinelA && break # column complete
                Cx, Ci::indtype(C) = 2, Ai
                Ak += oneunit(Ak); Ai = Ak < stopAk ? rowvals(A)[Ak] : rowsentinelA
                #Bk += oneunit(Bk); Bi = Bk < stopBk ? hinv[rowvals(B)[Bk]] : rowsentinelB
                Bk += oneunit(Bk); Bi = Bk < stopBk ? hinvrowvals[Bk] : rowsentinelB
            elseif Ai < Bi
                Cx, Ci = 1, Ai
                Ak += oneunit(Ak); Ai = Ak < stopAk ? rowvals(A)[Ak] : rowsentinelA
            else # Bi < Ai
                Cx, Ci = 1, Bi
                #Bk += oneunit(Bk); Bi = Bk < stopBk ? hinv[rowvals(B)[Bk]] : rowsentinelB
                Bk += oneunit(Bk); Bi = Bk < stopBk ? hinvrowvals[Bk] : rowsentinelB
            end
            # NOTE: The ordering of the conditional chain above impacts which matrices this
            # method performs best for. In the map situation (arguments have same shape, and
            # likely same or similar stored entry pattern), the Ai == Bi and termination
            # cases are equally or more likely than the Ai < Bi and Bi < Ai cases.
            if !iszero(Cx)
                if Ck > spaceC
                    error("Ck > spaceC")
                    spaceC = expandstorage!(C, Ck + (nnz(A) - (Ak - 1)) + (nnz(B) - (Bk - 1)))
                end
                rowvals(C)[Ck] = Ci
                nonzeros(C)[Ck] = Cx
                Ck += 1
            end
        end
    end
    @inbounds setcolptr!(C, size(C,2) + 1, Ck)
    # trimstorage!(C, Ck - 1)
    return C
end

function invperm!(b::AbstractVector, a::AbstractVector)
    # b = zero(a) # similar vector of zeros
    b .= 0
    n = length(a)
    (length(b) == n) || throw(ArgumentError("vector lengths"))
    @inbounds for (i, j) in enumerate(a)
        ((1 <= j <= n) && b[j] == 0) ||
            throw(ArgumentError("argument is not a permutation"))
        b[j] = i
    end
    b
end

# Auxiliary space to reduce allocations
function map_index_aux(A::SparseMatrixCSC, B::SparseMatrixCSC)
    S = size(A)
    IT = promote_type(SparseArrays.indtype(A), SparseArrays.indtype(B))
    ET = promote_type(eltype(A), eltype(B))

    maxnnz = nnz(A)+nnz(B)
    pointers = ones(IT, S[2] + 1)
    storedinds = Vector{IT}(undef, maxnnz)
    storedvals = Vector{ET}(undef, maxnnz)
    C = SparseMatrixCSC(S..., pointers, storedinds, storedvals)

    hinv = randperm(size(B,1))
    hinvrowvals = similar(hinv, nnz(B))

    (C, hinv, hinvrowvals)
end

function plus_getindex!(C::SparseMatrixCSC, A::SparseMatrixCSC, B::SparseMatrixCSC,
    h::AbstractVector, hinv::AbstractVector, hinvrowvals::AbstractVector)
    invperm!(hinv, h)
    rowvals(C) .= 0
    nonzeros(C) .= 0
    map_index_plus!(C, A, B, h, hinv, hinvrowvals)
end

end

module NetalGA

using SparseArrays
using Random
using ..GAFramework
import ..GAFramework: fitness, crossover!, mutation!, selection, randcreature, genauxga, printfitness
using ..CayleyCrossover
using ..Spmap

export NetalCreature, NetalModel

# Represents a permutation, i.e. 1-1 mapping from one set to another
struct NetalCreature
    f :: Vector{Int}
    score :: Float64 # alignment score
end

fitness(x::NetalCreature) where {T} = x.score

printfitness(curgen::Int, x::NetalCreature) =
    println("curgen: $curgen value: $(x.f) score: $(x.score)")

"""
Model to find network alignment by optimizing S3 (see MAGNA paper)
A network alignment is represented using a permutation
Faster implementation, slightly slower than the C++ code from the MAGNA++ paper
    using GAFramework, GAFramework.NetalGA
    using LightGraphs
    using Random
    mv = 100
    me = 400
    G1 = adjacency_matrix(LightGraphs.erdos_renyi(mv,me))
    perm = randperm(mv)
    G2 = G1[invperm(perm),invperm(perm)]
    # G1.|V| <= G2.|V|
    model = NetalModel(G1,G2)
    st = GAState(model, ngen=100, npop=6_000, elite_fraction=0.1,
        params=Dict(:print_fitness_iter=>1))
    ga!(st)
    println("accuracy = ", sum(st.pop[1].f .== perm)/mv)
    st.pop[1]
"""
struct NetalModel <: GAModel
    G1   :: SparseMatrixCSC{Int,Int}
    G2   :: SparseMatrixCSC{Int,Int}
end

"""
    using GAFramework, GAFramework.NetalGA, GAFramework.PermGA
    using LightGraphs, SparseArrays
    using Random
    GAFramework.setthreads(false)
    mv = 100
    me = 400
    G2 = adjacency_matrix(LightGraphs.erdos_renyi(mv,me))
    perm = randperm(mv)
    nv = 75
    G1 = G2[perm[1:nv], perm[1:nv]]
    model1 = NetalModel(G1,G2)
    model2 = NetalignModel(G1,G2)
    aux1 = genauxga(model1)
    f = randperm(size(G2,1));
    c1 = NetalCreature(f, model1, aux1)
    c2 = PermCreature(f, model2)
    c1.score, c2.score

    st1 = GAState(model1, ngen=10, npop=60_000, elite_fraction=0.1, rng=MersenneTwister(10));
    st2 = GAState(model2, ngen=10, npop=60_000, elite_fraction=0.1, rng=MersenneTwister(10));

    @time ga!(st1)
    @time ga!(st2)

    G2 = sparse([1, 1, 2, 3, 4], [2, 3, 3, 4, 5], ones(Int,5), 5, 5)
    G2 = min.(1, (G2 + G2'))
    G1 = G2[1:4,1:4]
    model1 = NetalModel(G1,G2)
    model2 = NetalignModel(G1,G2)

    aux1 = genauxga(model1)
    f = randperm(size(G2,1))
    c1 = NetalCreature(f, model1, aux1)
    c2 = PermCreature(f, model2)
    c1.score, c2.score

    st = GAState(model1, ngen=100, npop=6_000, elite_fraction=0.1,
        params=Dict(:print_fitness_iter=>1))
    ga!(st)
    println("accuracy = ", sum(st.pop[1].f[1:4] .== [1,2,3,4])/4)
    println("accuracy = ", sum(st.pop[1].f[1:4] .== [2,1,3,4])/4)

    # using Profile
    # @profile NetalCreature(f, model1, aux1)
    # open("/home/vvijayan/home/gaprof.txt", "w") do fd Profile.print(fd) end

"""

function genauxga(m::NetalModel)
    crossover_aux = randperm(size(m.G2,1)), fill(false, size(m.G2,1))
    Spmap.map_index_aux(m.G1, m.G2), crossover_aux
end

# aux = genauxga(..) can be used to reduce allocations here
# by using auxiliary space instead of reallocating
function NetalCreature(h, m::NetalModel, aux)
    # Assumes that all edge weights are 1
    ((K, hinv, hinvrowvals), _) = aux
    Spmap.plus_getindex!(K, m.G1, m.G2, h, hinv, hinvrowvals) # K = G1 + G2[h,h][1:m,1:m]

    w = nonzeros(K)
    Nc = count(x->x==2, w)
    Nn = count(x->x==1, w)

    Nc,cr = divrem(Nc,2)
    Nn,nr = divrem(Nn,2) # need this & above so it works well with sim. ann. code
    if cr!=0 || nr!=0 error("G1 and G2 need to be symmetric"); end
    score = Nc/(Nc + Nn)
    NetalCreature(h, score)
end

function randcreature(m::NetalModel, aux, rng::AbstractRNG) where {T}
    f = randperm(rng, size(m.G2,1))
    NetalCreature(f, m, aux)
end

function crossover!(z, x, y, st::GAState{NetalModel}, aux, rng::AbstractRNG) where {T}
    (_, (r, visited)) = aux
    visited .= false
    cayley_crossover!(z.f, x.f, y.f, rng, r, visited)
    NetalCreature(z.f, st.model, aux)
end

end # NetalGA


text:
using Test
using BasisFunctions, GridArrays

@testset "Basic PiecewiseConstants tests" begin
    s = zeros(10);s[1]=1
    H = PiecewiseConstants(10)
    @test BasisFunctions.dict_eval(H,0)≈s
    g = interpolation_grid(H)
    @test g≈MidpointEquispacedGrid(10,0,1)
    @test isorthogonal(H,discretemeasure(g))
    @test isorthogonal(H,FourierWeight())
    @test isorthogonal(H,discretemeasure(MidpointEquispacedGrid(10,0,1)))

    @test support(H) ≈ 0.0..1.0
    @test support(H,1) ≈ 0.0..0.1

    @test 0 ∈ support(H)
    @test 1 ∉ support(H)

    @test 0 ∈ support(H,1)
    @test .1 ∉ support(H,1)

    @test gram(H,discretemeasure(g)) ≈ IdentityOperator(H)
    @test evaluation(H,g)≈IdentityOperator(H)

    @test gram(H,FourierWeight())≈.1IdentityOperator(H)
end


text:
"""
Support for nondeterministic programming using iterable, backtracking Choices
and a special method of the |> operator for combining them with (anonymous)
functions (which stand in for continuations). This is as close as we can get to
McCarthy's amb operator in Julia, given that it doesn't have call/cc or
anything like it.
"""
module Choices

export choose, fail, cut

"""
    Choice{T}

Represents a call to choose() whose options (of type T) we have not yet fully
explored, as well as the overall nondeterministic computation that this choice
is at the end of, and the collection of its successful results.

# Example
```jldoctest
get_first_part() = choose("razzle", "root")
get_second_part() = choose("beer", "dazzle")
function check_pair(first_part, second_part)
  length(first_part) == length(second_part) || fail()
  return (first_part, second_part)
end
get_both_parts() =
  get_first_part() |> first_part ->
  get_second_part() |> second_part ->
  check_pair(first_part, second_part)
[get_both_parts()...]

# output

2-element Array{Tuple{String,String},1}:
 ("razzle", "dazzle")
 ("root", "beer")
```

See also [`choose`](@ref), [`|>`](@ref), [`fail`](@ref), and [`cut`](@ref).
"""
struct Choice{O,R}
  # an iterable collection of the options for this choice
  options::O
  # the Choice we should backtrack to (if any) if we run out of options
  back::Union{Choice,Nothing}
  # the function (or at least something callable) to pass the next option to
  receive::R
  # when we create a Choice, store it in task-local storage
  function Choice(options::O, back::Union{Choice,Nothing}, receive::R) where {O, R}
    c = new{O,R}(options, back, receive)
    task_local_storage(:previous_choice, c)
    return c
  end
end

function Base.show(io::IO, x::Choice)
  print(io, "(choose from ", x.options)
  identity === x.receive ||
    print(io, " and call ",
          (occursin(r"^\w+$", string(nameof(x.receive))) ?
	    x.receive : "a nontrivial function"))
  nothing === x.back || print(io, " or go back to ", x.back)
  print(io, ")")
end

"""
    choose(options[, back::Choice])

Make a Choice among the given options, backtracking to the given Choice if
those options run out. If back isn't given, backtrack to the previous Choice
made in this Task, if any.

This form accepts one iterator of options, and it will only be iterated as the
options are explored, so it's possible to use unbounded iterators.

See also [`Choice`](@ref).
"""
choose(options,
       back::Union{Choice,Nothing} =
	 get(task_local_storage(), :previous_choice, nothing)) =
  Choice(options, back, identity)

"""
    choose(options...)

Make a Choice among the given options, backtracking to the previous Choice made
in this Task, if any.

In this form, the options are given individually in the call, and there must be
at least two of them.

See also [`Choice`](@ref).
"""
choose(first_option, second_option, rest_of_options...) =
  choose((first_option, second_option, rest_of_options...))

"""
    choice |> function

Chain the given function onto the successful results of the given Choice. This
is slightly different from the usual meaning of |> in Julia, since the function
receives a result of the Choice, not the Choice itself, and the value of the
whole expression is another Choice, not the result of calling the function.

See also [`Choice`](@ref).
"""
function Base.:|>(x::Choice, f::Function)
  old_receive = x.receive # capture just the old receive function, not x
  return Choice(x.options, x.back, input -> (old_receive(input) |> f))
end

struct Failed <: Exception
end

"""
    fail()

Throw a Choices.Failed exception to indicate that the option currently being
considered won't work. Triggers backtracking while there are other options to
consider.

See also [`Choice`](@ref).
"""
fail() = throw(Failed())

"""
    cut()

Commit to any choices already made, preventing backtracking past this point.

See also [`Choice`](@ref).
"""
cut() = task_local_storage(:previous_choice, nothing)

# TODO mark(), make it so that a subsequent cut()-fail() sequence backtracks back to the marked point (skipping the region of the computation between mark() and cut())

"""
    Choices.IterationState

Track the state of iteration over successful results of a nondeterministic
computation.

See also [`Choice`](@ref).
"""
struct IterationState
  # the iteration state we should backtrack to (if any) if we run out of options
  back::Union{IterationState,Nothing}
  # the choice whose options we're currently iterating over
  currentChoice::Choice
  # the iteration state of that choice's options
  optionsState
end

Base.IteratorSize(::Type{Choice{T}}) where {T} = Base.SizeUnknown()
Base.IteratorEltype(::Type{Choice{T}}) where {T} = Base.EltypeUnknown()

Base.iterate(iter::Choice) =
  iterate(iter, IterationState(nothing, iter, nothing))

"""
    iterate(iter::Choice[, state::IterationState])

Iterate over successful results of the nondeterministic computation represented by x.

See also [`Choice`](@ref).
"""
function Base.iterate(x::Choice, state::IterationState)
  # basically this is "for opt ∈ state.currentChoice.options", except we expose
  # the iteration state of the options iterator as optionsState, and resume
  # from state.optionsState
  currentChoice = state.currentChoice
  optionsState = state.optionsState
  iterateResult =
    (nothing === optionsState ? iterate(currentChoice.options) :
	iterate(currentChoice.options, optionsState))
  while (nothing !== iterateResult)
    opt, optionsState = iterateResult
    @debug "considering option" opt
    try
      # try executing the receive function on the option
      result = currentChoice.receive(opt)
      @debug "got result" result
      # it succeeded!
      # prepare the next iteration state on this choice
      nextStateHere = IterationState(state.back, currentChoice, optionsState)
      if (isa(result, Choice)) # if the function returned another Choice, recur
	@debug "recursing..."
	return iterate(x, IterationState(nextStateHere, result, nothing))
      else # otherwise, just return the result with the next state
	@debug "returning!"
	return (result, nextStateHere)
      end
    catch ex
      if (!isa(ex, Failed))
	rethrow(ex)
      end # else move on to the next option for this Choice
    end
    @debug "iterating to next option"
    iterateResult =
      (nothing === optionsState ? iterate(currentChoice.options) :
	  iterate(currentChoice.options, optionsState))
  end
  # if we get here, we have exhausted the options for this Choice
  @debug "exhausted options" currentChoice.options
  if (nothing === state.back) # we can't backtrack anymore
    @debug "exhausted all options!"
    return nothing
  else # backtrack
    @debug "backtracking..."
    return iterate(x, state.back)
  end
end

end # module