text:
function parse_mask(mask::AbstractString)
    ones = parse(Int, replace(mask, 'X' => '0'), base=2)
    zeros = parse(Int,
    replace(replace(replace(mask, '1' => 'X'), '0' => '1'), 'X' => '0'),
    base=2)
    return ones, zeros
end

function apply_mask(val::Int, ones::Int, zeros::Int)
    return ones | val & ~zeros
end

function part1(lines)
    ones = zeros = 0
    mem = Dict{Int, Int}()
    for line in lines
        if startswith(line, "ma")
            # 8 because the line is "mask = XXXX...", so X starts at 8
            ones, zeros = parse_mask(line[8:end])
        else
            place, value = parse.(Int, (m.match for m in eachmatch(r"\d+", line)))
            mem[place] = apply_mask(value, ones, zeros)
        end
    end
    return sum(values(mem))
end

const DATA = readlines("day14.txt")

println("Part 1: $(part1(DATA))")

using Combinatorics

function parse_mask_p2(mask::AbstractString)
    ones = parse(Int, replace(mask, 'X' => '0'), base=2)
    Xindices = reduce(vcat, findall("X", reverse(mask)))
    # 1 << n is 2^n, and as the indexes start at 1, we need to substract 1
    Xvalues = 1 .<< (Xindices .- 1)
    return ones, Xvalues
end

function apply_mask_p2(val::Int, ones::Int, Xs::Array{Int})
    values = Int[]
    # set all the mask's 1 to 1 and set all the Xs to 0
    start_val = (val | ones) & ~(sum(Xs))
    push!(values, start_val)
    for comb in combinations(Xs)
        mask = sum(comb)
        push!(values, start_val | mask)
    end
    return values
end

function part2(lines)
    ones = 0
    Xs = Int[]
    mem = Dict{Int,Int}()
    for line in lines
        if startswith(line, "ma")
            # 8 because the line is "mask = XXXX...", so X starts at 8
            ones, Xs = parse_mask_p2(line[8:end])
        else
            place, value = parse.(Int, (m.match for m in eachmatch(r"\d+", line)))
            new_adresses = apply_mask_p2(place, ones, Xs)
            for address in new_adresses
                mem[address] = value
            end
        end
    end
    return sum(values(mem))
end

println("Part 2: $(part2(DATA))")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RoundRotorMachine <: Machine
        R::Float64
        Td0_p::Float64
        Td0_pp::Float64
        Tq0_p::Float64
        Tq0_pp::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Xd_pp::Float64
        Xl::Float64
        Se::Tuple{Float64, Float64}
        ext::Dict{String, Any}
        γ_d1::Float64
        γ_q1::Float64
        γ_d2::Float64
        γ_q2::Float64
        γ_qd::Float64
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of 4-states round-rotor synchronous machine with quadratic/exponential saturation:
IEEE Std 1110 §5.3.2 (Model 2.2). GENROU or GENROE model in PSSE and PSLF.

# Arguments
- `R::Float64`: Armature resistance, validation range: `(0, nothing)`
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: `(0, nothing)`
- `Td0_pp::Float64`: Time constant of sub-transient d-axis voltage, validation range: `(0, nothing)`
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: `(0, nothing)`
- `Tq0_pp::Float64`: Time constant of sub-transient q-axis voltage, validation range: `(0, nothing)`
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Xd_pp::Float64`: Sub-Transient reactance after EMF in d-axis per unit. Note: Xd_pp = Xq_pp, validation range: `(0, nothing)`
- `Xl::Float64`: Stator leakage reactance, validation range: `(0, nothing)`
- `Se::Tuple{Float64, Float64}`: Saturation factor at 1 and 1.2 pu flux: S(1.0) = B(|ψ_pp|-A)^2
- `ext::Dict{String, Any}`
- `γ_d1::Float64`: γ_d1 parameter
- `γ_q1::Float64`: γ_q1 parameter
- `γ_d2::Float64`: γ_d2 parameter
- `γ_q2::Float64`: γ_q2 parameter
- `γ_qd::Float64`: γ_qd parameter
- `states::Vector{Symbol}`: The states are:
	eq_p: q-axis generator voltage behind the transient reactance,
	ed_p: d-axis generator voltage behind the transient reactance,
	ψ_kd: flux linkage in the first equivalent damping circuit in the d-axis,
	ψ_kq: flux linkage in the first equivalent damping circuit in the d-axis
- `n_states::Int`: RoundRotorMachine has 4 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct RoundRotorMachine <: Machine
    "Armature resistance"
    R::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of sub-transient d-axis voltage"
    Td0_pp::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    "Time constant of sub-transient q-axis voltage"
    Tq0_pp::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Sub-Transient reactance after EMF in d-axis per unit. Note: Xd_pp = Xq_pp"
    Xd_pp::Float64
    "Stator leakage reactance"
    Xl::Float64
    "Saturation factor at 1 and 1.2 pu flux: S(1.0) = B(|ψ_pp|-A)^2"
    Se::Tuple{Float64, Float64}
    ext::Dict{String, Any}
    "γ_d1 parameter"
    γ_d1::Float64
    "γ_q1 parameter"
    γ_q1::Float64
    "γ_d2 parameter"
    γ_d2::Float64
    "γ_q2 parameter"
    γ_q2::Float64
    "γ_qd parameter"
    γ_qd::Float64
    "The states are:
	eq_p: q-axis generator voltage behind the transient reactance,
	ed_p: d-axis generator voltage behind the transient reactance,
	ψ_kd: flux linkage in the first equivalent damping circuit in the d-axis,
	ψ_kq: flux linkage in the first equivalent damping circuit in the d-axis"
    states::Vector{Symbol}
    "RoundRotorMachine has 4 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function RoundRotorMachine(R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext=Dict{String, Any}(), )
    RoundRotorMachine(R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext, (Xd_pp - Xl) / (Xd_p - Xl), (Xd_pp - Xl) / (Xq_p - Xl), (Xd_p - Xd_pp) / (Xd_p - Xl)^2, (Xq_p - Xd_pp) / (Xq_p - Xl)^2, (Xq - Xl) / (Xd - Xl), [:eq_p, :ed_p, :ψ_kd, :ψ_kq], 4, InfrastructureSystemsInternal(), )
end

function RoundRotorMachine(; R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext=Dict{String, Any}(), γ_d1=(Xd_pp - Xl) / (Xd_p - Xl), γ_q1=(Xd_pp - Xl) / (Xq_p - Xl), γ_d2=(Xd_p - Xd_pp) / (Xd_p - Xl)^2, γ_q2=(Xq_p - Xd_pp) / (Xq_p - Xl)^2, γ_qd=(Xq - Xl) / (Xd - Xl), states=[:eq_p, :ed_p, :ψ_kd, :ψ_kq], n_states=4, internal=InfrastructureSystemsInternal(), )
    RoundRotorMachine(R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext, γ_d1, γ_q1, γ_d2, γ_q2, γ_qd, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function RoundRotorMachine(::Nothing)
    RoundRotorMachine(;
        R=0,
        Td0_p=0,
        Td0_pp=0,
        Tq0_p=0,
        Tq0_pp=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Xd_pp=0,
        Xl=0,
        Se=(0.0, 0.0),
        ext=Dict{String, Any}(),
    )
end

"""Get [`RoundRotorMachine`](@ref) `R`."""
get_R(value::RoundRotorMachine) = value.R
"""Get [`RoundRotorMachine`](@ref) `Td0_p`."""
get_Td0_p(value::RoundRotorMachine) = value.Td0_p
"""Get [`RoundRotorMachine`](@ref) `Td0_pp`."""
get_Td0_pp(value::RoundRotorMachine) = value.Td0_pp
"""Get [`RoundRotorMachine`](@ref) `Tq0_p`."""
get_Tq0_p(value::RoundRotorMachine) = value.Tq0_p
"""Get [`RoundRotorMachine`](@ref) `Tq0_pp`."""
get_Tq0_pp(value::RoundRotorMachine) = value.Tq0_pp
"""Get [`RoundRotorMachine`](@ref) `Xd`."""
get_Xd(value::RoundRotorMachine) = value.Xd
"""Get [`RoundRotorMachine`](@ref) `Xq`."""
get_Xq(value::RoundRotorMachine) = value.Xq
"""Get [`RoundRotorMachine`](@ref) `Xd_p`."""
get_Xd_p(value::RoundRotorMachine) = value.Xd_p
"""Get [`RoundRotorMachine`](@ref) `Xq_p`."""
get_Xq_p(value::RoundRotorMachine) = value.Xq_p
"""Get [`RoundRotorMachine`](@ref) `Xd_pp`."""
get_Xd_pp(value::RoundRotorMachine) = value.Xd_pp
"""Get [`RoundRotorMachine`](@ref) `Xl`."""
get_Xl(value::RoundRotorMachine) = value.Xl
"""Get [`RoundRotorMachine`](@ref) `Se`."""
get_Se(value::RoundRotorMachine) = value.Se
"""Get [`RoundRotorMachine`](@ref) `ext`."""
get_ext(value::RoundRotorMachine) = value.ext
"""Get [`RoundRotorMachine`](@ref) `γ_d1`."""
get_γ_d1(value::RoundRotorMachine) = value.γ_d1
"""Get [`RoundRotorMachine`](@ref) `γ_q1`."""
get_γ_q1(value::RoundRotorMachine) = value.γ_q1
"""Get [`RoundRotorMachine`](@ref) `γ_d2`."""
get_γ_d2(value::RoundRotorMachine) = value.γ_d2
"""Get [`RoundRotorMachine`](@ref) `γ_q2`."""
get_γ_q2(value::RoundRotorMachine) = value.γ_q2
"""Get [`RoundRotorMachine`](@ref) `γ_qd`."""
get_γ_qd(value::RoundRotorMachine) = value.γ_qd
"""Get [`RoundRotorMachine`](@ref) `states`."""
get_states(value::RoundRotorMachine) = value.states
"""Get [`RoundRotorMachine`](@ref) `n_states`."""
get_n_states(value::RoundRotorMachine) = value.n_states
"""Get [`RoundRotorMachine`](@ref) `internal`."""
get_internal(value::RoundRotorMachine) = value.internal

"""Set [`RoundRotorMachine`](@ref) `R`."""
set_R!(value::RoundRotorMachine, val) = value.R = val
"""Set [`RoundRotorMachine`](@ref) `Td0_p`."""
set_Td0_p!(value::RoundRotorMachine, val) = value.Td0_p = val
"""Set [`RoundRotorMachine`](@ref) `Td0_pp`."""
set_Td0_pp!(value::RoundRotorMachine, val) = value.Td0_pp = val
"""Set [`RoundRotorMachine`](@ref) `Tq0_p`."""
set_Tq0_p!(value::RoundRotorMachine, val) = value.Tq0_p = val
"""Set [`RoundRotorMachine`](@ref) `Tq0_pp`."""
set_Tq0_pp!(value::RoundRotorMachine, val) = value.Tq0_pp = val
"""Set [`RoundRotorMachine`](@ref) `Xd`."""
set_Xd!(value::RoundRotorMachine, val) = value.Xd = val
"""Set [`RoundRotorMachine`](@ref) `Xq`."""
set_Xq!(value::RoundRotorMachine, val) = value.Xq = val
"""Set [`RoundRotorMachine`](@ref) `Xd_p`."""
set_Xd_p!(value::RoundRotorMachine, val) = value.Xd_p = val
"""Set [`RoundRotorMachine`](@ref) `Xq_p`."""
set_Xq_p!(value::RoundRotorMachine, val) = value.Xq_p = val
"""Set [`RoundRotorMachine`](@ref) `Xd_pp`."""
set_Xd_pp!(value::RoundRotorMachine, val) = value.Xd_pp = val
"""Set [`RoundRotorMachine`](@ref) `Xl`."""
set_Xl!(value::RoundRotorMachine, val) = value.Xl = val
"""Set [`RoundRotorMachine`](@ref) `Se`."""
set_Se!(value::RoundRotorMachine, val) = value.Se = val
"""Set [`RoundRotorMachine`](@ref) `ext`."""
set_ext!(value::RoundRotorMachine, val) = value.ext = val
"""Set [`RoundRotorMachine`](@ref) `γ_d1`."""
set_γ_d1!(value::RoundRotorMachine, val) = value.γ_d1 = val
"""Set [`RoundRotorMachine`](@ref) `γ_q1`."""
set_γ_q1!(value::RoundRotorMachine, val) = value.γ_q1 = val
"""Set [`RoundRotorMachine`](@ref) `γ_d2`."""
set_γ_d2!(value::RoundRotorMachine, val) = value.γ_d2 = val
"""Set [`RoundRotorMachine`](@ref) `γ_q2`."""
set_γ_q2!(value::RoundRotorMachine, val) = value.γ_q2 = val
"""Set [`RoundRotorMachine`](@ref) `γ_qd`."""
set_γ_qd!(value::RoundRotorMachine, val) = value.γ_qd = val

ext: jl


text:
"""
Support for nng library
"""
mutable struct nng_socket
    id::UInt32
end

mutable struct nng_listener
    id::UInt32
end

mutable struct nng_dialer
    id::UInt32
end

function nng_strerror(errno::Int32)
    err = ccall((:nn_strerror, "libnng"), Cstring, (Int32,), errno)
    if err == C_NULL
        error("Error nn_strerror errno: ", errno)
    end
    return unsafe_string(err)
end

function check_err(func_name, rv)
    if rv != 0
        err = nng_strerror(rv)
        error(func_name, err)
    end
end

nng_pull0_open(socket::nng_socket) = ccall((:nng_pull0_open, "libnng"), Int32, (Ref{nng_socket},), socket)
nng_push0_open(socket::nng_socket) = ccall((:nng_push0_open, "libnng"), Int32, (Ref{nng_socket},), socket)
nng_close(socket::nng_socket) = ccall((:nng_close, "libnng"), Int32, (nng_socket,), socket)

function nng_listen(sock::nng_socket, url, lp;  flags=0)
    return ccall((:nng_listen, "libnng"), Cint, (nng_socket, Ptr{UInt8}, Ref{nng_listener}, Cint), sock, url, lp, flags)
end

function nng_dial(socket::nng_socket, url, dp, flags)
    return ccall((:nng_dial, "libnng"), Cint, (nng_socket,  Ptr{UInt8}, Ref{nng_dialer}, Cint), socket, url, dp, flags)
end

function nng_send(socket::nng_socket, msg, flags)
    sz_msg = length(msg)
    bytes = ccall((:nng_send, "libnng"), Int32, (nng_socket, Ptr{UInt8}, Csize_t, Int32), socket, msg, sz_msg, flags)
    return bytes
end

function nng_recv(socket::nng_socket; blocking=false)
    flags = blocking ? 2 : 1 # NNG_FLAG_ALLOC=1 NNG_FLAG_NONBLOCK = 2
    buf = Vector{Ptr{Cchar}}(undef,1)
    sz = Csize_t(0)
    bytes = ccall(
        (:nng_recv, "libnng"),
        Int32,
        (Ref{nng_socket}, Ptr{Cchar}, Ref{Csize_t}, Cint),
        socket, pointer(buf), sz, 1)
    @show bytes
    msg = unsafe_string(buf[1])
    ccall((:nn_freemsg, "libnng"), Cint, (Ptr{Cchar},), buf[1])
    return msg
end

function node0(url)
    sock = nng_socket(0)
    rv = nng_pull0_open(sock)
    check_err("nng_pull0_open", rv)
    lp = nng_listener(0)
    rv = nng_listen(sock, url, lp)
    check_err("nng_listen", rv)
    while true
        msg = nng_recv(sock)
        check_err("nn_recv", msg)
        println("NODE0: RECEIVED ", msg)
    end
end

url = "ipc:///tmp/pipeline.ipc"
#node0(url)


text:
# Read in file of helper functions
# This file incluse: cleantext, countwords
include("../word_count_helpers.jl")

# Load the file names
dataLoc = "../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Create function for processing each file
function getwordcounts(fname)
    # Read in file and clean the text
    f = open(fname)
    text=cleantext(readlines(f))

    # Count number of times each word appears
    return countwords(text)
end

# Get counts for each book
allcounts = getwordcounts.(fnames)

# Calculate the overall word counts
overallcounts = merge(+,allcounts...)

# Sort and print the top 5 words with their counts
println(sort(collect(overallcounts), by=x->x[2], rev=true)[1:5,:])

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
using Flux, Zygote
using LinearAlgebra
using jInv.Mesh
using Printf
using Plots
using JLD
using MFGnet
include("viewers.jl")
include("runOMThelpers.jl")

# set default experimental parameters
@isdefined(R)  ? R = R : R = Float64                                  # element type. Use Float32 or Float64
@isdefined(d)  ? d = d : d = 2                                        # dimension of problem
@isdefined(ar) ? ar = ar : ar = x-> R.(x)                             # function to change element type of arrays  and device (e.g., to cuArray)
@isdefined(m)           ? m=m               : m=16                    # width of network
@isdefined(nTh)         ? nTh=nTh           : nTh=2                   # number of nodes in ResNet discretization
@isdefined(nTrain)      ? nTrain=nTrain     : nTrain=[32^2 64^2 128^2]    # number of training samples
@isdefined(nVal)        ? nVal=nVal         : nVal=minimum([32^2,nTrain[end]]) # number of validation samples
@isdefined(nt)          ? nt=nt             : nt=2                   # number of time steps for characteristics
@isdefined(stepper)     ? stepper=stepper   : stepper=RK4Step()       # time stepping for characteristics
@isdefined(T)           ? T=T               : T=1.0                   # final time for dynamical OT
@isdefined(alph)        ? α=alph             : α=[1.0;1.0;5.0;1.0;5.0] # weights for objective
@isdefined(sampleFreq)  ? sampleFreq = sampleFreq : sampleFreq = 25   # sample frequency
@isdefined(saveIter)    ? saveIter = saveIter     : saveIter   = 25   # iteration to save weights
@isdefined(maxIter)     ? maxIter  = maxIter      : maxIter = [200 200 200]     # max number of iters
@isdefined(optim)       ? optim  = optim      : optim = :bfgs     # max number of iters

@isdefined(saveStr)     ? saveStr=saveStr   : saveStr = "OMT-Multilevel-d-"*string(d)

sig = ar(.3*ones(R,d))
rho1 = Gaussian(d,sig,zeros(R,d),R(1.0))
Gs = Array{Gaussian{R,Vector{R}}}(undef,8)
ang = range(0,stop=2*pi, length=length(Gs)+1)[1:end-1]
for k=1:length(ang)
    μk = ar(4*([cos(ang[k]); sin(ang[k]);zeros(d-2)]))
    Gs[k] = Gaussian(d,sig,μk,R(1.0/length(Gs)))
end
rho0 = GaussianMixture(Gs)

################################################################################
# Note: validation from grid
################################################################################
# setup samples
domain = 5.0*[-1 1 -1 1];
# M      = getRegularMesh(domain,[sqrt(nTrain),sqrt(nTrain)])
# X0val     = Matrix(getCellCenteredGrid(M)')
# X0val   = ar([X0val; zeros(d-2,size(X0val,2))]); nVal = size(X0val,2)
# rho0xVal = rho0(X0val)
# rho1xVal = rho1(X0val)
#
# vol     = 10^d
# wVal    = (vol/nVal) * rho0xVal
# sum(wVal)

################################################################################
# validation from samples
################################################################################
# nVal    = minimum([64^2,nTrain])
X0val   = sample(rho0,nVal)
wVal    = ar(1/nVal * ones(nVal))
println("sum(wVal)=$(sum(wVal))")
if norm(sum(wVal) - R(1)) > 0.1
    error("Double check the weights in validation loss.")
end
################################################################################
################################################################################
Φ = getPotentialResNet(nTh,1.0,nTh,R)
(w0,ΘN,A0,b0,z0) = initializeWeights(d,m,nTh,ar)

α = R.(α)
tspan = R.([0.0, T])

Fv = F0()
Gv   = Gkl(rho0,rho1,rho0(X0val),rho1(X0val),R.(1.0))
Jv  = MeanFieldGame(Fv,Gv,X0val,rho0,wVal,Φ=Φ,stepper=stepper,nt=nt,α=α,tspan=tspan)

Θ = (w0,(ΘN),A0,b0,z0)

parms = MFGnet.myMap(x->param(x),Θ)
ps = Flux.params(parms)

global bestLoss = Inf
global Θtemp = initializeWeights(d,m,nTh,ar)
global Θbest = initializeWeights(d,m,nTh,ar)

cb = function(J,Jv,level,iter,His,parms,doPlots=true)

    global bestLoss,  Θtemp,  Θbest

    θc = parms
    Jvc = Jv(parms)

    if Jvc < bestLoss
        bestLoss = Jvc.data
        # println("loss:", bestLoss)
        Θtemp = copy(MFGnet.param2vec(parms))
        Θbest = MFGnet.vec2param!(Θtemp,Θbest)
        # println("Θopt:", Θopt[1][1:10])
    end

    Hisk = [J.cs[1].data J.cs[2].data J.cs[3].data J.cs[4].data J.cs[5].data Jv.cs[1].data Jv.cs[2].data Jv.cs[3].data Jv.cs[4].data Jv.cs[5].data]
    His[iter,:] = Hisk

    if iter==1
        str1 = @sprintf("iter\tobj\t\tcostL\t\tcostF\t\tcostG\t\tcostHJ\t\tcostHJfinal\t\t\tobjVal\t\tcostLVal\tcostFVal\tcostGVal\tcostHJVal\tcostHJfinalVal\n")
        str2 = @sprintf("%05d\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t\t\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e",
                iter,sum(His[iter,1:5]),His[iter,1],His[iter,2],His[iter,3],His[iter,4],His[iter,5],
                sum(His[iter,6:10]),His[iter,6],His[iter,7],His[iter,8],His[iter,9],His[iter,10])
        str = str1*str2
    else
        str = @sprintf("%05d\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t\t\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e",
                iter,sum(His[iter,1:5]),His[iter,1],His[iter,2],His[iter,3],His[iter,4],His[iter,5],
                sum(His[iter,6:10]),His[iter,6],His[iter,7],His[iter,8],His[iter,9],His[iter,10])
    end
    println(str)

    # if level>1 && iter==1
    #     println("saving...")
    #     currentSaveStr = saveStr * "-level-$level-iter-$iter.jld"
    #     settings = (m,α,nTrain,R,d,domain,nTh,m,nt,T)
    #     save(currentSaveStr,"Θbest",Θbest,"His",His,"settings",settings)
    # end

    # save every saveIter iterations
    if mod(iter,saveIter)==0
        println("saving...")
        currentSaveStr = saveStr * "-level-$level-iter-$iter.jld"
        settings = (m,α,nTrain,R,d,domain,nTh,m,nt,T)
        save(currentSaveStr,"Θbest",Θbest,"His",His,"settings",settings)
    end

    if mod(iter,sampleFreq)==0;
        println("resampling..")
        J.X0    = sample(rho0,nTrain[level])
        J.rho0x = rho0(J.X0)
        J.G.rho0x = rho0(J.X0)
        J.G.rho1x = rho1(J.X0)
    end

    if doPlots
        p1= viewImage2D(Jv.rho0x,M,aspect_ratio=:equal)
        title!("rho0(x)")
        p2= viewImage2D(Jv.G.rho1x,M,aspect_ratio=:equal)
        title!("rho1(x)")
        detDy = exp.(Jv.UN[d+1,:]).data
        rho1y = (Jv.G.rho1(Jv.UN[1:d,:]).*detDy).data
        p3= viewImage2D(rho1y ,M,aspect_ratio=:equal,clims=(minimum(Jv.rho0x),maximum(Jv.rho0x)))
        # p3= viewImage2D(rho1y ,M,aspect_ratio=:equal)
        title!("rho1(y).*det")

        res = Jv.rho0x - rho1y
        p4 = viewImage2D(log10.(abs.(res)),M,aspect_ratio=:equal)
        title!("rho0(x)-rho1(y).*det") # where's the determinant

        p5 = viewImage2D(detDy,M,aspect_ratio=:equal)
        title!("det")

        xc = [ Jv.X0; fill(0.0,1,size(Jv.X0,2))]
        Φ0 = Jv.Φ(xc,θc).data
        p6 = viewImage2D(Φ0,M,aspect_ratio=:equal)
        title!("Phi(x,0)")

        ################
        # p6 = viewImage2D(Q(Jv.X0), M, aspect_ratio=:equal)
        # title!("Potential Q(x)")

        xc = ar([ Jv.UN[1:d,:].data; fill(1.0,1,size(Jv.X0,2))])
        Φ1 = vec(Jv.Φ(xc,θc).data)
        p7 = viewImage2D(Φ1,M,aspect_ratio=:equal)
        title!("Phi(z,1)")

        deltaG = Jv.α[3]*MFGnet.getDeltaG(Jv.G,Jv.UN).data
        p8 = viewImage2D(deltaG,M,aspect_ratio=:equal)
        title!("deltaG(z)")
        p9 = viewImage2D(log10.(abs.(Φ1-deltaG)),M,aspect_ratio=:equal)
        title!("Phi1 - deltaG(z)")

        Xt    = sample(rho0,10)
        if d>2
            Xt[3:d,:] .= R(0.0)
        end
        wt     = 1/size(Xt,2) * ones(size(Xt,2))
        rho0xt = rho0(Xt)
        rho1xt = rho1(Xt)
        Ft = F0()
        Gt = Gkl(rho0,rho1,rho0xt,rho1xt,ar(1.0))
        Jt = MeanFieldGame(Ft,Gt,Xt,rho0,wt,Φ=Φ,stepper=stepper,nt=nt,α=α,tspan=tspan)
        ntv = Jt.nt
        Ut = MFGnet.integrate2(Jt.stepper,MFGnet.odefun,Jt,ar([Xt;zeros(4,size(Xt,2))]),θc,ar(tspan),ntv)
        # Y0  = ar([Ut[1:d,:,end];zeros(4,size(Xt,2))])
        # U0 = MFGnet.integrate2(Jt.stepper,MFGnet.odefun,Jt,Y0,Θopt,ar([T 0.0]),ntv)

        for j=1:size(Xt,2)
            uj = Ut[1:d,j,:]
            # vj = U0[1:d,j,:]
            plot!(p1,uj[1,:],uj[2,:],legend=false,linecolor=:white,linewidth=1.5)
            # plot!(p1,vj[1,:],vj[2,:],legend=false,linecolor=:red,linewidth=1.5)
        end

        pt = plot(p1,p2,p3,p4,p5,p6,p7,p8,p9)
        display(pt)
    end
end

for level=1:length(nTrain)
    nTr = nTrain[level]
    maxIt = maxIter[level]
    X0    = sample(rho0,nTr)
    w     = ar(1/nTr * ones(nTr))

    F = F0()
    G = Gkl(rho0,rho1,rho0(X0),rho1(X0),ar(1.0))
    J = MeanFieldGame(F,G,X0,rho0,w,Φ=Φ,stepper=stepper,nt=nt,α=α,tspan=tspan)

    println("\n\n ---------- Multilevel OMT Driver: level=$level of $(length(nTrain))-------------\n\n")
    println("results stored in: $(pwd())/$saveStr-level-$level.jld)")
    println("sampleFreq = $(sampleFreq), α = $(α), nTh = $(nTh), m = $(m), nt = $(J.nt)")
    println("DIMENSION = $(d), nTrain = $(nTr), nVal = $(nVal), saveIter = $(saveIter), maxIter = $(maxIt)\n\n")
    println("optim: $optim")

    His = zeros(maxIt,10)
    cbBFGS = (iter)-> cb(J,Jv,level,iter,His,parms,false)

    if optim==:bfgs
        f   =(Θ)-> evalObj(J,Θ,parms,ps)
        fdf =(Θ)-> evalObjAndGrad(J,Θ,parms,ps)
        Θ0 = MFGnet.param2vec(parms)
        f(Θ0)
        fdf(Θ0)
        runtime = @elapsed Θopt,flag,His,X,H = MFGnet.bfgs(f,fdf,Θ0,maxIter=size(His,1),
                                                   out=0,atol=1e-10,cb=cbBFGS)

    else
        println("using ADAM")
        opt = ADAM()
        runtime = @elapsed for k=1:maxIter
            if mod(k,sampleFreq)==0
                Xt    = sample(J.rho0,nTrain)
                w     = (1/nTrain) * ones(R,nTrain)
                J.X0  = Xt;
                J.w   = w;
                J.G.rho0x = J.G.rho0(Xt)
                J.G.rho1x = J.G.rho1(Xt)
                J.rho0x = J.rho0(Xt)
            end

            Jc,back = Zygote.Tracker.forward(() -> J(parms), ps)
            grads = back(1)
            ng = 0.0
            for p in ps
                ng += norm(grads[p]).^2
                Zygote.Tracker.update!(opt,p, grads[p])
            end
            cbBFGS(k)
        end
    end
    println("average time per iteration: $(runtime/maxIt)")

    Θopt = MFGnet.param2vec(parms)
    Θopt = MFGnet.vec2param!(Θopt,Θ)

    settings = (m,α,nTrain,R,d,domain,nTh,m,nt,T)

    save(saveStr*"-level-$level.jld","Θopt", Θopt,"Θbest",Θbest,"His",His,"settings",settings)
end


text:
#=
Author : Ng Yeesian
Modified : Guy Farmer
generates the Maximum Availability Location Problem model
=#
struct MALPDeployment <: DeploymentModel
    m::JuMP.Model
    x::Vector{JuMP.VariableRef}
end
deployment(m::MALPDeployment) = Int[round(Int,x) for x in JuMP.value.(m.x)]

function MALPDeployment(p::DeploymentProblem,
                        q::Float64; # busy fraction;
                        α::Float64 = 0.99, # reliability level
                        tol = params.δ,
                        solver = Gurobi.Optimizer(OutputFlag=0))
    demand = vec(mean(p.demand[p.train,:],dims = 1))
    @assert length(demand) == p.nregions

    I = 1:p.nlocations
    J = 1:p.nregions
    b = ceil(Int, log(1-α)/log(q))

    m = Model(GLPK.Optimizer)
    solver=Gurobi.Optimizer(OutputFlag=0)
    JuMP.@variable(m, x[1:p.nlocations] >= 0, Int)
    JuMP.@variable(m, z[1:p.nregions, 1:b], Bin)

    #possible breaking change with for j in J
    JuMP.@objective(m, Max, sum(demand[j]*z[j,b] for j in J))

    JuMP.@constraint(m, sum(x[i] for i in I) <= p.nambulances)
    for j in J # coverage over all regions
        JuMP.@constraint(m, sum(x[i] for i in filter(i->p.coverage[j,i], I)) >= 1)
        JuMP.@constraint(m, sum(x[i] for i in filter(i->p.coverage[j,i], I)) >= sum(z[j,k] for k in 1:b))
        for k in 2:b
            JuMP.@constraint(m, z[j,k] <= z[j,k-1])
        end
    end
    for i in I
        JuMP.@constraint(m, x[i] <= 5)
    end

    MALPDeployment(m, x)
end

optimize!(model::MALPDeployment) = JuMP.optimize!(model.m)
optimize!(model::MALPDeployment, p::DeploymentProblem) = JuMP.optimize!(model.m)

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
using DataStructures

CLEAN, WEAKENED, INFECTED, FLAGGED = 0, 1, 2, 3

d = readlines("input.txt")
grid = DefaultDict(0)
for y in reverse(1:length(d))
    for x in reverse(1:length(d))
        grid[[x, y]] = d[y][x] == '.' ? CLEAN : INFECTED
    end
end
dir = Dict("u" => [0, -1], "r" => [1, 0], "d" => [0, 1], "l" => [-1, 0])
dir_map = Dict("u" => ["l", "r"], "r" => ["u", "d"], "d" => ["r", "l"], "l" => ["d", "u"])
next_dir(cur_dir, v) = v == INFECTED ? dir_map[cur_dir][2] : dir_map[cur_dir][1]

function iter(grid, p, cur_dir, part_one, changes)
    if grid[p] == CLEAN
        cur_dir = next_dir(cur_dir, CLEAN)
        if part_one
            changes += 1
        end
    elseif grid[p] == WEAKENED && !part_one
        changes += 1
    elseif grid[p] == INFECTED
        cur_dir = next_dir(cur_dir, INFECTED)
    else
        cur_dir = next_dir(next_dir(cur_dir, INFECTED), INFECTED)
    end
    grid[p] = (grid[p] + 1 + part_one) % 4
    p += dir[cur_dir]
    return p, cur_dir, changes
end

function solve(grid, part_one, n)
    p = [Int(ceil(length(d)/2)), Int(ceil(length(d)/2))]
    cur_dir = "u"
    changes = 0
    for _ in 1:n
        p, cur_dir, changes = iter(grid, p, cur_dir, part_one, changes)
    end
    println(changes)
end
solve(copy(grid), true, 10_000)
solve(grid, false, 10_000_000)


text:
using Test
using ConvOptDL: _get_block_mask, DropBlock, dropblock
using Flux
using Zygote: pullback

@testset "DropBlock Test" begin
    mask = reshape(
        [
            [0, 0, 0, 0, 0]
            [0, 1, 0, 0, 0]
            [0, 0, 0, 0, 0]
            [0, 1, 0, 0, 0]
            [0, 0, 0, 0, 0]
        ],
        5,
        5,
        1,
        1,
    )
    expected = reshape(
        [
            [0, 0, 0, 1, 1]
            [0, 0, 0, 1, 1]
            [0, 0, 0, 1, 1]
            [0, 0, 0, 1, 1]
            [0, 0, 0, 1, 1]
        ],
        5,
        5,
        1,
        1,
    )
    block = DropBlock(3, 0.1)
    block_mask = _get_block_mask(block, rand(5, 5, 1, 1), mask)
    @test block_mask == expected

    testmode!(block)
    x = rand(5, 5, 3, 10)
    @test block(x) == x

    (_, block_mask), back = pullback(dropblock, x, block)
    Δy = ones(size(x))
    Δx, _  = back((Δy, nothing))
    @test Δx == Δy .* block_mask

    block = DropBlock(3, 0)
    testmode!(block)
    Δy, = gradient((x) -> sum(block(x)), x)
    Δŷ, = gradient(sum, x)
    @test Δy == Δŷ

    trainmode!(block)
    Δy, = gradient((x) -> sum(block(x)), x)
    Δŷ, = gradient(sum, x)
    @test Δy != Δŷ
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
    flexible_demand(EP::Model, inputs::Dict, CapacityReserveMargin::Int)
This function defines the operating constraints for flexible demand resources. As implemented, flexible demand resources ($y \in \mathcal{DF}$) are characterized by: a) maximum deferrable demand as a fraction of available capacity in a particular time step $t$, $\rho^{max}_{y,z,t}$, b) the maximum time this demand can be advanced and delayed, defined by parameters, $\tau^{advance}_{y,z}$ and $\tau^{delay}_{y,z}$, respectively and c) the energy losses associated with shifting demand, $\eta_{y,z}^{dflex}$.
**Tracking total deferred demand**
The operational constraints governing flexible demand resources are as follows.
The first two constraints model keep track of inventory of deferred demand in each time step.  Specifically, the amount of deferred demand remaining to be served ($\Gamma_{y,z,t}$) depends on the amount in the previous time step minus the served demand during time step $t$ ($\Theta_{y,z,t}$) while accounting for energy losses associated with demand flexibility, plus the demand that has been deferred during the current time step ($\Pi_{y,z,t}$). Note that variable $\Gamma_{y,z,t} \in \mathbb{R}$, $\forall y \in \mathcal{DF}, t  \in \mathcal{T}$. Similar to hydro inventory or storage state of charge constraints, for the first time step of the year (or each representative period), we define the deferred demand level based on level of deferred demand in the last time step of the year (or each representative period).
```math
\begin{aligned}
\Gamma_{y,z,t} = \Gamma_{y,z,t-1} -\eta_{y,z}^{dflex}\Theta_{y,z,t} +\Pi_{y,z,t} \hspace{4 cm}  \forall y \in \mathcal{DF}, z \in \mathcal{Z}, t \in \mathcal{T}^{interior} \\
\Gamma_{y,z,t} = \Gamma_{y,z,t +\tau^{period}-1} -\eta_{y,z}^{dflex}\Theta_{y,z,t} +\Pi_{y,z,t} \hspace{4 cm}  \forall y \in \mathcal{DF}, z \in \mathcal{Z}, t \in \mathcal{T}^{start}
\end{aligned}
```
**Bounds on available demand flexibility**
At any given time step, the amount of demand that can be shifted or deferred cannot exceed the maximum deferrable demand, defined by product of the availability factor ($\rho^{max}_{y,t}$) times the available capacity($\Delta^{total}_{y,z}$).
```math
\begin{aligned}
\Pi_{y,t} \leq \rho^{max}_{y,z,t}\Delta_{y,z} \hspace{4 cm}  \forall y \in \mathcal{DF}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
**Maximum time delay and advancements**
Delayed demand must then be served within a fixed number of time steps. This is done by enforcing the sum of demand satisfied ($\Theta_{y,z,t}$) in the following $\tau^{delay}_{y,z}$ time steps (e.g., t + 1 to t + $\tau^{delay}_{y,z}$) to be greater than or equal to the level of energy deferred during time step $t$.
```math
\begin{aligned}
\sum_{e=t+1}^{t+\tau^{delay}_{y,z}}{\Theta_{y,z,e}} \geq \Gamma_{y,z,t}
    \hspace{4 cm}  \forall y \in \mathcal{DF},z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
A similar constraints maximum time steps of demand advancement. This is done by enforcing the sum of demand deferred ($\Pi_{y,t}$) in the following $\tau^{advance}_{y}$ time steps (e.g., t + 1 to t + $\tau^{advance}_{y}$) to be greater than or equal to the total level of energy deferred during time $t$ (-$\Gamma_{y,t}$). The negative sign is included to account for the established sign convention that treat demand deferred in advance of the actual demand is defined to be negative.
```math
\begin{aligned}
\sum_{e=t+1}^{t+\tau^{advance}_{y,z}}{\Pi_{y,z,e}} \geq -\Gamma_{y,z,t}
    \hspace{4 cm}  \forall y \in \mathcal{DF}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
If $t$ is first time step of the year (or the first time step of the representative period), then the above two constraints are implemented to look back over the last n time steps, starting with the last time step of the year (or the last time step of the representative period). This time-wrapping implementation is similar to the time-wrapping implementations used for defining the storage balance constraints for hydropower reservoir resources and energy storage resources.
"""
function flexible_demand(EP::Model, inputs::Dict, CapacityReserveMargin::Int)
    ## Flexible demand resources available during all hours and can be either delayed or advanced (virtual storage-shiftable demand) - DR ==1

    println("Flexible Demand Resources Module")

    dfGen = inputs["dfGen"]

    T = inputs["T"]     # Number of time steps (hours)
    Z = inputs["Z"]     # Number of zones
    FLEX = inputs["FLEX"] # Set of flexible demand resources

    START_SUBPERIODS = inputs["START_SUBPERIODS"]
    INTERIOR_SUBPERIODS = inputs["INTERIOR_SUBPERIODS"]

    hours_per_subperiod = inputs["hours_per_subperiod"] # Total number of hours per subperiod

    END_HOURS = START_SUBPERIODS .+ hours_per_subperiod .- 1 # Last subperiod of each representative period

    ### Variables ###

    # Variable tracking total advanced (negative) or deferred (positive) demand for demand flex resource y in period t
    @variable(EP, vS_FLEX[y in FLEX, t = 1:T])

    # Variable tracking demand deferred by demand flex resource y in period t
    @variable(EP, vCHARGE_FLEX[y in FLEX, t = 1:T] >= 0)

    ### Expressions ###

    ## Power Balance Expressions ##
    @expression(
        EP,
        ePowerBalanceDemandFlex[t = 1:T, z = 1:Z],
        sum(
            -EP[:vP][y, t] + EP[:vCHARGE_FLEX][y, t] for
            y in intersect(FLEX, dfGen[(dfGen[!, :Zone].==z), :][!, :R_ID])
        )
    )

    EP[:ePowerBalance] += ePowerBalanceDemandFlex

    # Capacity Reserves Margin policy
    if CapacityReserveMargin > 0
        @expression(
            EP,
            eCapResMarBalanceFlex[res = 1:inputs["NCapacityReserveMargin"], t = 1:T],
            sum(
                dfGen[y, Symbol("CapRes_$res")] * (EP[:vCHARGE_FLEX][y, t] - EP[:vP][y, t])
                for y in FLEX
            )
        )
        EP[:eCapResMarBalance] += eCapResMarBalanceFlex
    end

    ## Objective Function Expressions ##

    # Variable costs of "charging" for technologies "y" during hour "t" in zone "z"
    @expression(
        EP,
        eCVarFlex_in[y in FLEX, t = 1:T],
        inputs["omega"][t] * dfGen[!, :Var_OM_Cost_per_MWh_In][y] * vCHARGE_FLEX[y, t]
    )

    # Sum individual resource contributions to variable charging costs to get total variable charging costs
    @expression(EP, eTotalCVarFlexInT[t = 1:T], sum(eCVarFlex_in[y, t] for y in FLEX))
    @expression(EP, eTotalCVarFlexIn, sum(eTotalCVarFlexInT[t] for t = 1:T))
    EP[:eObj] += eTotalCVarFlexIn

    ### Constraints ###

    ## Flexible demand is available only during specified hours with time delay or time advance (virtual storage-shiftable demand)
    for z = 1:Z
        # NOTE: Flexible demand operates by zone since capacity is now related to zone demand
        FLEX_Z = intersect(FLEX, dfGen[dfGen[!, :Zone].==z, :][!, :R_ID])

        @constraints(
            EP,
            begin
                # State of "charge" constraint (equals previous state + charge - discharge)
                # NOTE: no maximum energy "stored" or deferred for later hours
                # NOTE: Flexible_Demand_Energy_Eff corresponds to energy loss due to time shifting
                [y in FLEX_Z, t in INTERIOR_SUBPERIODS],
                EP[:vS_FLEX][y, t] ==
                EP[:vS_FLEX][y, t-1] -
                dfGen[!, :Flexible_Demand_Energy_Eff][y] * (EP[:vP][y, t]) +
                (EP[:vCHARGE_FLEX][y, t])
                # Links last time step with first time step, ensuring position in hour 1 is within eligible change from final hour position
                [y in FLEX_Z, t in START_SUBPERIODS],
                EP[:vS_FLEX][y, t] ==
                EP[:vS_FLEX][y, t+hours_per_subperiod-1] -
                dfGen[!, :Flexible_Demand_Energy_Eff][y] * (EP[:vP][y, t]) +
                (EP[:vCHARGE_FLEX][y, t])

                # Maximum charging rate
                # NOTE: the maximum amount that can be shifted is given by hourly availability of the resource times the maximum capacity of the resource
                [y in FLEX_Z, t = 1:T],
                EP[:vCHARGE_FLEX][y, t] <= inputs["pP_Max"][y, t] * EP[:eTotalCap][y]
                # NOTE: no maximum discharge rate unless constrained by other factors like transmission, etc.
            end
        )

        for y in FLEX_Z

            # Require deferred demands to be satisfied within the specified time delay
            max_flexible_demand_delay = Int(floor(dfGen[!, :Max_Flexible_Demand_Delay][y]))
            FLEXIBLE_DEMAND_DELAY_HOURS = [] # Set of hours in the summation term of the maximum demand delay constraint for the first subperiod of each representative period
            for s in START_SUBPERIODS
                flexible_demand_delay_start =
                    s + hours_per_subperiod - max_flexible_demand_delay
                FLEXIBLE_DEMAND_DELAY_HOURS = union(
                    FLEXIBLE_DEMAND_DELAY_HOURS,
                    flexible_demand_delay_start:(s+hours_per_subperiod-2),
                )
            end

            @constraints(
                EP,
                begin
                    # cFlexibleDemandDelay: Constraints looks back over last n hours, where n = dfGen[!,:Max_Flexible_Demand_Delay][y]
                    [t in setdiff(1:T, FLEXIBLE_DEMAND_DELAY_HOURS, END_HOURS)],
                    sum(
                        EP[:vP][y, e] for
                        e = (t+1):(t+Int(floor(dfGen[!, :Max_Flexible_Demand_Delay][y])))
                    ) >= EP[:vS_FLEX][y, t]

                    # cFlexibleDemandDelayWrap: If n is greater than the number of subperiods left in the period, constraint wraps around to first hour of time series
                    # cFlexibleDemandDelayWrap constraint is equivalant to: sum(EP[:vP][y,e] for e=(t+1):(hours_per_subperiod_max)+sum(EP[:vP][y,e] for e=hours_per_subperiod_min:(hours_per_subperiod_min-1+dfGen[!,:Max_Flexible_Demand_Delay][y]-(hours_per_subperiod-(t%hours_per_subperiod)))) >= EP[:vS_FLEX][y,t]
                    [t in FLEXIBLE_DEMAND_DELAY_HOURS],
                    sum(
                        EP[:vP][y, e] for
                        e = (t+1):(t+hours_per_subperiod-(t%hours_per_subperiod))
                    ) + sum(
                        EP[:vP][y, e] for e =
                            (hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1):((hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1)-1+Int(
                                floor(dfGen[!, :Max_Flexible_Demand_Delay][y]),
                            )-(hours_per_subperiod-(t%hours_per_subperiod)))
                    ) >= EP[:vS_FLEX][y, t]

                    # cFlexibleDemandDelayEnd: cFlexibleDemandDelayEnd constraint is equivalant to: sum(EP[:vP][y,e] for e=hours_per_subperiod_min:(hours_per_subperiod_min-1+dfGen[!,:Max_Flexible_Demand_Delay][y])) >= EP[:vS_FLEX][y,t]
                    [t in END_HOURS],
                    sum(
                        EP[:vP][y, e] for e =
                            (hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1):((hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1)-1+Int(floor(dfGen[!, :Max_Flexible_Demand_Delay][y])))
                    ) >= EP[:vS_FLEX][y, t]

                    # NOTE: Expression (hours_per_subperiod*Int(floor((t-1)/hours_per_subperiod))+1) is equivalant to "hours_per_subperiod_min"
                    # NOTE: Expression t+hours_per_subperiod-(t%hours_per_subperiod) is equivalant to "hours_per_subperiod_max"
                end
            )

            # Require advanced demands to be satisfied within the specified time period
            max_flexible_demand_advance =
                Int(floor(dfGen[!, :Max_Flexible_Demand_Advance][y]))
            FLEXIBLE_DEMAND_ADVANCE_HOURS = [] # Set of hours in the summation term of the maximum advance demand constraint for the first subperiod of each representative period
            for s in START_SUBPERIODS
                flexible_demand_advance_start =
                    s + hours_per_subperiod - max_flexible_demand_advance
                FLEXIBLE_DEMAND_ADVANCE_HOURS = union(
                    FLEXIBLE_DEMAND_ADVANCE_HOURS,
                    flexible_demand_advance_start:(s+hours_per_subperiod-2),
                )
            end

            @constraints(
                EP,
                begin
                    # cFlexibleDemandAdvance: Constraint looks back over last n hours, where n = dfGen[!,:Max_Flexible_Demand_Advance][y]
                    [t in setdiff(1:T, FLEXIBLE_DEMAND_ADVANCE_HOURS, END_HOURS)],
                    sum(
                        EP[:vCHARGE_FLEX][y, e] for
                        e = (t+1):(t+Int(dfGen[!, :Max_Flexible_Demand_Advance][y]))
                    ) >= -EP[:vS_FLEX][y, t]

                    # cFlexibleDemandAdvanceWrap: If n is greater than the number of subperiods left in the period, constraint wraps around to first hour of time series
                    # cFlexibleDemandAdvanceWrap constraint is equivalant to: sum(EP[:vCHARGE_FLEX][y,e] for e=(t+1):hours_per_subperiod_max)+sum(EP[:vCHARGE_FLEX][y,e] for e=hours_per_subperiod_min:(hours_per_subperiod_min-1+dfGen[!,:Max_Flexible_Demand_Advance][y]-(hours_per_subperiod-(t%hours_per_subperiod)))) >= -EP[:vS_FLEX][y,t]
                    [t in FLEXIBLE_DEMAND_ADVANCE_HOURS],
                    sum(
                        EP[:vCHARGE_FLEX][y, e] for
                        e = (t+1):(t+hours_per_subperiod-(t%hours_per_subperiod))
                    ) + sum(
                        EP[:vCHARGE_FLEX][y, e] for e =
                            (hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1):((hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1)-1+Int(
                                dfGen[!, :Max_Flexible_Demand_Advance][y],
                            )-(hours_per_subperiod-(t%hours_per_subperiod)))
                    ) >= -EP[:vS_FLEX][y, t]

                    # cFlexibleDemandAdvanceEnd: cFlexibleDemandAdvanceEnd constraint is equivalant to: sum(EP[:vCHARGE_FLEX][y,e] for e=hours_per_subperiod_min:(hours_per_subperiod_min-1+dfGen[!,:Max_Flexible_Demand_Advance][y])) >= -EP[:vS_FLEX][y,t]
                    [t in END_HOURS],
                    sum(
                        EP[:vCHARGE_FLEX][y, e] for e =
                            (hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1):((hours_per_subperiod*Int(
                                floor((t - 1) / hours_per_subperiod),
                            )+1)-1+Int(dfGen[!, :Max_Flexible_Demand_Advance][y]))
                    ) >= -EP[:vS_FLEX][y, t]
                    # NOTE: Expression (hours_per_subperiod*Int(floor((t-1)/hours_per_subperiod))+1) is equivalant to "hours_per_subperiod_min"
                    # NOTE: Expression t+hours_per_subperiod-(t%hours_per_subperiod) is equivalant to "hours_per_subperiod_max"
                end
            )
        end
    end

    return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using JuLIPMaterials.CLE: onb3D, euclidean, spherical
using JuLIPMaterials: Vec3, Mat3, Ten33, Ten43, ForceConstantMatrix1

using Einsum, StaticArrays, ForwardDiff, LinearAlgebra

export GreenFunctionCorrector

abstract type AbstractGreenFunctionCorrector{T} end

struct GreenFunctionCorrector3D{T} <: AbstractGreenFunctionCorrector{T}
   Nquad::Int
   C
   FCM::ForceConstantMatrix1{T}
   remove_singularity::Bool
end

function GreenFunctionCorrector(ℂ, FCM::ForceConstantMatrix1; Nquad = nothing, remove_singularity = true)
   if Nquad == nothing
      error("still need to implement auto-tuning, please provide Nquad")
   end
   return GreenFunctionCorrector3D(Nquad, ℂ, FCM, remove_singularity)
end

function (Gcorr::GreenFunctionCorrector3D{T})(x) where T
   if Gcorr.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 3, 3)
   end
   return eval_corrector(Vec3(x), Gcorr.C, Gcorr.FCM, Gcorr.Nquad)
end

# Utility functions for computing Fourier multipliers
function _C4(k::Vec3{TT},FCM::ForceConstantMatrix1{SS}) where {TT, SS}
   H0 = zeros(TT,3,3)
   for i=1:length(FCM.R)
      ρ = FCM.R[i]
      Hmat = FCM.H[i]
      H0 += (dot(k, ρ)^4)/24*Hmat
   end
   return H0
end

function _dC4(k::Vec3{TT},FCM::ForceConstantMatrix1{SS}) where {TT, SS}
   H0 = zeros(TT,3,3,3)
   for i=1:length(FCM.R)
      ρ = FCM.R[i]
      Hmat = FCM.H[i]
      @einsum H0[a,b,c] += (dot(k, ρ)^3)/6 * Hmat[a,b] * ρ[c]
   end
   return H0
end

function _ddC4(k::Vec3{TT},FCM::ForceConstantMatrix1{SS}) where {TT, SS}
   H0 = zeros(TT,3,3,3,3)
   for i=1:length(FCM.R)
      ρ = FCM.R[i]
      Hmat = FCM.H[i]
      @einsum H0[a,b,c,d] += dot(k, ρ)^2/2 * Hmat[a,b] * ρ[c] * ρ[d]
   end
   return H0
end

function _C2(k::Vec3{TT},ℂ::Array{SS,4}) where {TT, SS}
   C2 = zeros(TT,3,3)
   @einsum C2[i,j] = k[a] * ℂ[i,a,j,b] * k[b]
   return C2;
end

function _dC2(k::Vec3{TT},ℂ::Array{SS,4}) where {TT, SS}
   dC2 = zeros(TT,3,3,3)
   @einsum dC2[i,j,b] = ℂ[i,a,j,b] * k[a] + ℂ[i,b,j,a] * k[a]
   return dC2;
end

function _C2inv(k::Vec3{TT},ℂ::Array{SS,4}) where {TT, SS}
   C2 = _C2(k,ℂ);
   return inv(C2);
end

function _dC2inv(k::Vec3{TT},ℂ::Array{SS,4}) where {TT, SS}
   dC2inv = zeros(TT,3,3,3)
   dC2 = _dC2(k,ℂ)
   C2inv = _C2inv(k,ℂ)
   @einsum dC2inv[a,b,c] = -C2inv[a,d]*dC2[d,e,c]*C2inv[e,b]
   return dC2inv;
end

function _dC2inv1(k::Vec3{TT},ℂ::Array{SS,4}) where {TT, SS}
   dC2inv = zeros(TT,3,3,3)
   C2 = _C2(k,ℂ)
   dC2 = _dC2(k,ℂ)
   for c=1:3
      dC2inv[:,:,c] = -C2 \ dC2[:,:,c] / C2
   end
   return dC2inv;
end

function _ddC2inv(k::Vec3{TT},ℂ::Array{SS,4}) where {TT, SS}
   ddC2inv = zeros(TT,3,3,3,3)
   # Get matrices
   dC2 = _dC2(k,ℂ)
   dC2inv = _dC2inv1(k,ℂ)
   C2inv = _C2inv(k,ℂ)
   @einsum ddC2inv[a,b,c,n] = -dC2inv[a,d,n]*dC2[d,e,c]*C2inv[e,b]-
                                 2*C2inv[a,d]*ℂ[d,e,c,n]*C2inv[e,b]-
                                    C2inv[a,d]*dC2[d,e,c]*dC2inv[e,b,n]
   return ddC2inv;
end

function _ddC2inv1(k::Vec3{TT},ℂ::Array{SS,4}) where {TT, SS}
   ddC2inv = zeros(TT,3,3,3,3)
   # Get matrices
   dC2 = _dC2(k,ℂ)
   dC2inv = _dC2inv1(k,ℂ)
   C2 = _C2(k,ℂ)
   for c=1:3, n=1:3
      ddC2inv[:,:,c,n] = -(dC2inv[:,:,n] * dC2[:,:,c]) / C2 -
                           2*(C2 \ ℂ[:,:,c,n] / C2) -
                           C2 \ (dC2[:,:,c] * dC2inv[:,:,n]);
   end
   return ddC2inv;
end

# Multiplier for lattice elasticity
function _full_lattice_multiplier(k::Vec3{TT},FCM::ForceConstantMatrix1{SS}) where {TT, SS}
   D = @SMatrix zeros(TT,3,3)
   for i=1:length(FCM.R)
      D += -2*sin( 0.5*dot(k,FCM.R[i]) )^2 * FCM.H[i]
   end
   return D
end

# 0--homogeneous multiplier
function _corrector_multiplier(k::Vec3{TT},ℂ::Array{SS,4},FCM::ForceConstantMatrix1{SS}) where {TT, SS}
   C2inv = _C2inv(k,ℂ)
   C4 = _C4(k,FCM)
   return C2inv * C4 * C2inv
end

function _d_corrector_multiplier(k::Vec3{TT},ℂ::Array{SS,4},FCM::ForceConstantMatrix1{SS}) where {TT, SS}
   dH2 = zeros(TT,3,3,3)
   C2inv = _C2inv(k,ℂ)
   dC2inv = _dC2inv(k,ℂ)
   C4 = _C4(k,FCM)
   dC4 = _dC4(k,FCM)
   @einsum dH2[a,b,c] = dC2inv[a,i,c]*C4[i,j]*C2inv[j,b] + C2inv[a,i]*dC4[i,j,c]*C2inv[j,b] + C2inv[a,i]*C4[i,j]*dC2inv[j,b,c]
   return dH2
end

function _dd_corrector_multiplier(k::Vec3{TT},ℂ::Array{SS,4},FCM::ForceConstantMatrix1{SS}) where {TT, SS}
   ddH2 = zeros(TT,3,3,3,3)
   C2inv = _C2inv(k,ℂ)
   dC2inv = _dC2inv1(k,ℂ)
   ddC2inv = _ddC2inv1(k,ℂ)
   C4 = _C4(k,FCM)
   dC4 = _dC4(k,FCM)
   ddC4 = _ddC4(k,FCM)
   for a=1:3, b=1:3, c=1:3, d=1:3
       ddH2[a,b,c,d] = ddC2inv[a,:,c,d]' * C4 * C2inv[:,b] +
                       dC2inv[a,:,d]' * dC4[:,:,c] * C2inv[:,b] +
                       dC2inv[a,:,d]' * C4 * dC2inv[:,b,c] +
                       dC2inv[a,:,c]' * dC4[:,:,d] * C2inv[:,b] +
                       C2inv[a,:]' * ddC4[:,:,c,d] * C2inv[:,b] +
                       C2inv[a,:]' * dC4[:,:,d] * dC2inv[:,b,c] +
                       dC2inv[a,:,c]' * C4 * dC2inv[:,b,d] +
                       C2inv[a,:]' * dC4[:,:,c] * dC2inv[:,b,d] +
                       C2inv[a,:]' * C4 * ddC2inv[:,b,c,d];
   end
   return ddH2
end

"eval_corrector(x::Vec3, ℂ::Ten43, Nquad::Int)"
function eval_corrector(x::Vec3{TT}, ℂ::Array{SS,4},FCM::ForceConstantMatrix1{SS},Nquad::Int) where {TT, SS}
   # allocate
   Gcorr = @SMatrix zeros(TT,3,3)
   zz = @MMatrix zeros(TT,3,3)
   # Initialise tensors.
   x̂ = x/norm(x)
   # two vectors orthogonal to x.
   x1, x2 = onb3D(x̂)
   # Integrate
   DH2 = zeros(TT,3,3,3)
   ∂H2 = zeros(TT,3,3)
   D²H2 = zeros(TT,3,3,3,3)
   ∂²H2 = zeros(TT,3,3)
   for ω in range(0.0, step=pi/Nquad, length=Nquad)
      z = cos(ω) * x1 + sin(ω) * x2
      D²H2 = _dd_corrector_multiplier(z,ℂ,FCM)
      for i=1:3, j=1:3
         ∂²H2[i,j] = x̂' * D²H2[i,j,:,:] * x̂
      end
      DH2 = _d_corrector_multiplier(z,ℂ,FCM)
      @einsum ∂H2[i,j] = DH2[i,j,a]*z[a]
      Gcorr += -∂H2 + ∂²H2
   end
   # Normalise appropriately
   return Gcorr / (4*pi*norm(x)^3*Nquad)
end


text:
using ManifoldsBase
using Test

using ManifoldsBase: TraitList, merge_traits

struct PassthoughTrait <: AbstractTrait end

struct PassthroughDecorator{𝔽,MT<:AbstractManifold{𝔽}} <: AbstractDecoratorManifold{𝔽}
    manifold::MT
end

function ManifoldsBase.active_traits(f, ::PassthroughDecorator, ::Any...)
    return merge_traits(PassthoughTrait())
end

function ManifoldsBase.log!(
    ::TraitList{PassthoughTrait},
    M::AbstractDecoratorManifold,
    X,
    p,
    q,
)
    return log!(M.manifold, X, p, q)
end
function ManifoldsBase.exp!(
    ::TraitList{PassthoughTrait},
    M::AbstractDecoratorManifold,
    q,
    p,
    X,
)
    return log!(M.manifold, q, p, X)
end

@testset "PassthroughDecorator" begin
    M = PassthroughDecorator(ManifoldsBase.DefaultManifold(2))

    q = [0.0, 0.0]
    p = [0.0, 0.0]
    X = [1.0, 2.0]
    Y = [0.0, 0.0]
    @test inverse_retract!(M, q, p, X) == [1.0, 2.0]
    @test retract(M, p, q) == [1.0, 2.0]
    @test retract!(M, Y, p, q) == [1.0, 2.0]
end


text:
using ConstraintSolver, JuMP, MathOptInterface

if !@isdefined CS
    const CS = ConstraintSolver
end
const MOI = MathOptInterface
const MOIU = MOI.Utilities
include("../../test/sudoku_fcts.jl")

function from_file(filename, sep = '\n')
    s = open(filename) do file
        read(file, String)
    end
    str_sudokus = split(strip(s), sep)
    grids = AbstractArray[]
    for str_sudoku in str_sudokus
        str_sudoku = replace(str_sudoku, "." => "0")
        one_line_grid = parse.(Int, split(str_sudoku, ""))
        grid = reshape(one_line_grid, 9, 9)
        push!(grids, grid)
    end
    return grids
end

function solve_all(grids; benchmark = false, single_times = true)
    ct = time()
    grids = grids
    for (i, grid) in enumerate(grids)
        m = CS.Optimizer(logging = [])

        x = [[MOI.add_constrained_variable(m, MOI.Integer()) for i in 1:9] for j in 1:9]
        for r in 1:9, c in 1:9
            MOI.add_constraint(m, x[r][c][1], MOI.GreaterThan(1.0))
            MOI.add_constraint(m, x[r][c][1], MOI.LessThan(9.0))
        end

        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                sat = [MOI.ScalarAffineTerm(1.0, x[r][c][1])]
                MOI.add_constraint(
                    m,
                    MOI.ScalarAffineFunction{Float64}(sat, 0.0),
                    MOI.EqualTo(convert(Float64, grid[r, c])),
                )
            end
        end
        # sudoku constraints
        moi_add_sudoku_constr!(m, x)

        if single_times
            GC.enable(false)
            MOI.optimize!(m)
            status = MOI.get(m, MOI.TerminationStatus())
            GC.enable(true)
            # println(i - 1, ", ", MOI.get(m, MOI.SolveTime()))
        else
            GC.enable(false)
            MOI.optimize!(m)
            status = MOI.get(m, MOI.TerminationStatus())
            GC.enable(true)
        end
        if !benchmark
            println("Status: ", status)
            solution = zeros(Int, 9, 9)
            for r in 1:9
                solution[r, :] = [MOI.get(m, MOI.VariablePrimal(), x[r][c][1]) for c in 1:9]
            end
            @assert jump_fulfills_sudoku_constr(solution)
        end
    end
    # println("")
    tt = time() - ct
    # println("total time: ", tt)
    # println("avg: ", tt / length(grids))
end

function main(; benchmark = false, single_times = true)
    solve_all(
        from_file("data/top95.txt");
        benchmark = benchmark,
        single_times = single_times,
    )
    # solve_all(from_file("hardest.txt"), "hardest")
end


text:
# Fuzzy C means algorithm

## Interface

struct FuzzyCMeansResult{T<:AbstractFloat} <: ClusteringResult
    centers::Matrix{T}          # cluster centers (d x C)
    weights::Matrix{Float64}    # assigned weights (n x C)
    iterations::Int             # number of elasped iterations
    converged::Bool             # whether the procedure converged
end

## Utility functions

function update_weights!(weights, data, centers, fuzziness, dist_metric)
    pow = 2.0/(fuzziness-1)
    nrows, ncols = size(weights)
    dists = pairwise(dist_metric, data, centers, dims=2)
    for i in 1:nrows
        for j in 1:ncols
            den = 0.0
            for k in 1:ncols
                den += (dists[i,j]/dists[i,k])^pow
            end
            weights[i,j] = 1.0/den
        end
    end
end

function update_centers!(centers, data, weights, fuzziness)
    nrows, ncols = size(weights)
    T = eltype(centers)
    for j in 1:ncols
        num = zeros(T, size(data,1))
        den = zero(T)
        for i in 1:nrows
            δm = weights[i,j]^fuzziness
            num += δm * data[:,i]
            den += δm
        end
        centers[:,j] = num/den
    end
end

const _fcmeans_default_maxiter = 100
const _fcmeans_default_tol = 1.0e-3
const _fcmeans_default_display = :none

function fuzzy_cmeans(
    data::AbstractMatrix{T},
    C::Int,
    fuzziness::Real;
    maxiter::Int = _fcmeans_default_maxiter,
    tol::Real = _fcmeans_default_tol,
    dist_metric::Metric = Euclidean(),
    display::Symbol = _fcmeans_default_display
    ) where T<:Real

    nrows, ncols = size(data)
    2 <= C < ncols || error("C must have 2 <= C < n")
    1 < fuzziness || error("fuzziness must be greater than 1")

    _fuzzy_cmeans(data, C, fuzziness, maxiter, tol, dist_metric, display_level(display))

end

## Core implementation

function _fuzzy_cmeans(
    data::AbstractMatrix{T},                        # data matrix
    C::Int,                                         # total number of classes
    fuzziness::Real,                                # fuzziness
    maxiter::Int,                                   # maximum number of iterations
    tol::Real,                                      # tolerance
    dist_metric::Metric,                            # metric to calculate distance
    displevel::Int                                  # the level of display
    ) where T<:Real

    nrows, ncols = size(data)

    # Initialize weights randomly
    weights = rand(Float64, ncols, C)
    weights ./= sum(weights, dims=2)

    centers = zeros(T, (nrows, C))
    prev_centers = identity.(centers)

    δ = Inf
    iter = 0

    if displevel >= 2
        @printf "%7s %18s\n" "Iters" "center-change"
        println("----------------------------")
    end

    while iter < maxiter && δ > tol
        update_centers!(centers, data, weights, fuzziness)
        update_weights!(weights, data, centers, fuzziness, dist_metric)
        δ = maximum(colwise(dist_metric, prev_centers, centers))
        copyto!(prev_centers, centers)
        iter += 1
        if displevel >= 2
            @printf("%7d %18.6e\n", iter, δ)
        end
    end

    if displevel >= 1
        if δ <= tol
            println("Fuzzy C-means converged with $iter iterations (δ = $δ)")
        else
            println("Fuzzy C-means terminated without convergence after $iter iterations (δ = $δ)")
        end
    end

    FuzzyCMeansResult(centers, weights, iter, δ <= tol)
end

function Base.show(io::IO, result::FuzzyCMeansResult)
    d, C = size(result.centers)
    n, iter = size(result.weights, 1), result.iterations
    print(io, "FuzzyCMeansResult: $C clusters for $n points in $d dimensions ")
    if result.converged
        print(io, "(converged in $iter iterations)")
    else
        print(io, "(failed to converge in $iter iterations)")
    end
end


text:
using Zeros
using Test
using Test: @inferred

@testset "ambiguities" begin
    ambiguities = detect_ambiguities(Base, Zeros)
    for a in ambiguities
        println(a[1], "\n", a[2], "\n")
    end
    @test length(ambiguities) <= 10
end

const Z = Zero()
const I = One()

@testset "Promotion" begin
    @test promote_type(One, Bool) === Bool
    @test promote_type(One, Int32) === Int32
    @test promote_type(One, Int64) === Int64
    @test promote_type(One, Float32) === Float32
    @test promote_type(One, Float64) === Float64
    @test promote_type(One, Complex{Bool}) === Complex{Bool}
    @test promote_type(One, Complex{Float64}) === Complex{Float64}
    @test promote_type(Zero, Bool) === Bool
    @test promote_type(Zero, Int32) === Int32
    @test promote_type(Zero, Int64) === Int64
    @test promote_type(Zero, Float32) === Float32
    @test promote_type(Zero, Float64) === Float64
    @test promote_type(Zero, Complex{Bool}) === Complex{Bool}
    @test promote_type(Zero, Complex{Float64}) === Complex{Float64}
end

@testset "Real" begin
    @test iszero(Z) === true
    @test isone(I) === true
    @test isone(Z) === false
    @test iszero(I) === false
    @test Int(Z) === 0
    @test Int(I) === 1
    @test Float64(Z) === 0.0
    @test Float64(I) === 1.0
    @test float(Z) === 0.0
    @test float(I) === 1.0
    @test Z + 1 === 1
    @test I + 1 === 2
    @test Z + I === I
    @test Z - I === -1
    @test Zero(0) === Z
    @test One(1) === I
    @test one(Z) === I
    @test one(I) === I
    @test zero(Z) === Z
    @test zero(I) === Z
    @test Z === Zero()
    @test I === One()
    @test Z == Z
    @test I == I
    @test Z == 0
    @test I == 1
    @test sizeof(Z) == 0
    @test float(Z) === 0.0
    @test float(I) === 1.0
    @test -Z === Z
    @test -I === -1
    @test +Z === Z
    @test +I === I
    @test 2*Z === Z
    @test Z*Z === Z
    @test I*I === I
    @test 2*I === 2
    @test 2/I === 2.0
    @test I*2 === 2
    @test 2.0*Z === Z
    @test 2.0*I === 2.0
    @test Z*3 === Z
    @test Z/2 === Z
    @test Z/I === Z
    @test Z-Z === Z
    @test Z+Z === Z
    @test Z*Z === Z
    @test 1-Z === 1
    @test 1.0-Z === 1.0
    @test Z-1 == -1
    @test Z+1 == 1
    @test 2+Z == 2
    @test (Z < Z) == false
    @test (Z > Z) == false
    @test Z <= Z
    @test Z >= Z
    @test Z < 3
    @test Z > -2.0
    @test ldexp(Z, 3) === Z
    @test copysign(Z, 3) === Z
    @test copysign(Z, -1) === Z
    @test flipsign(Z, -1) === Z
    @test sign(Z) === Z
    @test round(Z) === Z
    @test floor(Z) === Z
    @test ceil(Z) === Z
    @test trunc(Z) === Z
    @test significand(Z) === Z
    @test sign(I) === I
    @test round(I) === I
    @test floor(I) === I
    @test ceil(I) === I
    @test trunc(I) === I
    @test significand(I) === I
    @test !isodd(Z)
    @test iseven(Z)
    @test isodd(I)
    @test !iseven(I)
    @test string(Z) == "𝟎"
    @test string(I) == "𝟏"
    @test mod(Z, 3) === Z
    @test mod(Z, 3.0) === Z
    @test mod(I, I) == mod(1, 1)
    @test mod(I, I) === Z
    @test rem(I, I) == rem(1, 1)
    @test rem(I, I) === Z
    @test rem(Z, 3) === Z
    @test rem(Z, 3.0) === Z
    @test modf(Z) === (Z, Z)
    @test I/I === I
    @test 1.0/I === 1.0
    @test 1/I === 1.0

    @test log(I) == log(1)
    @test log(I) === Z
    @test log2(I) == log2(1)
    @test log2(I) === Z
    @test log10(I) == log10(1)
    @test log10(I) === Z
    @test log1p(Z) == log1p(0)
    @test log1p(Z) === Z
    @test exp(Z) == exp(0)
    @test exp(Z) === I
    @test exp2(Z) == exp2(0)
    @test exp2(Z) === I
    @test exp10(Z) == exp10(0)
    @test exp10(Z) === I
    @test expm1(Z) == expm1(0)
    @test expm1(Z) === Z
    @test sin(Z) == sin(0)
    @test sin(Z) === Z
    @test cos(Z) == cos(0)
    @test cos(Z) === I
    @test sincos(Z) == sincos(0)
    @test sincos(Z) === (Z, I)
    @test tan(Z) == tan(0)
    @test tan(Z) === Z
    @test asin(Z) == asin(0)
    @test asin(Z) === Z
    @test atan(Z) == atan(0)
    @test atan(Z) === Z
    @test sinpi(Z) == sinpi(0)
    @test sinpi(Z) === Z
    @test sinpi(I) == sinpi(1)
    @test sinpi(I) === Z
    @test cospi(Z) == cospi(0)
    @test cospi(Z) === I
    @test sinh(Z) == sinh(0)
    @test sinh(Z) === Z
    @test cosh(Z) == cosh(0)
    @test cosh(Z) === I
    @test tanh(Z) == tanh(0)
    @test tanh(Z) === Z
    @test sqrt(Z) == sqrt(0)
    @test sqrt(Z) === Z
    @test sqrt(I) == sqrt(1)
    @test sqrt(I) === I
    @test cbrt(Z) == cbrt(0)
    @test cbrt(Z) === Z
    @test cbrt(I) == cbrt(1)
    @test cbrt(I) === I

    @test 3.0^Z === I
    @test 3.0f0^Z === I
    @test 3^Z === I
    @test 3^Z === I
    @test Z^Z === I
    @test I^Z === I

    @test I//I isa Rational{One}
    @test I//I + I//I == 2

    @test sum(fill(I, 23)) == 23
    @inferred sum(fill(I, 23))
    @test sum(One[]) === 0
    @test sum([One()]) === 1

    @test sum(fill(Z, 23)) === Z
    @inferred sum(fill(Z, 23))
    @test sum(Zero[]) === Z
    @test sum([Zero()]) === Z

    @test prod(fill(I, 23)) == 1
    @inferred prod(fill(I, 23))
    @test prod(One[]) === I
    @test prod([One()]) === I

    @test prod(fill(Z, 23)) === false
    @inferred prod(fill(Zero(), 23))
    @test prod(Zero[]) === true
    @test prod([Zero()]) === false

    @test false + One() === 1
    @test One() + false === 1

    @test false + Zero() === false
    @test Zero() + false === false

    @test Integer(Z) === 0
    @test Integer(I) === 1

    @test AbstractFloat(Z) === 0.0
    @test AbstractFloat(I) === 1.0

    @test Zero(Z) === Z
    @test One(I) === I

    @test inv(I) === I

    @test div(3, I) === 3
    @test div(Int32(3), I) === Int32(3)

    @test Bool(Z) === false
    @test Bool(I) === true
    @test Int32(Z) === Int32(0)
    @test Int32(I) === Int32(1)
    @test Int64(Z) === 0
    @test Int64(I) === 1
    @test Float32(Z) === 0.0f0
    @test Float32(I) === 1.0f0
    @test Float64(Z) === 0.0
    @test Float64(I) === 1.0

    @test Zero()^0 === true
    @test Zero()^1 === false

    # Bypass litteral_pow
    n0 = 0
    @test Zero()^n0 === true
    n1 = 1
    @test Zero()^n1 === false

    @test hypot(-7, Z) === 7
    @test hypot(Z, -7) === 7
    @test hypot(-7.0, Z) === 7.0
    @test hypot(Z, -7.0) === 7.0
    @test hypot(Z, Z) === Z

    @test I:3 === Base.OneTo(3)
end

# @testset "muladd" begin
#     @test fma(Z,1,Z) === Z
#     @test fma(Z,1.0,Z) === Z
#     @test muladd(Z,1,Z) === Z
#     @test muladd(Z,1.0,Z) === Z
#     @test fma(Z,1,3) === 3
#     @test fma(Z,1.0,3) === 3.0
#     @test muladd(Z,1,3) === 3
#     @test muladd(Z,1.0,3) === 3.0
#     @test fma(Z,Z,Z) === Z
#     @test muladd(Z,Z,Z) === Z
# end

@testset "Complex" begin
    @test Z*im === Z
    @test im*Z === Z
    @test real(Z) === Z
    @test imag(Z) === Z
    @test Z+(2+3im) === 2+3im
    @test (2+3im)+Z === 2+3im
    @test Z-(2+3im) === -2-3im
    @test (2+3im)-Z === 2+3im
    @test Z*(2+3im) === Z
    @test (2.0+3.0im)*Z === Z
    @test Z/(2+3im) === Z

    @test Complex(Z,Z) == Z
    @test Complex(1,Z) == 1
    @test Complex(1.0,Z) == 1.0
    @test Complex(true,Z) == true
end

@testset "testzero" begin
    @test testzero(3) === 3
    @test testzero(3+3im) === 3+3im
    @test testzero(0) === Z
    @test testzero(0+0im) === Z
end

@testset "error handling" begin
    @test_throws InexactError Zero(1)
    @test_throws InexactError convert(Zero, 0.1)
    @test_throws DivideError 1.0/Z
    @test_throws DivideError (1.0+2.0im)/Z
    @test_throws DivideError Z/Z
    @test_throws DivideError I/Z
    @test_throws DivideError inv(Z)
    @test_throws InexactError Zero(One())
    @test_throws InexactError One(Zero())
    @test_throws DivideError mod(Z, Z)
    @test_throws DivideError rem(Z, Z)
    @test_throws MethodError convert(Irrational{:π}, Zero())
    @test_throws MethodError convert(Irrational{:π}, One())
end

if VERSION ≥ v"1.2" # no test on 1.0 because of BenchmarkTools dependency
    # Test `MyComplex` example type
    include("complex2_example.jl")
    @testset "Complex example" begin
        @test Complex2(Zero(), 2)*Complex2(Zero(), 3) === Complex2(-6, Zero())
        @test abs(i) === One()
    end
end

Zeros.@pirate Base

@testset "piracy" begin
    @test +() === Z
    @test *() === I
    @test zero() === Z
    @test one() === I
    @test sum([]) == 0
    @test prod([]) == 1
    @test sin(π) == 0
    @test tan(π) == 0
end


text:
import UnicodePlots

plot2!(::Val{:UnicodePlots}, p::Union{Plot,Nothing}, X, Ys; kind, title, xlab, ylab, legend, legendpos, names, colors, xlim, ylim, ycat, xunit, yunit, aspect, width=40, height=15) = begin
    canvas = if get(ENV, "CI", nothing) == "true"
        UnicodePlots.DotCanvas
    else
        UnicodePlots.BrailleCanvas
    end

    if kind == :line || kind == :scatterline
        plot! = UnicodePlots.lineplot!
    elseif kind == :scatter
        plot! = UnicodePlots.scatterplot!
    elseif kind == :step
        plot! = UnicodePlots.stairs!
        #HACK: convert symbol array to integer index array
        ylim = (1, length(ycat))
        Ys = indexin.(Ys, Ref(ycat))
        height = length(ycat)
    else
        error("unrecognized plot kind = $kind")
    end

    !isnothing(legendpos) && @warn "unsupported legend position = $legendpos"

    create_colors(colors; n0=0) = begin
        n = length(colors)
        C = collect(Iterators.take(Iterators.cycle(UnicodePlots.color_cycle), n+n0))[n0+begin:end]
        f(c::Int, _) = p.opt[:colors][c]
        f(c, _) = c in keys(UnicodePlots.color_encode) ? c : :normal
        f(::Nothing, i) = C[i]
        [f(c, i) for (i, c) in enumerate(colors)]
    end

    #HACK: handle x-axis type of Date/DateTime adapted from UnicodePlots.lineplot()
    xlim_value(v::Dates.TimeType) = Dates.value(v)
    xlim_value(v) = v
    xlimval = xlim_value.(xlim)

    annotate_x_axis!(obj) = begin
        #HACK: override xlim string (for Date/DateTime)
        UnicodePlots.label!(obj, :bl, string(xlim[1]), color=:light_black)
        UnicodePlots.label!(obj, :br, string(xlim[2]), color=:light_black)
    end

    annotate_y_axis!(obj) = begin
        if kind == :step
            #HACK: show original string of symbol instead of index number
            for (i, y) in zip(height:-1:1, ycat)
                UnicodePlots.label!(obj, :l, i, string(y), color=:light_black)
            end
        end
    end

    if isnothing(p)
        a = Float64[]
        !isnothing(aspect) && (width = round(Int, aspect * 2height))
        obj = UnicodePlots.Plot(a, a, canvas; title, xlabel=xlab, ylabel=ylab, xlim=xlimval, ylim, width, height)
        UnicodePlots.label!(obj, :r, legend)
        annotate_x_axis!(obj)
        annotate_y_axis!(obj)
        p = Plot(obj; Xs=[], Ys=[], kinds=[], colors=[], title, xlab, ylab, legend, names, xlim, ylim, xunit, yunit, aspect, width, height)
    end
    colors = create_colors(colors; n0=length(p.opt[:Ys]))
    for (i, (Y, name)) in enumerate(zip(Ys, names))
        #HACK: UnicodePlots can't handle missing
        Y = coalesce.(Y, NaN)
        color = colors[i]
        plot!(p.obj, X, Y; name, color)
        #TODO: remember colors
        update!(p; Xs=[X], Ys=[Y], kinds=[kind], colors=[color])
    end
    p
end

plot3!(::Val{:UnicodePlots}, X, Y, Z; kind, title, legend, legendpos, xlab, ylab, zlab, xlim, ylim, zlim, zgap, zlabgap, aspect, width=0, height=30) = begin
    if kind == :heatmap
        ;
    elseif kind == :contour
        @warn "unsupported plot kind = $kind"
    else
        error("unrecognized plot kind = $kind")
    end

    !legend && @warn "unsupported legend = $legend"
    !isnothing(legendpos) && @warn "unsupported legend position = $legendpos"

    !isnothing(zgap) && @warn "unsupported contour interval = $zgap"
    !isnothing(zlabgap) && @warn "unsupported countour label interval = $zlabgap"

    arr(a) = sort(unique(a))
    x = arr(X)
    y = arr(Y)
    M = reshape(Z, length(y), length(x))

    offset(a) = a[1]
    xoffset = offset(x)
    yoffset = offset(y)
    fact(a) = (a[end] - offset(a)) / (length(a) - 1)
    xfact = fact(x)
    yfact = fact(y)

    !isnothing(aspect) && (width = round(Int, aspect * height))

    #TODO: support zlim (minz/maxz currentyl fixed in UnicodePlots)
    obj = UnicodePlots.heatmap(M; title, xlabel=xlab, ylabel=ylab, zlabel=zlab, xfact, yfact, xlim, ylim, xoffset, yoffset, width, height)
    Plot(obj; X, Y, Z, kind, title, xlab, ylab, zlab, xlim, ylim, zlim, aspect, width, height)
end


text:
import Yao: probs
export QCBM, QCBMGo!, psi, mmdgrad

include("Kernels.jl")

##### QCBM methods #####
struct QCBM{BT<:AbstractBlock, KT<:AbstractKernel}
    circuit::BT
    kernel::KT
    ptrain::Vector{Float64}

    dbs
end
function QCBM(circuit::AbstractBlock, kernel::AbstractKernel, ptrain::Vector)
    QCBM(circuit, kernel, ptrain, collect_blocks(Diff, circuit))
end

# INTERFACES
circuit(qcbm::QCBM) = qcbm.circuit
diff_blocks(qcbm::QCBM) = qcbm.dbs
"""
    loss(qcbm::QCBM, [p=qcbm|>probs]) -> Float64

loss function, optional parameter `p` is the probability distribution.
"""
function loss(qcbm::QCBM, p=qcbm|>probs)
    p -= qcbm.ptrain
    kernel_expect(qcbm.kernel, p, p)
end
"""
    gradient(qcbm::QCBM, p0=qcbm |> probs) -> Vector

gradient of MMD two sample test loss, `db` must be contained in qcbm.
optional parameter `p0` is current probability distribution.
"""
gradient(qcbm::QCBM, p0=qcbm |> probs) = mmdgrad.(Ref(qcbm), qcbm.dbs, p0=p0)

"""generated wave function"""
psi(qcbm) = zero_state(qcbm.circuit |> nqubits) |> qcbm.circuit
"""generated probability distribution"""
probs(qcbm::QCBM) = qcbm |> psi |> probs
"""
    mmdgrad(qcbm::QCBM, db::Diff; p0::Vector) -> Float64

gradient of MMD two sample test loss, `db` must be contained in qcbm.
`p0` is current probability distribution.
"""
function mmdgrad(qcbm::QCBM, db::Diff; p0::Vector)
    statdiff(()->probs(qcbm) |> as_weights, db, StatFunctional(kmat(qcbm.kernel)), initial=p0 |> as_weights) -
        2*statdiff(()->probs(qcbm) |> as_weights, db, StatFunctional(kmat(qcbm.kernel)*qcbm.ptrain))
end

"""
quantum circuit born machine trainer.
"""
struct QCBMGo!{QT<:QCBM, OT}
    qcbm::QT
    optimizer::OT
    niter::Int
end
QCBMGo!(qcbm::QCBM, optimizer) = QCBMGo!(qcbm, optimizer, typemax(Int))
Base.length(qo::QCBMGo!) = qo.niter

function Base.iterate(qo::QCBMGo!, state=(1, parameters(qo.qcbm.circuit)))
    state[1] > qo.niter && return nothing

    # initialize the parameters
    p0 = qo.qcbm |> probs
    grad = gradient(qo.qcbm, p0)
    QuAlgorithmZoo.update!(state[2], grad, qo.optimizer)
    dispatch!(qo.qcbm.circuit, state[2])
    Dict("probs"=>p0, "step"=>state[1], "gradient"=>grad), (state[1]+1, state[2])
end

"""
    gaussian_pdf(x, μ::Real, σ::Real)

gaussian probability density function.
"""
function gaussian_pdf(x, μ::Real, σ::Real)
    pl = @. 1 / sqrt(2pi * σ^2) * exp(-(x - μ)^2 / (2 * σ^2))
    pl / sum(pl)
end

@testset "qcbm" begin
    # problem setup
    n = 6
    depth = 6

    N = 1<<n
    kernel = rbf_kernel(0:N-1, 0.25)
    pg = gaussian_pdf(1:N, N/2-0.5, N/4)
    circuit = variational_circuit(n, depth, pair_ring(n)) |> autodiff(:QC)
    dispatch!(circuit, :random)
    qcbm = QCBM(circuit, kernel, pg)

    # training
    niter = 100
    optim = QuAlgorithmZoo.Adam(lr=0.1)
    for info in QCBMGo!(qcbm, optim, niter) end
    @test qcbm |> loss < 1e-4
end


text:
### ARSState

# ARSState holds the internal state ("local variables") of the acceptance-rejection sampler

type ARSState{S<:ValueSupport, F<:VariateForm} <: MCSamplerState{F}
  pstate::ParameterState{S, F} # Parameter state used internally by ARS
  tune::MCTunerState
  logproposal::Real
  weight::Real
end

ARSState{S<:ValueSupport, F<:VariateForm}(
  pstate::ParameterState{S, F},
  tune::MCTunerState,
  logproposal::Real,
  weight::Real
) =
  ARSState{S, F}(pstate, tune, logproposal, weight)

ARSState{S<:ValueSupport, F<:VariateForm}(
  pstate::ParameterState{S, F},
  tune::MCTunerState=BasicMCTune()
) =
  ARSState(pstate, tune, NaN, NaN)

### Acceptance-rejection sampler (ARS)

immutable ARS <: MCSampler
  logproposal::Function # Possibly unnormalized log-proposal (envelope)
  proposalscale::Real # Scale factor for logproposal (envelope) to cover target
  jumpscale::Real # Scale factor for adapting jump size

  function ARS(logproposal::Function, proposalscale::Real, jumpscale::Real)
    @assert jumpscale > 0 "Scale factor for adapting jump size is not positive"
    new(logproposal, proposalscale, jumpscale)
  end
end

ARS(logproposal::Function; proposalscale::Real=1., jumpscale::Real=1.) = ARS(logproposal, proposalscale, jumpscale)

### Initialize ARS sampler

## Initialize parameter state

function initialize!{F<:VariateForm}(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::ARS,
  outopts::Dict
)
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array(Any, length(pstate.diagnostickeys))
  end
end

## Initialize ARSState

sampler_state{S<:ValueSupport, F<:VariateForm}(
  parameter::Parameter{S, F},
  sampler::ARS,
  tuner::MCTuner,
  pstate::ParameterState{S, F},
  vstate::VariableStateVector
) =
  ARSState(generate_empty(pstate), tuner_state(parameter, sampler, tuner))

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Univariate},
  x::Real,
  parameter::Parameter{Continuous, Univariate},
  sampler::ARS
)
  pstate.value = x
  parameter.logtarget!(pstate)
end

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::ARS
)
  pstate.value = copy(x)
  parameter.logtarget!(pstate)
end

reset!{S<:ValueSupport, F<:VariateForm}(
  sstate::ARSState{S, F},
  pstate::ParameterState{S, F},
  parameter::Parameter{S, F},
  sampler::MCSampler,
  tuner::MCTuner
) =
  reset!(sstate.tune, sampler, tuner)

Base.show(io::IO, sampler::ARS) =
  print(io, "ARS sampler: proposal scale = $(sampler.proposalscale), jump scale = $(sampler.jumpscale)")

Base.writemime(io::IO, ::MIME"text/plain", sampler::ARS) = show(io, sampler)


text:
"""
Composition-Rejection with Rejection sampling method (RSSA-CR)
"""

const MINJUMPRATE = 2.0^exponent(1e-12)

mutable struct RSSACRJumpAggregation{F,U,S,F1,F2,RNG,VJMAP,JVMAP,BD,T2V,P<:PriorityTable,W<:Function} <: AbstractSSAJumpAggregator
    next_jump::Int
    prev_jump::Int
    next_jump_time::F
    end_time::F
    cur_rate_low::Vector{F}
    cur_rate_high::Vector{F}
    sum_rate::F
    ma_jumps::S
    rates::F1
    affects!::F2
    save_positions::Tuple{Bool,Bool}
    rng::RNG
    vartojumps_map::VJMAP
    jumptovars_map::JVMAP
    bracket_data::BD
    ulow::T2V
    uhigh::T2V
    minrate::F
    maxrate::F   # initial maxrate only, table can increase beyond it!
    rt::P #rate table
    ratetogroup::W
    cur_u_bnds::Matrix{U} # current bounds on state u
  end

function RSSACRJumpAggregation(nj::Int, njt::F, et::F, crs::Vector{F}, sum_rate::F, maj::S, rs::F1,
                               affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG; u::U,
                               vartojumps_map=nothing, jumptovars_map=nothing,
                               bracket_data=nothing, minrate=convert(F,MINJUMPRATE),
                               maxrate=convert(F,Inf), kwargs...) where {F,S,F1,F2,RNG,U}
    # a dependency graph is needed and must be provided if there are constant rate jumps
    if vartojumps_map === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use the RSSACR algorithm a map from variables to dependent jumps must be supplied.")
        else
            vtoj_map = var_to_jumps_map(length(u), maj)
        end
    else
        vtoj_map = vartojumps_map
    end

    if jumptovars_map === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use the RSSACR algorithm a map from jumps to dependent variables must be supplied.")
        else
            jtov_map = jump_to_vars_map(maj)
        end
    else
        jtov_map = jumptovars_map
    end

    # vectors to store bracketing intervals for jump rates
    crl_bnds = similar(crs)
    crh_bnds = similar(crs)

    # a bracket data structure is needed for updating species populations
    bd = (bracket_data === nothing) ? BracketData{F,eltype(U)}() : bracket_data

    # matrix to store bracketing interval for species and the relative interval width
    # first row is Xlow, second is Xhigh
    cs_bnds = Matrix{eltype(U)}(undef, 2, length(u))
    ulow    = @view cs_bnds[1,:]
    uhigh   = @view cs_bnds[2,:]

    # mapping from jump rate to group id
    minexponent = exponent(minrate)

    # use the largest power of two that is <= the passed in minrate
    minrate = 2.0^minexponent
    ratetogroup = rate -> priortogid(rate, minexponent)

    # construct an empty initial priority table -- we'll reset this in init
    rt = PriorityTable(ratetogroup, zeros(F, 1), minrate, 2*minrate)

    RSSACRJumpAggregation{typeof(njt),eltype(U),S,F1,F2,RNG,typeof(vtoj_map),typeof(jtov_map),typeof(bd),typeof(ulow),typeof(rt),typeof(ratetogroup)}(
                            nj, nj, njt, et, crl_bnds, crh_bnds, sum_rate, maj, rs, affs!, sps, rng, vtoj_map,
                            jtov_map, bd, ulow, uhigh, minrate, maxrate, rt, ratetogroup, cs_bnds)
end

############################# Required Functions ##############################

# creating the JumpAggregation structure (function wrapper-based constant jumps)
function aggregate(aggregator::RSSACR, u, p, t, end_time, constant_jumps,
                   ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using function wrappers
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)

    build_jump_aggregation(RSSACRJumpAggregation, u, p, t, end_time, ma_jumps,
                           rates, affects!, save_positions, rng; u=u, kwargs...)
end

# set up a new simulation and calculate the first jump / jump time
function initialize!(p::RSSACRJumpAggregation, integrator, u, params, t)
    set_bracketing!(p,u,params,t)

    # setup PriorityTable
    reset!(p.rt)
    for (pid,priority) in enumerate(p.cur_rate_high)
        insert!(p.rt, pid, priority)
    end

    generate_jumps!(p, integrator, u, params, t)
    nothing
end

# execute one jump, changing the system state
function execute_jumps!(p::RSSACRJumpAggregation, integrator, u, params, t)
    # execute jump
    u = update_state!(p, integrator, u)

    # update rates
    update_dependent_rates!(p, u, params, t)
    nothing
end

# calculate the next jump / jump time
function generate_jumps!(p::RSSACRJumpAggregation, integrator, u, params, t)
    sum_rate = p.sum_rate
    # if no more events possible there is nothing to do
    if nomorejumps!(p, sum_rate)
        return nothing
    end

    @unpack rt, ma_jumps, rates, cur_rate_high, cur_rate_low, rng = p
    num_majumps = get_num_majumps(ma_jumps)
    rerl = zero(sum_rate)

    jidx = sample(rt, cur_rate_high, rng)
    if iszero(jidx)
        p.next_jump_time = Inf
        return nothing
    end
    rerl += randexp(rng)
    while rejectrx(ma_jumps, num_majumps, rates, cur_rate_high, cur_rate_low, rng, u, jidx, params, t)
        # sample candidate reaction
        jidx  = sample(rt, cur_rate_high, rng)
        rerl += randexp(rng)
    end
    p.next_jump = jidx

    # update time to next jump
    p.next_jump_time = t + rerl / sum_rate
    nothing
end

######################## SSA specific helper routines #########################
"update bracketing for species that depend on the just executed jump"
@inline function update_dependent_rates!(p::RSSACRJumpAggregation, u, params, t)
    # update bracketing intervals
    ubnds       = p.cur_u_bnds
    crhigh      = p.cur_rate_high

    @inbounds for uidx in p.jumptovars_map[p.next_jump]
        uval = u[uidx]
        # if new u value is outside the bracketing interval
        if uval == zero(uval) || uval < ubnds[1,uidx] || uval > ubnds[2,uidx]
            # update u bracketing interval
            ubnds[1,uidx], ubnds[2,uidx] = get_spec_brackets(p.bracket_data, uidx, uval)

            # for each dependent jump, update jump rate brackets
            for jidx in p.vartojumps_map[uidx]
                oldrate = crhigh[jidx]
                p.cur_rate_low[jidx], crhigh[jidx] = get_jump_brackets(jidx, p, params, t)

                # update the priority table
                update!(p.rt, jidx, oldrate, crhigh[jidx])
            end
        end
    end

    p.sum_rate = groupsum(p.rt)
    nothing
  end


text:
module PorousModels

using .._Diff: derivative

using ArgCheck: @argcheck

include("base.jl")
include("common.jl")
include("brooksandcorey.jl")
include("vangenuchten.jl")

export UnsaturatedFlowModel, hθ, θh, Ch, Kh, Cθ, Kθ, Dθ
export BrooksAndCorey
export VanGenuchten

end


text:
using TestSetExtensions
using Suppressor

using Compat.Test

orig_color = Base.have_color

eval(Base, :(have_color = true))
output_color = @capture_out begin
    @testset ExtendedTestSet "top-level tests" begin
        @testset "2nd-level tests 1" begin
            @test true
            @test 1 == 1
        end
        @testset "2nd-level tests 2" begin
            @test true
            @test 1 == 1
        end
    end
end

eval(Base, :(have_color = false))
output_nocolor = @capture_out begin
    @testset ExtendedTestSet "top-level tests" begin
        @testset "2nd-level tests 1" begin
            @test true
            @test 1 == 1
        end
        @testset "2nd-level tests 2" begin
            @test true
            @test 1 == 1
        end
    end
end

eval(Base, :(have_color = $orig_color))

try
    info("You should see 3 failing tests with pretty diffs...")
    include(joinpath("..", "diffdemo.jl"))
catch
end
try
    info("These 4 failing tests don't have pretty diffs to display")
    @testset ExtendedTestSet "not-pretty" begin
        @testset "No pretty diff for matrices" begin
            @test [1 2; 3 4] == [1 4; 3 4]
        end
        @testset "don't diff non-equality" begin
            @test 1 > 2
        end
        @testset "don't diff non-comparisons" begin
            @test iseven(7)
        end
        @testset "errors don't have diffs either" begin
            throw(ErrorException("This test is supposed to throw an error"))
        end
    end
catch
end

@testset ExtendedTestSet "TextSetExtensions Tests" begin
    @testset "check output" begin
        if VERSION <= v"0.6.0-"
            @test split(output_color, '\n')[1] == "\e[1m\e[32m.\e[0m\e[1m\e[32m.\e[0m\e[1m\e[32m.\e[0m\e[1m\e[32m.\e[0m"
        else
            @test split(output_color, '\n')[1] == "\e[32m.\e[39m\e[32m.\e[39m\e[32m.\e[39m\e[32m.\e[39m"
        end
        @test split(output_nocolor, '\n')[1] == "...."
    end

    @testset "DottedTestSet is deprecated" begin
        @test_warn "DottedTestSet is deprecated, use ExtendedTestSet instead." @testset DottedTestSet "testing" begin
            @test true
        end
    end

    @testset "Auto-run test files" begin
        global file1_run = false
        global file2_run = false
        global file3_run = false

        @includetests

        @test file1_run
        @test file2_run
        @test file3_run
    end

    @testset "run selected test files" begin
        global file1_run = false
        global file2_run = false
        global file3_run = false

        @includetests ["file1", "file3"]

        @test file1_run
        @test !file2_run
        @test file3_run
    end

    @testset "more than one arg to @includetests is an error" begin
        ex = macroexpand(:(@includetests one two))
        @test ex.head == :error
    end
end


text:
# This file is a part of SIS3316Digitizers.jl, licensed under the MIT License (MIT).

using StruckVMEDevices.MemRegisters
using Test

using Tables, TypedTables

@testset "registers" begin
    @enum SomeIntEnum::UInt8 begin
        SIE_FOO = 0
        SIE_BAR = 1
        SIE_BAZ = 2
    end

    @enum SomeBoolEnum::Bool begin
        SBE_FOO = false
        SBE_BAR = true
    end

    shift_by_one = (x -> signed(x + 1), x -> unsigned(x) - 1)
    bool_to_3_4 = (x -> (x ? 4 : 3), x -> Bool(x - 3))

    @test bits_to_val(0x42, LiteralBits()) == 0x42
    @test val_to_bits(0x42, LiteralBits()) == 0x42

    @test bits_to_val(0x42, TransformedBits(shift_by_one...)) == 67
    @test val_to_bits(67, TransformedBits(shift_by_one...)) == 0x42

    @test bits_to_val(0x01, EnumBits(SomeIntEnum)) == SIE_BAR
    @test val_to_bits(SIE_BAR, EnumBits(SomeIntEnum)) == 0x01

    @test bits_to_val(true, EnumBits(SomeBoolEnum)) == SBE_BAR
    @test val_to_bits(SBE_BAR, EnumBits(SomeBoolEnum)) == true

    @test getraw(0x5a, Bit{RW}(1)) == true
    @test setraw(0x50, Bit{RW}(1), true) == 0x52

    @test getval(0x5a, Bit{RO}(1)) == true
    @test getval(0x5a, Bit{RO}(1, bool_to_3_4)) == 4
    @test getval(0x5a, Bit{RO}(1, SomeBoolEnum)) == SBE_BAR

    @test setval(0x50, Bit{WO}(1), true) == 0x52
    @test setval(0x50, Bit{WO}(1, bool_to_3_4), 4) == 0x52
    @test setval(0x50, Bit{WO}(1, SomeBoolEnum), SBE_BAR) == 0x52

    @test getraw(0xa5, BitRange{RO}(1:2)) == 0x02
    @test setraw(0xa0, BitRange{WO}(1:2), 0x02) == 0xa4

    @test getval(0xa5, BitRange{RW}(1:2)) == 0x02
    @test getval(0xa5, BitRange{RW}(1:2, shift_by_one)) == 3
    @test getval(0xa5, BitRange{RW}(1:2, SomeIntEnum)) == SIE_BAZ

    @test setval(0x50, BitRange{RW}(1:2), 0x02) == 0x54
    @test setval(0x50, BitRange{RW}(1:2, shift_by_one), 3) == 0x54
    @test setval(0x50, BitRange{RW}(1:2, SomeIntEnum), SIE_BAZ) == 0x54

    @test setraw(BitWriteOperation{UInt32}(0, 0, true), BitRange{JKRW}(4:7), 0x05) == BitWriteOperation{UInt32}(0x00a00050, 0xffffffff, true)
    @test setraw(BitWriteOperation{UInt32}(0x00a00050, 0xffffffff, true), Bit{JKRW}(1), true) == BitWriteOperation{UInt32}(0x00a00052, 0xffffffff, true)
    @test setraw(BitWriteOperation{UInt32}(0x00a00050, 0xffffffff, true), Bit{JKRW}(1), false) == BitWriteOperation{UInt32}(0x00a20050, 0xffffffff, true)
    @test_throws ArgumentError setraw(BitWriteOperation{UInt32}(0x0, 0x1, true), BitRange{JKRW}(4:7), 0x05)
    @test_throws ArgumentError setraw(BitWriteOperation{UInt32}(0, 0, false), BitRange{JKRW}(4:7), 0x05)

    @test @inferred(merge(
        BitWriteOperation{UInt32}(0x000000ff, 0x0f00ffff, false),
        BitWriteOperation{UInt32}(0xf000f000, 0xf000f00f, false)
    )) == BitWriteOperation{UInt32}(0xf000f0f0, 0xff00ffff, false)

    @test @inferred(merge(
        BitWriteOperation{UInt32}(0xff0000ff, 0xffffffff, true),
        BitWriteOperation{UInt32}(0x000ff000, 0xffffffff, true)
    )) == BitWriteOperation{UInt32}(0x0f0ff0f0, 0xffffffff, true)

    nbits = (a = Bit{RO}(5), b = BitRange{WO}(4:5, SomeIntEnum))

    @test @inferred(Register{RW,UInt32}(0x08, nbits)) isa Register
    @test isbits(Register{RW,UInt32}(0x08, nbits))

    reg = Register{RW,UInt32}(0x08, nbits)

    @test @inferred((r -> (;r...))(reg)) == nbits

    @test @inferred((x -> x.a)(reg)) == BitSelRef{UInt32}(0x08, Bit{RO}(5))

    @test @inferred(keys(reg)) == (:a, :b)
    @test reg[:b] == BitSelRef{UInt32}(0x08, BitRange{WO}(4:5, SomeIntEnum))

    @test @inferred(filter(Readable, getlayout(reg))) == (a = Bit{RO}(5),)
    @test @inferred(filter(Writeable, getlayout(reg))) == (b = BitRange{WO}(4:5, SomeIntEnum),)

    reg2 = Register{RW,UInt16}(UInt32(0x08), (a = Bit{RW}(3, SomeBoolEnum), b = BitRange{RW}(4:7, SomeIntEnum)));
end


text:
using MLJ
import RDatasets: dataset
import DataFrames: DataFrame, select, Not, describe
using Random

data = dataset("datasets", "USArrests")
names(data)

describe(data, :mean, :std)

X = select(data, Not(:State))
X = coerce(X, :UrbanPop=>Continuous, :Assault=>Continuous);

PCA = @load PCA pkg=MultivariateStats

pca_mdl = PCA(pratio=1)
pca = machine(pca_mdl, X)
fit!(pca)
PCA
W = MLJ.transform(pca, X);

schema(W).names

r = report(pca)
cumsum(r.principalvars ./ r.tvar)

data = dataset("ISLR", "OJ")

feature_names = [
    :PriceCH, :PriceMM, :DiscCH, :DiscMM, :SalePriceMM, :SalePriceCH,
    :PriceDiff, :PctDiscMM, :PctDiscCH
]

X = select(data, feature_names);

Random.seed!(1515)

SPCA = Pipeline(
    Standardizer(),
    PCA(pratio=1-1e-4)
)

spca = machine(SPCA, X)
fit!(spca)
W = MLJ.transform(spca, X)
names(W)

rpca = collect(values(report(spca).report_given_machine))[2]
cs = cumsum(rpca.principalvars ./ rpca.tvar)

using PyPlot

figure(figsize=(8,6))

PyPlot.bar(1:length(cs), cs)
plot(1:length(cs), cs, color="red", marker="o")

xlabel("Number of PCA features", fontsize=14)
ylabel("Ratio of explained variance", fontsize=14)

Random.seed!(1515)

KMeans = @load KMeans pkg=Clustering
SPCA2 = Pipeline(
    Standardizer(),
    PCA(),
    KMeans(k=3)
)

spca2 = machine(SPCA2, X)
fit!(spca2)

assignments = collect(values(report(spca2).report_given_machine))[3].assignments
mask1 = assignments .== 1
mask2 = assignments .== 2
mask3 = assignments .== 3;

using PyPlot

figure(figsize=(8, 6))
for (m, c) in zip((mask1, mask2, mask3), ("red", "green", "blue"))
    plot(W[m, 1], W[m, 2], ls="none", marker=".", markersize=10, color=c)
end

xlabel("PCA-1", fontsize=13)
ylabel("PCA-2", fontsize=13)
legend(["Group 1", "Group 2", "Group 3"], fontsize=13)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# ml.jl
# in future


text:
module EnvConfig

export greet2, greet1

using Pkg;
# Pkg.resolve()
# Pkg.add(PackageSpec(url="https://github.com/DennisRutjes/Binance.jl",rev="master"))

packages=["Dates","DataFrames","JSON"]
Pkg.add(packages)

using Dates
using DataFrames
using JSON

function greet1()
        a = 3
    println("hello from greet1")
end

function greet2()

    function greet3()
        println("greet3 start")
        greet1()
        println("hello from greet3")
    end

    a = 3
    println("hello from greet2")
    greet3()
end

# greet() = print("Hello World!")

println("greetings from module env_config")
end # module

println("greetings from top level env_config")


text:
function beta_equation(cur_reactor::AbstractReactor)
  cur_G_T = K_TB(cur_reactor)

  cur_G_T *= cur_reactor.T_bar

  cur_equation = Equation(1, 1, 0, cur_G_T)

  cur_equation
end


text:
"""
    function mech_and_cow()

# Example

```jldoctest
julia> cowsay("Do you ever get that feeling...?", cow=Cowsay.mech_and_cow)
 __________________________________
< Do you ever get that feeling...? >
 ----------------------------------
   \\                               ,-----.
   /                               |     |
   \\                            ,--|     |-.
   /                     __,----|  |     | |
   \\                   ,;::     |  `_____' |
   /                   `._______|    i^i   |
   \\                            `----| |---'| .
   /                       ,-------._| |== ||//
   \\                       |       |_|P`.  /'/
   /                       `-------' 'Y Y/'/'
   \\                                 .==\\ /_\\
   ^__^                             /   /'|  `i
   (oo)\\_______                   /'   /  |   |
   (__)\\       )\\/\\            /'    /   |   `i
      ||----w||           ___,;`----'.___L_,-'`\\__
      ||     ||          i_____;----\\.____i""\\____\\
```
"""
function mech_and_cow(;eyes="oo", tongue="  ", thoughts="\\")
sthguoht = thoughts;
if thoughts == "\\"
    sthguoht = "/"
end #if

the_cow =
"""
   $thoughts                               ,-----.
   $sthguoht                               |     |
   $thoughts                            ,--|     |-.
   $sthguoht                     __,----|  |     | |
   $thoughts                   ,;::     |  `_____' |
   $sthguoht                   `._______|    i^i   |
   $thoughts                            `----| |---'| .
   $sthguoht                       ,-------._| |== ||//
   $thoughts                       |       |_|P`.  /'/
   $sthguoht                       `-------' 'Y Y/'/'
   $thoughts                                 .==\\ /_\\
   ^__^                             /   /'|  `i
   ($eyes)\\_______                   /'   /  |   |
   (__)\\       )\\/\\            /'    /   |   `i
    $tongue||----w||           ___,;`----'.___L_,-'`\\__
      ||     ||          i_____;----\\.____i""\\____\\
"""
return the_cow
end


text:
"""
    aov1(x, Y)
    Univariate anova test.
* `x` : Univariate categorical X-data.
* `Y` : Y-data.
"""
function aov1(x, Y)
    Y = ensure_mat(Y)
    n = length(x)
    q = size(Y, 2)
    A = length(unique(x))
    Xdummy = dummy(x).Y
    zY = center(Y, colmeans(Y))
    fm = mlr(Xdummy, zY) ;
    pred = predict(fm, Xdummy).pred
    SSF = sum((pred.^2), dims = 1)   # = colvars(pred) * n
    SSR = ssr(pred, zY)
    df_fact = A - 1
    df_res = n - A
    MSF = SSF / df_fact
    MSR = SSR / df_res
    F = MSF ./ MSR
    d = Distributions.FDist(df_fact, df_res)
    pval = Distributions.ccdf(d, F)
    (SSF = SSF, SSR = SSR, df_fact = df_fact, df_res = df_res, F = F, pval = pval)
end


text:
"""
    ZeroOneBridge{T}

The `ZeroOneBridge` splits a `MOI.VariableIndex`-in-`MOI.ZeroOne` constraint
into a `MOI.VariableIndex`-in-`MOI.Integer` constraint
and a `MOI.VariableIndex`-in-`MOI.Interval(0, 1)` constraint.
"""
struct ZeroOneBridge{T} <: AbstractBridge
    interval_index::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{T}}
    integer_index::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}
end

function bridge_constraint(
    ::Type{ZeroOneBridge{T}},
    model::MOI.ModelLike,
    f::MOI.VariableIndex,
    ::MOI.ZeroOne,
) where {T<:Real}
    interval_index =
        MOI.add_constraint(model, f, MOI.Interval{T}(zero(T), one(T)))
    integer_index = MOI.add_constraint(model, f, MOI.Integer())
    return ZeroOneBridge{T}(interval_index, integer_index)
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:ZeroOneBridge{T}},
) where {T}
    return Tuple{Type,Type}[
        (MOI.VariableIndex, MOI.Interval{T}),
        (MOI.VariableIndex, MOI.Integer),
    ]
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:ZeroOneBridge})
    return Tuple{Type}[]
end

function concrete_bridge_type(
    ::Type{<:ZeroOneBridge{T}},
    ::Type{MOI.VariableIndex},
    ::Type{MOI.ZeroOne},
) where {T}
    return ZeroOneBridge{T}
end

function MOI.supports_constraint(
    ::Type{<:ZeroOneBridge},
    ::Type{MOI.VariableIndex},
    ::Type{MOI.ZeroOne},
)
    return true
end

MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, ::ZeroOneBridge) = MOI.ZeroOne()

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::ZeroOneBridge,
)
    return MOI.get(model, attr, bridge.interval_index)
end

function MOI.delete(model::MOI.ModelLike, bridge::ZeroOneBridge)
    MOI.delete(model, bridge.interval_index)
    MOI.delete(model, bridge.integer_index)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal},
    bridge::ZeroOneBridge,
)
    return MOI.get(model, attr, bridge.interval_index)
end

function MOI.get(
    ::ZeroOneBridge{T},
    ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.Interval{T}},
)::Int64 where {T}
    return 1
end

function MOI.get(
    ::ZeroOneBridge,
    ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.Integer},
)::Int64
    return 1
end

function MOI.get(
    bridge::ZeroOneBridge,
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Interval{T}},
) where {T}
    return [bridge.interval_index]
end

function MOI.get(
    bridge::ZeroOneBridge,
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer},
)
    return [bridge.integer_index]
end