text:
#### LeveSpec types for various types of level specification

## TODO: Heads => True

abstract type LevelSpec end

struct LevelSpecToDepth  <: LevelSpec
    level::Int
end

struct LevelSpecAtDepth  <: LevelSpec
    level::Int
end

struct LevelSpecRange  <: LevelSpec
    start::Int
    stop::Int
end

struct LevelSpecAll  <: LevelSpec
    includezero::Bool
end

LevelSpecAll() = LevelSpecAll(false)

#### Create LevelSpec instances from Mxpr expressions

## Interpretation of negative index is different here than from Mma
## In Mma, a negative index counts up from each leaf.
## TODO: or FIXME
## 1. Implement both depth-first and breadth-first traversal. Only
## breadth-first is implmented now.
## 2. For negative indicies. Traverse the whole tree and build a flat
## Array of of the depth of each leaf. Then count negative levels from
## each leaf.

# Form n.  levels 1 through n
function make_level_specification(expr,n::Integer)
    if n < 0 n += depth(expr) end
    LevelSpecToDepth(n)
end

# Form [n] or [m,n]  level n only, or levels m through n
function make_level_specification(expr,spec::Mxpr{:List})
    _make_level_specification(expr,spec,margs(spec)...)
end

# Form Infinity or else error. levels 1 through Infinity
function make_level_specification(expr,x)
    if x == Infinity     # could use singlton type
        return LevelSpecRange(1,depth(expr))
#        return LevelSpecAll(false)
    end
    symerror("Level specification $x is not of the form n, {n}, or {m, n}.")
end

# Form [n]
function _make_level_specification(expr,spec::Mxpr{:List}, n::Integer)
    if n < 0 n += depth(expr) end
    LevelSpecAtDepth(n)
end

# Form [m,n]
function _make_level_specification(expr, spec::Mxpr{:List}, start::Integer, stop::Integer)
    if start < 0 start += depth(expr) end
    if stop < 0 stop += depth(expr) end
    LevelSpecRange(start, stop)
end

function _make_level_specification(expr, spec::Mxpr{:List}, start::Integer, stop)
    stop != Infinity && symerror("Level specification $spec is not of the form n, [n], or [m, n].")
    d = depth(expr)
    if start < 0 start += d end
    LevelSpecRange(start,d)
end

function _make_level_specification(expr, spec::Mxpr{:List},args...)
    symerror("Level specification $spec is not of the form n, [n], or [m, n].")
end

#### Structure for performing action at levels

mutable struct LevelAction
    data::Any            # application-specific data
    doaction::Function   # function takes arguments (data,expr)
    levelind::Int
    subind::Int
    parent::Any
    levelbreak::Bool
end

LevelAction(data,doaction) = LevelAction(data,doaction,0,0,Null,false)

has_level_zero(spec::LevelSpecAtDepth) = spec.level == 0
has_level_zero(spec::LevelSpecToDepth) = spec.level == 0
has_level_zero(spec::LevelSpecRange) = spec.start == 0
has_level_zero(spec::LevelSpecAll) = spec.includezero

# This is length for mapping purposes
function length_for_level(x)
    return 0
end

function length_for_level(mx::Mxpr)
    return length(mx)
end

function checkbreak(action)
    if action.levelbreak
        true
    else
        false
    end
end

# macro travcode()
#     code = quote
#         elen = length_for_level(expr)
#         if elen > 0
#             action.parent = expr
#             for i in 1:elen
#                 action.parent = expr
#                 action.levelind += 1
#                 action.subind = i

#                 traverse_levels!(action,spec,expr[i])
#                 if  checkbreak(action)
#                     return
#                 end
#                 action.levelind -= 1
#             end
#             #            action.parent = Null
#         end
#     end
#     return code
# end

const traversecode =
    quote
        elen = length_for_level(expr)
        if elen > 0
            action.parent = expr
            for i in 1:elen
                action.parent = expr
                action.levelind += 1
                action.subind = i

                traverse_levels!(action,spec,expr[i])
                if  checkbreak(action)
                    return
                end
                action.levelind -= 1
            end
            #            action.parent = Null
        end
    end

@eval function traverse_levels!(action::LevelAction, spec::LevelSpecAtDepth, expr)
    if action.levelind == spec.level
        action.doaction(action.data, expr)
        action.levelbreak && return
    elseif action.levelind < spec.level
        $traversecode
    end
end

@eval function traverse_levels!(action::LevelAction, spec::LevelSpecToDepth, expr)
    if action.levelind <= spec.level && action.levelind > 0
        action.doaction(action.data, expr)
        checkbreak(action) && return
    end
    if action.levelind < spec.level
        $traversecode
    end
end

@eval function traverse_levels!(action::LevelAction, spec::LevelSpecRange, expr)
    if action.levelind <= spec.stop && action.levelind >= spec.start
        action.doaction(action.data, expr)
        action.levelbreak && return
    end
    if action.levelind < spec.stop
        $traversecode
    end
end

@eval function traverse_levels!(action::LevelAction, spec::LevelSpecAll, expr)
    action.doaction(action.data, expr)
    action.levelbreak && return
    $traversecode
end


text:
using Test
using Forecast

@testset "sma" begin

    @test_throws AssertionError sma(1:5,10)

    x = sma(1:5,3)
    @test length(x) == 3
    @test x == [2.0, 3.0, 4.0]

    x = sma(1:5,3,true)
    @test length(x) == 5
    @test length(collect(skipmissing(x))) == 3
    @test findlast(!ismissing,x) == 4
    @test collect(skipmissing(x)) == [2.0, 3.0, 4.0]

    x = sma(1:5,3,false)
    @test length(x) == 5
    @test length(collect(skipmissing(x))) == 3
    @test findlast(!ismissing,x) == 3
    @test collect(skipmissing(x)) == [2.0, 3.0, 4.0]

    r = rand(100)
    x = sma(r,10)
    y = collect(skipmissing(sma(r,10,false)))
    z = collect(skipmissing(sma(r,10,true)))
    @test x == y == z

end

ext: jl


text:
# module MySQL_Octo.SQL

function table_from(adapter::QueryBuilders.AbstractAdapter, tablename::String)
    meta = SQL.prepare(adapter, QueryBuilders.FromClause(tablename))
    fields = tuple(map(Symbol, meta.names)...)
    Octo.Adapters.Base.build_table(adapter, tablename, nothing, fields)
end


text:
module PackageTry

export f

function f()
x = 1
y = 2
return x + y
end

end


text:
"""
    Master

"""
mutable struct Master

    #===========================================================================
        RMP data
    ===========================================================================#

    rmp::MOI.ModelLike  # Restricted Master Problem

    # Constraints information
    con_cvx::Vector{MOI.ConstraintIndex}  # Convexity constraints in RMP
    con_link::Vector{MOI.ConstraintIndex}  # Linking constraints in RMP

    rhs_link::Vector{Float64}  # Right-hand side of linking constraints

    π::Vector{Float64}  # Shadow prices (linking constraints)
    σ::Vector{Float64}  # Shadow marginal costs (convexity constraints)

    # Columns information
    var_link::Vector{MOI.VariableIndex}  # Linking variables (may be artificial)
    columns::Vector{MOI.VariableIndex}  # Indices of RMP variables

    #===========================================================================
        Column-Generation
    ===========================================================================#

    #=
        Master Problem Status

        The Master Problem is:
        - Feasible if a feasible solution has been found in the RMP.
        - Optimal if solved to optimality, i.e. there exists a primal
            feasible solution, and pricing fails to identify a column of
            negative reduced cost.
        - Infeasible if proven infeasible, i.e. RMP is infeasible and
            Farkas pricing fails to cut the infeasibility ray.
        - Unbounded if the RMP is unbounded
    =#
    mp_status::MOI.TerminationStatusCode  # Status of MasterProblem
    mp_gap::Float64

    primal_lp_bound::Float64  # Primal bound for the linear (DW) relaxation
    primal_ip_bound::Float64  # Primal bound for the integer problem
    dual_bound::Float64   # Lagrange dual bound
    dual_bound_estimate::Float64  # Estimate of the Lagrange dual bound

    #===========================================================================
        Constructor
    ===========================================================================#

    function Master(
        rmp::MOI.ModelLike,
        con_cvx::Vector{<:MOI.ConstraintIndex},
        con_link::Vector{<:MOI.ConstraintIndex},
        rhs_link::Vector{Float64},
        var_link::Vector{MOI.VariableIndex},
        columns::Vector{MOI.VariableIndex}
    )
        # Check indices
        mcvx = length(con_cvx)
        mlink = length(con_link)
        for cidx in con_cvx
            MOI.is_valid(rmp, cidx) || error("Invalid constraint index $cidx")
        end
        for cidx in con_link
            MOI.is_valid(rmp, cidx) || error("Invalid constraint index $cidx")
        end
        length(rhs_link) == mlink || error(
            "There are $mlink constraints but RHS has length $(length(rhs_link))"
        )

        # Check variable indices
        for vidx in var_link
            MOI.is_valid(rmp, vidx) || error("Invalid variable index $vidx")
        end
        for vidx in columns
            MOI.is_valid(rmp, vidx) || error("Invalid variable index $vidx")
        end

        # Instanciate model
        mp = new(
            rmp,
            con_cvx, con_link, rhs_link,
            zeros(mlink), zeros(mcvx),
            var_link, columns,
            MOI.OPTIMIZE_NOT_CALLED, Inf, Inf, Inf, -Inf, -Inf
        )

        return mp
    end
end

"""
    add_column!(mp, col)

Add column `col` to Master Problem `mp`.

All columns are non-negative by default.
"""
function add_column!(mp::Master, col::Column)

    # add column to RMP
    vidx, _ = MOI.add_constrained_variable(mp.rmp, MOI.GreaterThan(0.0))

    # Set objective coefficient
    MOI.modify(mp.rmp,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarCoefficientChange(vidx, col.cost)
    )

    # Set column coefficient for convexity constraint
    if col.is_vertex
        MOI.modify(mp.rmp, mp.con_cvx[col.spind], MOI.ScalarCoefficientChange(vidx, 1.0))
    end
    # Set column coefficients for linking constraints
    for (i, val) in zip(col.rind, col.rval)
        MOI.modify(mp.rmp, mp.con_link[i], MOI.ScalarCoefficientChange(vidx, val))
    end

    # Book-keeping
    push!(mp.columns, vidx)

    # Done
    return vidx
end


text:
using Parameters, Random
import Distributions: logpdf, loglikelihood

struct Paired{T1,T2} <: ContinuousUnivariateDistribution
    d::T1
    parms::T2
end

Paired(;d, parms) = Paired(d, parms)

loglikelihood(d::Paired, data::Array{<:NamedTuple,1}) = logpdf(d, data)

function logpdf(d::Paired, data::Array{<:NamedTuple,1})
    LL = computeLL(data, d.parms; d=d.d)
    return LL
end

function sample_stimuli(stimuli, N)
    return stimuli[1:N]
end

function simulate(all_stimuli, fixed_parms, n_blocks, n_trials, deadline=5.0, isi=5.0; d)
    # sample unique stimului for the simulation
    stimuli = sample_stimuli(all_stimuli, n_trials)
    # initialize declarative memory
    chunks = [Chunk(word=:_,number=1)] |> empty!
    memory = Declarative(;memory=chunks)
    # create act-r object
    actr = ACTR(;declarative=memory, fixed_parms..., d)
    # preallocate data
    data = Array{NamedTuple,1}()
    # iterate through all blocks
    for block in 1:n_blocks
        # randomize stimulus order
        shuffle!(stimuli)
        temp = simulate_block(actr, stimuli, block, deadline, isi)
        push!(data, temp...)
    end
    return vcat(data...)
end

function simulate_block(actr, stimuli, block, deadline=5.0, isi=5.0)
    # initialize data for the current block
    data = Array{NamedTuple,1}()
    start_time = get_time(actr)
    for stimulus in stimuli
        temp = simulate_trial(actr, stimulus, block, deadline)
        push!(data, temp)
        # 5 second response deadline + additional 5 seconds
        start_time += (deadline + isi)
        set_time!(actr, start_time)
    end
    return data
end

function simulate_trial(actr, stimulus, block, deadline=5.0)
    N=0;L=0.0;recent=Float64[];retrieved=:failed;time_created=0.0
    trial_start = get_time(actr)
    # encode and conflict resolution time
    e_time = 0.085 + 0.05 + 0.05
    add_time!(actr, e_time)
    # retrieve paired associate based on word
    chunk = retrieve(actr; word=stimulus.word)
    # retrieval time
    r_time = compute_RT(actr, chunk)
    # add encoding time and retrieval time to reaction time
    rt = e_time + r_time + .05 + 0.300
    # interrupt retrieval if exceeds deadline
    if rt > deadline
        # set to truncated value
        rt = deadline
        # advance current time to deadline
        set_time!(actr, trial_start + deadline)
        retrieved = :truncated
        # retrieval info if no chunk is retrieved
        chunk = get_chunks(actr; stimulus...)
        if !isempty(chunk)
            # get chunk state at time of retrieval
            @unpack N,L,recent,time_created = chunk[1]
            recent = copy(recent)
        end
    elseif !isempty(chunk)
        # add retreival time to current time
        add_time!(actr, r_time)
        # next conflict resolution
        add_time!(actr, .05)
        retrieved = :retrieved
        # get chunk state at time of retrieval
        @unpack N,L,recent,time_created = chunk[1]
        recent = copy(recent)
        # add chunk if new, or update N and time stamps
        add_chunk!(actr; stimulus...)
    else
        # add retreival failure time to current time
        add_time!(actr, r_time)
        # next conflict resolution
        add_time!(actr, .05)
        # retrieval info if no chunk is retrieved
        chunk = get_chunks(actr; stimulus...)
        if !isempty(chunk)
            # get chunk state at time of retrieval
            @unpack N,L,recent,time_created = chunk[1]
            recent = copy(recent)
        end
    end
    # feedback is always presented after the deadline
    set_time!(actr, trial_start + deadline)
    # encode feedback
    add_time!(actr, 0.05 + 0.05 + 0.085)
    # add chunk: create new chunk if does not exist, otherwise increment N and add time stamp
    add_chunk!(actr; stimulus...)
    data = (stimulus=stimulus,N=N,L=L,block=block,time_created=time_created,
        recent=recent,rt=rt,retrieved=retrieved)
    return data
end

function computeLL(data, fixed_parms; d)
    chunk = Chunk(act=zero(d))
    # initialize declarative memory
    memory = Declarative(;memory=[chunk])
    # create ACTR object with declarative memory and parameters
    actr = ACTR(;declarative=memory, fixed_parms..., d)
    # extract parameters
    @unpack s,lf,τ,ter=actr.parms
    # do not add noise activation values
    actr.parms.noise = false
    σ = s * pi / sqrt(3)
    # initialize log likelihood
    LL = 0.0
    for k in data
        cur_time = k.L + k.time_created
        if k.retrieved == :retrieved
            # reproduce the memory state right before retrieval
            modify!(chunk; N=k.N, recent=k.recent, time_created=k.time_created)
            # compute activation
            compute_activation!(actr, chunk, cur_time)
            # adjust activation values for latency factor
            v = [chunk.act,τ] .- log(lf)
            # log normal race distribution object
            dist = LNR(-v, σ, ter)
            # compute log likelihood of retrieving chunk at time k.rt
            LL += logpdf(dist, 1, k.rt)
        elseif k.retrieved == :truncated
            # special case in which chunk has not been created yet
            if k.N == 0
                v = τ .- log(lf)
                # compute log likelihood of failing to respond within deadline
                LL += logccdf(LogNormal(-v, σ), 5)
            else
                # reproduce the memory state right before retrieval
                modify!(chunk; N=k.N, recent=k.recent, time_created=k.time_created)
                # compute activation
                compute_activation!(actr, chunk, cur_time)
                # adjust activation values for latency factor
                v = [chunk.act,τ] .- log(lf)
                # compute log likelihood of failing to respond within deadline
                LL += logccdf(LogNormal(-v[1], σ), 5) + logccdf(LogNormal(-v[2], σ), 5)
            end
        else
            # retrieval failure for novel pair
            if k.N == 0
                dist = LNR(-[τ-log(lf)], σ, ter)
                LL += logpdf(dist, 1, k.rt)
            else
                # reproduce the memory state immediately prior to retrieval
                modify!(chunk; N=k.N, recent=k.recent, time_created=k.time_created)
                # compute activation
                compute_activation!(actr, chunk, cur_time)
                # mean activation values
                v = [chunk.act,τ] .- log(lf)
                # log normal race distribution object
                dist = LNR(-v, σ, ter)
                # compute log likelihood of failing to retrieve chunk at time k.rt
                LL += logpdf(dist, 2, k.rt)
            end
        end
    end
    return LL
end

function show_learning(all_stimuli, parms, n_blocks=8, n_trials=20; d)
    data = simulate(all_stimuli, parms, n_blocks, n_trials; d)
    df = DataFrame(data)
    groups = groupby(df, [:block,:retrieved])
    return combine(groups, :rt=>mean=>:rt)
end


text:
## Kalman filter for univariate linear Gaussian models

struct KalmanOut
  predictionMeans::Vector{Float64}
  predictionVariances::Vector{Float64}
  filteringMeans::Vector{Float64}
  filteringVariances::Vector{Float64}
  smoothingMeans::Vector{Float64}
  smoothingVariances::Vector{Float64}
  logZhats::Vector{Float64}
end

function lognormpdf(y::Float64, μ::Float64, σ²::Float64)
  v::Float64 = y-μ
  lnc::Float64 = -0.5*log(2*π*σ²)
  return lnc - 1/(2*σ²)*v*v
end

function kalman(theta::LGTheta, ys::Vector{Float64})
  n = length(ys)
  A = theta.A
  C = theta.C
  Q = theta.Q
  R = theta.R
  x0 = theta.x0
  v0 = theta.v0
  predictionMeans = Vector{Float64}(undef, n)
  predictionVariances = Vector{Float64}(undef, n)
  filteringMeans = Vector{Float64}(undef, n)
  filteringVariances = Vector{Float64}(undef, n)
  logZhats = Vector{Float64}(undef, n)
  mutt1 = 0.0
  mutt = 0.0
  sigmatt1 = 0.0
  sigmatt = 0.0
  lZ = 0.0
  for p = 1:n
    if p == 1
      mutt1 = x0
      sigmatt1 = v0
    else
      mutt1 = A*mutt
      sigmatt1 = A*sigmatt*A + Q
    end
    predictionMeans[p] = mutt1
    predictionVariances[p] = sigmatt1
    lZ += lognormpdf(ys[p], C*mutt1, C*sigmatt1*C + R)
    logZhats[p] = lZ
    K = sigmatt1 * C / (C * sigmatt1 * C + R)
    mutt = mutt1 + K * (ys[p]-C*mutt1)
    sigmatt = sigmatt1 - K * C * sigmatt1
    filteringMeans[p] = mutt
    filteringVariances[p] = sigmatt
  end
  smoothingMeans = Vector{Float64}(undef, n)
  smoothingVariances = Vector{Float64}(undef, n)
  smoothingMeans[n] = filteringMeans[n]
  smoothingVariances[n] = filteringVariances[n]
  for p = n:-1:2
    J = filteringVariances[p-1] * A * inv(predictionVariances[p])
    smoothingMeans[p-1] = filteringMeans[p-1] +
      J * (smoothingMeans[p] - predictionMeans[p])
    smoothingVariances[p-1] = filteringVariances[p-1] +
      J * (smoothingVariances[p] - predictionVariances[p]) * J
  end
  return KalmanOut(predictionMeans, predictionVariances, filteringMeans,
        filteringVariances, smoothingMeans, smoothingVariances, logZhats)
end

function kalmanlogZ(theta::LGTheta, ys::Vector{Float64})
  n = length(ys)
  A = theta.A
  C = theta.C
  Q = theta.Q
  R = theta.R
  x0 = theta.x0
  v0 = theta.v0
  mutt1 = 0.0
  mutt = 0.0
  sigmatt1 = 0.0
  sigmatt = 0.0
  lZ = 0.0
  for p = 1:n
    if p == 1
      mutt1 = x0
      sigmatt1 = v0
    else
      mutt1 = A*mutt
      sigmatt1 = A*sigmatt*A + Q
    end
    lZ += lognormpdf(ys[p], C*mutt1, C*sigmatt1*C + R)
    K = sigmatt1 * C / (C * sigmatt1 * C + R)
    mutt = mutt1 + K * (ys[p]-C*mutt1)
    sigmatt = sigmatt1 - K * C * sigmatt1
  end
  return lZ
end


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

mutable struct TrafficSelectorPolicy <: SwaggerModel
    localAddressRanges::Any # spec type: Union{ Nothing, Vector{String} } # spec name: localAddressRanges
    remoteAddressRanges::Any # spec type: Union{ Nothing, Vector{String} } # spec name: remoteAddressRanges

    function TrafficSelectorPolicy(;localAddressRanges=nothing, remoteAddressRanges=nothing)
        o = new()
        validate_property(TrafficSelectorPolicy, Symbol("localAddressRanges"), localAddressRanges)
        setfield!(o, Symbol("localAddressRanges"), localAddressRanges)
        validate_property(TrafficSelectorPolicy, Symbol("remoteAddressRanges"), remoteAddressRanges)
        setfield!(o, Symbol("remoteAddressRanges"), remoteAddressRanges)
        o
    end
end # type TrafficSelectorPolicy

const _property_map_TrafficSelectorPolicy = Dict{Symbol,Symbol}(Symbol("localAddressRanges")=>Symbol("localAddressRanges"), Symbol("remoteAddressRanges")=>Symbol("remoteAddressRanges"))
const _property_types_TrafficSelectorPolicy = Dict{Symbol,String}(Symbol("localAddressRanges")=>"Vector{String}", Symbol("remoteAddressRanges")=>"Vector{String}")
Base.propertynames(::Type{ TrafficSelectorPolicy }) = collect(keys(_property_map_TrafficSelectorPolicy))
Swagger.property_type(::Type{ TrafficSelectorPolicy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TrafficSelectorPolicy[name]))}
Swagger.field_name(::Type{ TrafficSelectorPolicy }, property_name::Symbol) =  _property_map_TrafficSelectorPolicy[property_name]

function check_required(o::TrafficSelectorPolicy)
    (getproperty(o, Symbol("localAddressRanges")) === nothing) && (return false)
    (getproperty(o, Symbol("remoteAddressRanges")) === nothing) && (return false)
    true
end

function validate_property(::Type{ TrafficSelectorPolicy }, name::Symbol, val)
end


text:
using Distributions
using StatsBase
using Interpolations
using ATA
using DataFrames
using DelimitedFiles
using JuMP
using Cbc
using CSV
# @everywhere using Pkg; Pkg.activate(".")

# @everywhere using Psychometrics

#using Pkg; Pkg.activate(".")
@everywhere using Psychometrics

cd("examples/bootstrap")
function bs()
    I_total = 250
    N = 3_000
    metric = [0.0, 1.0]

    #pre-test settings
    T = 6
    test_length = 50
    nDifficulty = 3 #levels of guessed difficulty (low, medium, high)
    alphas = [0.25, 0.75] #difficulty cuts

    ### generate true items

    a_dist = LogNormal(0, 0.25); # mu= exp(0.4 + (0.5*0.1^2)), sigma^2= exp(0.4 + 0.1^2)*(exp(0.1^2) - 1)
    a_bounds = [0.5, 5.0];
    b_dist = Normal(0, 1);
    b_bounds = [-4.0, 4.0];
    items = [Item(i, string("item_", i), ["math"], Parameters2PL(Product([a_dist, b_dist]), a_bounds, b_bounds)) for i = 1:I_total];

    ## generate true examinees

    latent_dist = Normal(0.0, 1.0)
    latent_bounds = [-4.0, 4.0]
    examinees = Psychometrics.Examinee[]
    for n in 1:N
        push!(examinees,  Examinee(n, string("examinee_", n), Latent1D(latent_dist, latent_bounds)))
    end

    ## pre-test assembly  (needs ATA package)
    # generate guessed difficulties @everywhere using true difficulties and quantiles
    bQles = quantile(map(i -> i.parameters.b, items), alphas)
    guessed_difficulty = zeros(I_total)
    for i=1:I_total
        if items[i].parameters.b < bQles[1]
            guessed_difficulty[i]=1
        elseif items[i].parameters.b < bQles[2]
            guessed_difficulty[i]=2
        else
            guessed_difficulty[i]=3
        end
    end
    #save true bank
    item_bank = DataFrame(
        a = map(i -> i.parameters.a, items),
        b = map(i -> i.parameters.b, items),
        guessed_difficulty = string.(Int.(guessed_difficulty))
    )
    CSV.write("item_bank.csv", item_bank)

    op = [ 0 10 0 0 0 0;
          10 0 10 0 0 0;
          0 10 0 10 0 0;
          0 0 10 0 10 0;
          0 0 0 10 0 10;
          0 0 0 0 10 0 ]
    DelimitedFiles.writedlm("op.csv", op, ',')

    #write constraints and save in csv
    constraints = DataFrame(
        group = [1, 1, 1],
        var = ["guessed_difficulty", "guessed_difficulty", "guessed_difficulty"],
        value = [1, 2, 3],
        min = [
            Int(floor(alphas[1]*test_length - 2)),
            Int(floor((alphas[2] - alphas[1])*test_length - 2)),
            Int(floor((1 - alphas[2])*test_length - 2))
            ],
        max = [
            Int(ceil(alphas[1]*test_length + 2)),
            Int(ceil((alphas[2] - alphas[1])*test_length + 2)),
            Int(ceil((1 - alphas[2])*test_length + 2))
            ],
        weight = [1.0, 1.0, 1.0]
    )
    CSV.write("constraints.csv", constraints)

    ata_model = compact_ata(
        settings_file = "pre_test_settings.jl",
        bank_file = "item_bank.csv",
        bank_delim = ",",
        add_friends = false,
        add_enemies = false,
        add_constraints = true,
        constraints_file = "constraints.csv",
        constraints_delim = ",",
        add_overlap = true,
        overlap_file = "op.csv",
        overlap_delim = ",",
        add_exp_score = false,
        group_by_friends = false,
        add_obj_fun = false,
        solver = "jumpATA",
        print_it = true,
        print_folder = "RESULTS",
        plot_it = false,
        optimizer_constructor = "Cbc",
        optimizer_attributes = [("seconds", 500), ("logLevel", 1)]
    )

    ## generate responses, use ata_model.output.design to assign items to the examinees
    # 1:T examinees get the 1:T test forms,
    responses = Response[]
    t = 1
    for n = 1:N
        for i in items[findall(ata_model.output.design[:, t] .> 0)]
            push!(responses, answer(examinees[n], i))
        end
        t += 1
        if t > T
            t = 1
        end
    end

    #start points and probs
    prob = Distributions.pdf(Normal(metric[1], metric[2]), collect(-6.0:0.3:6.0))
    prob = prob ./ sum(prob)
    dist = Distributions.DiscreteNonParametric(collect(-6.0:0.3:6.0), prob)

    # initialize examinees
    latent_est_prior = dist
    latent_est_bounds = [-4.0, 4.0]
    examinees_est = [Examinee(e, string("examinee_",e), Latent1D(latent_est_prior, latent_est_bounds)) for e = 1 : N];
    #set starting value taking a random value in the interval -0.5:0.5
    map(e ->
        begin
            e.latent.val = 0.0
        end,
    examinees_est)

    #initalize items estimates
    a_est_bounds = [1e-5, 5.0];
    b_est_bounds = [-4.0, 4.0];

    items_est = [Item(i, string("item_", i), ["math"], Parameters2PL(a_est_bounds, b_est_bounds)) for i = 1 : I_total];
    map( i -> begin
        i.parameters.calibrated = false
    end, items_est);
    # starting values
    map( i -> begin
        i.parameters.a = 1.0
    end, items_est);
    map( i -> begin
        i.parameters.b = 0.0
    end, items_est);

    #start calibration
    joint_estimate!(
        items_est,
        examinees_est,
        responses;
        method = "mmle",
        dist = dist,
        metric = metric,
        max_iter = 500,
        max_time = 500,
        x_tol_rel = 0.0001,
        super_fast_2pl_1d = true,
    );

        println("a RMSE")
        println(sqrt(StatsBase.mean(map( (i, i_est) -> (i.parameters.a-i_est.parameters.a)^2, items, items_est))))
        println("b RMSE")
        println(sqrt(StatsBase.mean(map( (i, i_est) -> (i.parameters.b-i_est.parameters.b)^2, items, items_est))))
        println("latent RMSE")
        println(sqrt(StatsBase.mean(map( (i, i_est) -> (i.latent.val-i_est.latent.val)^2, examinees, examinees_est))))

    items_est_bs = map(i -> i, items_est)
    examinees_est_bs = map(e -> e, examinees_est)

    bootstrap!(
        items_est_bs,
        examinees_est_bs,
        responses;
        method = "mmle",
        dist = dist,
        metric = metric,
        max_iter = 500,
        max_time = 100,
        x_tol_rel = 0.01,
        f_tol_rel = 0.0001,
        replications = 500,
        type = "nonparametric",
        sample_fraction = 1.0,
        super_fast_2pl_1d = true,
        )
    return examinees, examinees_est, items, items_est, responses, items_est_bs, examinees_est_bs
end
examinees, examinees_est, items, items_est, responses, items_est_bs, examinees_est_bs = bs();
@everywhere using JLD2
@save items_est_bs "items_est_bs.jld2"


text:
struct FunctionMap{T, F1, F2} <: LinearMap{T}
    f::F1
    fc::F2
    M::Int
    N::Int
    _ismutating::Bool
    _issymmetric::Bool
    _ishermitian::Bool
    _isposdef::Bool
end
function FunctionMap{T}(f::F1, fc::F2, M::Int, N::Int;
    ismutating::Bool  = _ismutating(f),
    issymmetric::Bool = false,
    ishermitian::Bool = (T<:Real && issymmetric),
    isposdef::Bool    = false) where {T, F1, F2}
    FunctionMap{T, F1, F2}(f, fc, M, N, ismutating, issymmetric, ishermitian, isposdef)
end

# additional constructors
FunctionMap{T}(f, M::Int; kwargs...) where {T}         = FunctionMap{T}(f, nothing, M, M; kwargs...)
FunctionMap{T}(f, M::Int, N::Int; kwargs...) where {T} = FunctionMap{T}(f, nothing, M, N; kwargs...)
FunctionMap{T}(f, fc, M::Int; kwargs...) where {T}     = FunctionMap{T}(f, fc, M, M; kwargs...)

# show
function Base.show(io::IO, A::FunctionMap{T, F, Nothing}) where {T, F}
    print(io, "LinearMaps.FunctionMap{$T}($(A.f), $(A.M), $(A.N); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))")
end
function Base.show(io::IO, A::FunctionMap{T}) where {T}
    print(io, "LinearMaps.FunctionMap{$T}($(A.f), $(A.fc), $(A.M), $(A.N); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))")
end

# properties
Base.size(A::FunctionMap) = (A.M, A.N)
LinearAlgebra.issymmetric(A::FunctionMap) = A._issymmetric
LinearAlgebra.ishermitian(A::FunctionMap) = A._ishermitian
LinearAlgebra.isposdef(A::FunctionMap)    = A._isposdef
ismutating(A::FunctionMap) = A._ismutating
_ismutating(f) = first(methods(f)).nargs == 3

# multiplication with vector
function Base.:(*)(A::FunctionMap, x::AbstractVector)
    length(x) == A.N || throw(DimensionMismatch())
    if ismutating(A)
        y = similar(x, promote_type(eltype(A), eltype(x)), A.M)
        A.f(y, x)
    else
        y = A.f(x)
    end
    return y
end

function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
    (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch())
    ismutating(A) ? A.f(y, x) : copyto!(y, A.f(x))
    return y
end

function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
    issymmetric(A) && return A_mul_B!(y, A, x)
    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch())
    if A.fc !== nothing
        if !isreal(A)
            x = conj(x)
        end
        (ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x)))
        if !isreal(A)
            conj!(y)
        end
        return y
    else
        error("transpose not implemented for $A")
    end
end

function Ac_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
    ishermitian(A) && return A_mul_B!(y, A, x)
    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch())
    if A.fc !== nothing
        return (ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x)))
    else
        error("adjoint not implemented for $A")
    end
end


text:
include("../examples/tent.jl")
using JLD
σ_o = 0.1
function p_y_g_x(a)
		return log(1/sqrt(2π)/σ_o) -0.5*a*a/σ_o/σ_o
end
function obs(x)
    return x
end
function sir(y,x,s,N_thr=10)
	K = size(y)[1]
	N_p = size(x)[1]
    w = ones(N_p)./N_p
	new_pts = similar(x)
	for k = 1:K
    	for i = 1: N_p
				log(w[i]) = log(w[i]) + p_y_g_x(y[k] - obs(x[i]))
    	end
		@show sum(w)
		w ./= sum(w)
	    N_eff = 1.0/sum(w.*w)
        if (N_eff < N_thr)
	        #resample
            # multinomial resampling
			cdfw = cumsum(w)
			new_pts .= similar(x)
            for j = 1:N_p
			    r = rand()
			    for i = 1:N_p
		    	    if cdfw[i] >= r
					    new_pts[j] = x[i]
			            break
			        end
		        end
		    end
			x .= new_pts
	    end
		x .= next.(x,s)
    end
	return x, w
end
function assimilate()
    # Often much larger.
	Np = 100
	K = 10
	x = rand(Np)
	s = 0.001
	x_true = ones(K)
	x0_true = 2*pi*rand()
	x_true[1] = x0_true
	y = zeros(K)
	for i = 2:K
	    x_true[i] = next(x_true[i-1], s)
		y[i] = x_true[i] + σ_o*randn()
	end
	N_thr = 20
	# store trajectory of w and x.
	x, w = sir(y, x, s, N_thr)
	return x, w, x_true[K]
end


text:
__parse(::Type{T}, data::AbstractString) where {T} = convert(T, data)

__parse(::Type{T}, data::AbstractString) where {T<:Number} = parse(T, data)

__parse(::Type{T}, data::AbstractString) where {T<:Symbol} = T(data)

__parse(::Type{Vector{Symbol}}, data::Vector) = Symbol.(data)

__parse(::Type{Dict{Symbol, T1}}, data::Dict{String, T2}
       ) where {T1<:AbstractFloat, T2} =
    Dict(Symbol(k)=>T1(v) for (k,v) in data)

__parse(::Type{T}, data::S) where{T,S} = try
    T(data)
catch
    convert(T, data)
end

safe_symbol_eval(input_symbol, default_symbol) = begin
    if isdefined(@__MODULE__, input_symbol)
        return eval(input_symbol)
    else
        @warn "$input_symbol is not defined, defaulting to $default_symbol."
        return eval(default_symbol)
    end
end

unzip(it; n=length(it), offset=0, ndims=1) = begin
    from = offset+1
    to = min(length(it), from+n-1)
    map(i->getindex(getindex.(it, i), from:to), 1:ndims)
end


text:
mutable struct BitBuffer{B} <: AbstractVector{UInt}
    len::Int
    buffer::Vector{UInt8}

    function BitBuffer{B}(len, buffer) where B
        if !(B isa Int)
            throw(TypeError(:BitBuffer, Int, typeof(B)))
        end
        if !((B > 0) & (B < 8*sizeof(UInt) - 6))
            throw(ArgumentError("B must be in 1:$(8*sizeof(UInt) - 7)"))
        end
        new(len, buffer)
    end
end

bitsize(n::BitBuffer{B}) where B = B
bytes(n, b::Int) = cld((n-1)*b, 8) + sizeof(UInt)
BitBuffer{B}(::UndefInitializer, n) where B = BitBuffer{B}(n, Vector{UInt8}(undef, bytes(n, B)))
BitBuffer{B}(n) where B = BitBuffer{B}(n, zeros(UInt8, bytes(n, B)))
Base.copy(x::BitBuffer) = typeof(x)(x.len, copy(x.buffer))
Base.size(x::BitBuffer) = (x.len,)
Base.lastindex(x::BitBuffer) = x.len
Base.IndexStyle(x::BitBuffer) = IndexLinear()
@inline Base.checkbounds(::Type{Bool}, x::BitBuffer, i) = (i ≥ firstindex(x)) & (i ≤ lastindex(x))

function byteindex(ui::Unsigned, x::BitBuffer)
    bits = unsigned((ui-1) * bitsize(x))
    return Core.Intrinsics.udiv_int(bits, UInt(8)) + 1
end

function bitoffset(ui::Unsigned, x::BitBuffer)
    bits = unsigned((ui-1) * bitsize(x))
    return Core.Intrinsics.urem_int(bits, UInt(8))
end

bitmask(x::BitBuffer) = (UInt(1) << bitsize(x)) - 1

@inline function loadbits(x::BitBuffer, i::Unsigned)
    buffer = x.buffer
    GC.@preserve buffer begin
        bits = unsafe_load(Ptr{UInt}(pointer(buffer, i)))
    end
    return bits
end

@inline function unsafe_getindex(x::BitBuffer, i)
    index = byteindex(unsigned(i), x)
    bits = loadbits(x, index)
    offset = bitoffset(unsigned(i), x)
    return (bits >>> offset) & bitmask(x)
end

function Base.getindex(x::BitBuffer, i::Int)
    @boundscheck checkbounds(x, i)
    unsafe_getindex(x, i)
end

function unsafe_setindex!(x::BitBuffer, v::UInt, i)
    index = byteindex(unsigned(i), x)
    bits = loadbits(x, index)
    offset = bitoffset(unsigned(i), x)
    mask = ~(bitmask(x) << offset)
    bits &= mask
    bits |= v << offset
    buffer = x.buffer
    GC.@preserve buffer begin
        unsafe_store!(Ptr{UInt}(pointer(buffer, index)), bits)
    end
end

@noinline function throw_indexerr(::Type{BitBuffer{B}}, v::Integer) where B
    throw(ArgumentError("$v has more bits than maximum $B"))
end

function Base.setindex!(x::BitBuffer, v::Integer, i::Int)
    @boundscheck checkbounds(x, i)
    vu = convert(UInt, unsigned(v))
    vu > bitmask(x) && throw_indexerr(typeof(x), v)
    unsafe_setindex!(x, vu, i)
end


text:
# # Wiring diagrams in Compose.jl
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/graphics/composejl_wiring_diagrams.ipynb)
#
# Catlab can draw wiring diagrams using the Julia package
# [Compose.jl](https://github.com/GiovineItalia/Compose.jl).
#
# For best results, it is recommended to load the packages
# [Convex.j](https://github.com/JuliaOpt/Convex.jl) and
# [SCS.jl](https://github.com/JuliaOpt/SCS.jl). When available they are used to
# optimize the layout of the outer ports.

using Catlab.WiringDiagrams, Catlab.Graphics

import Convex, SCS

# ## Examples

# ### Symmetric monoidal category

using Catlab.Doctrines

A, B, C, D = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C, :D)
f, g = Hom(:f, A, B), Hom(:g, B, A);

# To start, here are a few very simple examples.

to_composejl(f)
#-
to_composejl(f⋅g)
#-
to_composejl(f⊗g)

# Here is a more complex example, involving generators with compound domains and
# codomains.

h, k = Hom(:h, C, D),  Hom(:k, D, C)
m, n = Hom(:m, B⊗A, A⊗B), Hom(:n, D⊗C, C⊗D)
q = Hom(:l, A⊗B⊗C⊗D, D⊗C⊗B⊗A)

to_composejl((f⊗g⊗h⊗k)⋅(m⊗n)⋅q⋅(n⊗m)⋅(h⊗k⊗f⊗g))

# Identities and braidings appear as wires.

to_composejl(id(A))
#-
to_composejl(braid(A,B))
#-
to_composejl(braid(A,B) ⋅ (g⊗f) ⋅ braid(A,B))

# The isomorphism $A \otimes B \otimes C \to C \otimes B \otimes A$ induced by
# the permutation $(3\ 2\ 1)$ is a composite of braidings and identities.

σ = (braid(A,B) ⊗ id(C)) ⋅ (id(B) ⊗ braid(A,C) ⋅ (braid(B,C) ⊗ id(A)))

to_composejl(σ)

# By default, anchor points are added along identity and braiding wires to
# reproduce the expression structure in the layout. The anchors can be disabled
# to get a more "unbiased" layout.

to_composejl(σ, anchor_wires=false)

# ### Biproduct category

A, B, C = Ob(FreeBiproductCategory, :A, :B, :C)
f = Hom(:f, A, B)

to_composejl(mcopy(A))
#-
to_composejl(delete(A))
#-
to_composejl(mcopy(A)⋅(f⊗f)⋅mmerge(B))
#-
to_composejl(mcopy(A⊗B), orientation=TopToBottom)
#-
to_composejl(mcopy(A⊗B⊗C), orientation=TopToBottom)

# ### Compact closed category

# The unit and co-unit of a compact closed category appear as caps and cups.

A, B = Ob(FreeCompactClosedCategory, :A, :B)

to_composejl(dunit(A))
#-
to_composejl(dcounit(A))

# In a self-dual compact closed category, such as a bicategory of relations,
# every morphism $f: A \to B$ has a transpose $f^\dagger: B \to A$ given by
# bending wires:

A, B = Ob(FreeBicategoryRelations, :A, :B)
f = Hom(:f, A, B)

to_composejl((dunit(A) ⊗ id(B)) ⋅ (id(A) ⊗ f ⊗ id(B)) ⋅ (id(A) ⊗ dcounit(B)))

# ### Abelian bicategory of relations

# In an abelian bicategory of relations, such as the category of linear
# relations, the duplication morphisms $\Delta_X: X \to X \otimes X$ and
# addition morphisms $\blacktriangledown_X: X \otimes X \to X$ belong to a
# bimonoid. Among other things, this means that the following two morphisms are
# equal.

X = Ob(FreeAbelianBicategoryRelations, :X)

to_composejl(plus(X) ⋅ mcopy(X))
#-
to_composejl((mcopy(X)⊗mcopy(X)) ⋅ (id(X)⊗braid(X,X)⊗id(X)) ⋅ (plus(X)⊗plus(X)))

# ## Custom styles

# The visual appearance of wiring diagrams can be customized by passing Compose
# [properties](http://giovineitalia.github.io/Compose.jl/latest/gallery/properties/).

using Compose: fill, stroke

A, B, = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

to_composejl(f⋅g, props=Dict(
  :box => [fill("lavender"), stroke("black")],
))
#-
X = Ob(FreeAbelianBicategoryRelations, :X)

to_composejl(plus(X) ⋅ mcopy(X), props=Dict(
  :junction => [fill("red"), stroke("black")],
  :variant_junction => [fill("blue"), stroke("black")],
))

# By default, the boxes are rectangular (`:rectangle`). Other available shapes
# include circles (`:circle`) and ellipses (`:ellipse`).

to_composejl(f⋅g, box_shape=:circle)

# ## Output formats

# The function `to_composejl` returns a `ComposePicture` object, which contains
# a Compose.jl context as well as a recommended width and height. When displayed
# interactively, this object is rendered using Compose's SVG backend.

# Any backend can be used by calling Compose's `draw` function. The SVG and
# [PGF](https://ctan.org/pkg/pgf) (LaTeX) backends are always available. To use
# the PNG or PDF backends, the extra packages
# [Cairo.jl](https://github.com/JuliaGraphics/Cairo.jl) and
# [Fontconfig.jl](https://github.com/JuliaGraphics/Fontconfig.jl) must be
# installed.

# For example, here is how to use the PGF backend.

using Compose: draw, PGF

pic = to_composejl(f⋅g, rounded_boxes=false)
pgf = sprint() do io
  pgf_backend = PGF(io, pic.width, pic.height,
    false, # emit_on_finish
    true,  # only_tikz
    texfonts=true)
  draw(pgf_backend, pic.context)
end
println(pgf)


text:
# Autogenerated wrapper script for Imath_jll for i686-w64-mingw32
export libImath

JLLWrappers.@generate_wrapper_header("Imath")
JLLWrappers.@declare_library_product(libImath, "libImath-3_1.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libImath,
        "bin\\libImath-3_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools

@testset "basic functionality" begin
    A = 2 * rand(ComplexF64, (20, 10)) .- 1
    v = rand(ComplexF64, 10)
    u = rand(ComplexF64, 20, 1)
    V = rand(ComplexF64, 10, 3)
    W = rand(ComplexF64, 20, 3)
    α = rand()
    β = rand()
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(M)

    @testset "LinearMaps.jl" begin
        @test eltype(M) == eltype(A)
        @test size(M) == size(A)
        @test size(N) == size(A)
        @test !isreal(M)
        @test ndims(M) == 2
        @test_throws ErrorException size(M, 3)
        @test length(M) == length(A)
        @test occursin(Base.dims2string(size(M)) * " LinearMaps.WrappedMap", summary(M))
    end

    @testset "dimension checking" begin
        w = vec(u)
        @test_throws DimensionMismatch M * similar(v, length(v) + 1)
        @test_throws DimensionMismatch mul!(similar(w, length(w) + 1), M, v)
        @test_throws DimensionMismatch similar(w, length(w) + 1)' * M
        @test_throws DimensionMismatch mul!(copy(v)', similar(w, length(w) + 1)', M)
        @test_throws DimensionMismatch mul!(similar(W, size(W).+(0,1)), M, V)
        @test_throws DimensionMismatch mul!(copy(W), M, similar(V, size(V).+(0,1)))
    end
end

# new type
struct SimpleFunctionMap <: LinearMap{Float64}
    f::Function
    N::Int
end
struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
    f::Function
    N::Int
end
Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::AbstractVector) = A.f(v)
LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::AbstractVector) = copyto!(y, *(A, x))

@testset "new LinearMap type" begin
    F = SimpleFunctionMap(cumsum, 10)
    @test parent(F) === F
    FC = SimpleComplexFunctionMap(cumsum, 10)
    @test @inferred ndims(F) == 2
    @test @inferred size(F, 1) == 10
    @test @inferred length(F) == 100
    @test @inferred !issymmetric(F)
    @test @inferred !ishermitian(F)
    @test @inferred !ishermitian(FC)
    @test @inferred !isposdef(F)
    @test occursin("10×10 SimpleFunctionMap{$(eltype(F))}", sprint((t, s) -> show(t, "text/plain", s), F))
    @test occursin("10×10 SimpleComplexFunctionMap{$(eltype(FC))}", sprint((t, s) -> show(t, "text/plain", s), FC))
    α = rand(ComplexF64); β = rand(ComplexF64)
    v = rand(ComplexF64, 10); V = rand(ComplexF64, 10, 3)
    w = rand(ComplexF64, 10); W = rand(ComplexF64, 10, 3)
    @test mul!(w, F, v) === w == F * v
    @test_throws ErrorException F' * v
    @test_throws ErrorException transpose(F) * v
    @test_throws ErrorException mul!(w, adjoint(FC), v)
    @test_throws ErrorException mul!(w, transpose(F), v)
    FM = convert(AbstractMatrix, F)
    L = LowerTriangular(ones(10, 10))
    @test FM == L
    @test F * v ≈ L * v
    # generic 5-arg mul! and matrix-mul!
    @test mul!(copy(w), F, v, α, β) ≈ L*v*α + w*β
    @test mul!(copy(w), F, v, 0, β) ≈ w*β
    @test mul!(copy(W), F, V) ≈ L*V
    @test mul!(copy(W), F, V, α, β) ≈ L*V*α + W*β
    @test mul!(copy(W), F, V, 0, β) ≈ W*β

    Fs = sparse(F)
    @test SparseMatrixCSC(F) == Fs == L
    @test Fs isa SparseMatrixCSC
end

struct MyFillMap{T} <: LinearMaps.LinearMap{T}
    λ::T
    size::Dims{2}
    function MyFillMap(λ::T, dims::Dims{2}) where {T}
        all(d -> d >= 0, dims) || throw(ArgumentError("dims of MyFillMap must be non-negative"))
        promote_type(T, typeof(λ)) == T || throw(InexactError())
        return new{T}(λ, dims)
    end
end
Base.size(A::MyFillMap) = A.size
function LinearAlgebra.mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector)
    LinearMaps.check_dim_mul(y, A, x)
    return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
end
function LinearAlgebra.mul!(y::AbstractVecOrMat, transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap}, x::AbstractVector)
    LinearMaps.check_dim_mul(y, transA, x)
    λ = transA.lmap.λ
    return fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x))
end

@testset "transpose of new LinearMap type" begin
    A = MyFillMap(5.0, (3, 3))
    x = ones(3)
    @test A * x == fill(15.0, 3)
    @test mul!(zeros(3), A, x) == mul!(zeros(3), A, x, 1, 0) == fill(15.0, 3)
    @test mul!(ones(3), A, x, 2, 2) == fill(32.0, 3)
    @test mul!(ones(3,3), A, reshape(collect(1:9), 3, 3), 2, 2)[:,1] == fill(62.0, 3)
    @test A'x == mul!(zeros(3), A', x)
    for α in (true, false, 0, 1, randn()), β in (true, false, 0, 1, randn())
        @test mul!(ones(3), A', x, α, β) == fill(β, 3) + fill(15α, 3)
        @test mul!(ones(3, 2), A', [x x], α, β) == fill(β, 3, 2) + fill(15α, 3, 2)
    end
end


text:
include("rob_lib.jl")

function mark_pyramid(r::Robot)
    num_verticale = digit_movements!(r, Sud)
    num_horizont = digit_movements!(r, West)

    side = Ost
    all_steps = digit_movements!(r, Ost)
    max_steps = all_steps
    digit_movements!(r, West)
# while we arent in the corner we cant reach thr satisfaction
    while !(isborder(r, Nord) && (isborder(r, Ost) || isborder(r, West)))
        now_steps = digit_movements!(r, Ost, max_steps)
        digit_movements!(r, West)
        putmarker!(r)
        putmarkers!(r, side, now_steps)
        Back!(r, inverse(side), now_steps)
        move!(r, Nord)
        max_steps -= 1
    end

    now_steps = digit_movements!(r, Ost, max_steps)
    digit_movements!(r, West)
    putmarker!(r)
    putmarkers!(r, side, now_steps)

    digit_movements!(r, Sud)
    digit_movements!(r, West)

    Back!(r, Nord, num_verticale)
    Back!(r, Ost, num_horizont)
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
import NNlibCUDA: batchnorm, ∇batchnorm

function (BN::Flux.BatchNorm)(x::Union{CuArray{T,2},CuArray{T,4},CuArray{T,5}},
                              cache=nothing) where T<:Union{Float32, Float64}

  @assert BN.affine "BatchNorm: only affine=true supported on gpu"
  @assert BN.track_stats "BatchNorm: only track_stats=true supported on gpu"
  @assert length(BN.β) == size(x, ndims(x)-1) "BatchNorm: input has wronng number of channels"
  return BN.λ.(batchnorm(BN.γ, BN.β, x, BN.μ, BN.σ², BN.momentum;
                  cache=cache, alpha=1, beta=0, eps=BN.ϵ,
                  training=Flux._isactive(BN)))
end

function ChainRulesCore.rrule(::typeof(batchnorm), g, b, x, running_mean, running_var, momentum; kw...)
  y = batchnorm(g, b, x, running_mean, running_var, momentum; kw...)
  function batchnorm_pullback(Δ)
    grad = ∇batchnorm(g, b, x, Δ, running_mean, running_var, momentum; kw...)
    (NoTangent(), grad..., NoTangent(), NoTangent(), NoTangent())
  end
  y, batchnorm_pullback
end


text:
"""
    quaddata(operator, test_refspace, trial_refspace, test_elements, trial_elements)

Returns an object cashing data required for the computation of boundary element
interactions. It is up to the client programmer to decide what (if any) data is
cached. For double numberical quadrature, storing the integration points for
example can significantly speed up matrix assembly.

- `operator` is an integration kernel.
- `test_refspace` and `trial_refspace` are reference space objects. `quadata`
is typically overloaded on the type of these local spaces of shape functions.
(See the implementation in `maxwell.jl` for an example).
- `test_elements` and `trial_elements` are iterable collections of the geometric
elements on which the finite element space are defined. These are provided to
allow computation of the actual integrations points - as opposed to only their
coordinates.
"""
function quaddata end

"""
    quadrule(operator, test_refspace, trial_refspace, test_index, test_chart, trial_index, trial_chart, quad_data)

Based on the operator kernel and the test and trial elements, this function builds
an object whose type and data fields specify the quadrature rule that needs to be
used to accurately compute the interaction integrals. The `quad_data` object created
by `quaddata` is passed to allow reuse of any precomputed data such as quadrature
points and weights, geometric quantities, etc.

The type of the returned quadrature rule will help in deciding which method of
`momintegrals` to dispatch to.
"""
function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd)
    # defines coincidence of points
    dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))

    # decides on whether to use singularity extraction
    xtol = 0.2
    k = norm(op.gamma)

    hits = 0
    xmin = xtol
    for t in τ.vertices
      for s in σ.vertices
        d = norm(t-s)
        xmin = min(xmin, k*d)
        if d < dtol
          hits +=1
          break
        end
      end
    end

  hits == 3   && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
  hits == 2   && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
  hits == 1   && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
  xmin < xtol && return DoubleQuadRule(
      qd.tpoints[2,i],
      qd.bpoints[2,j],)
  return DoubleQuadRule(
    qd.tpoints[1,i],
    qd.bpoints[1,j],)
end


text:
const Passport = Dict{String,String}

function readinput(filepath)

    function parsepassport(str)
        p = Passport()
        ss = split(str, r"\s")
        for s in ss
            m = match(r"([a-z]{3}):(.+)$", s)
            p[m[1]] = m[2]
        end
        p
    end

    strs = open(filepath) do file
        split(strip(read(file, String)), "\n\n")
    end
    map(parsepassport, strs)
end

# First Half

allrqdflds(p) = length(keys(p)) == (haskey(p, "cid") ? 8 : 7) # Treat cid as optional.

passports = readinput("inputs/day04.txt")
count(allrqdflds, passports)
# Count valid passports - those that have all required fields.

# Second Half

function allvldflds(p)
    # True if all required fields are both present and valid

    if !allrqdflds(p)

        return false
    end

    if !(parse(Int, p["byr"]) ∈ 1920:2002)
        # byr (Birth Year) - four digits; at least 1920 and at most 2002

        return false
    end
    if !(parse(Int, p["iyr"]) ∈ 2010:2020)
        # iyr (Issue Year) - four digits; at least 2010 and at most 2020

        return false
    end
    if !(parse(Int, p["eyr"]) ∈ 2020:2030)
        # eyr (Expiration Year) - four digits; at least 2020 and at most 2030

        return false
    end
    m = match(r"^(\d+)(cm|in)$", p["hgt"])
    if isnothing(m)
        # hgt (Height) - a number followed by either cm or in

        return false
    elseif m[2] == "cm"
        if !(parse(Int, m[1]) ∈ 150:193)
            # If cm, the number must be at least 150 and at most 193

            return false
        end
    else
        if !(parse(Int, m[1]) ∈ 59:76)
            # If in, the number must be at least 59 and at most 76

            return false
        end
    end
    if isnothing(match(r"^#[0-9|a-f]{6}$", p["hcl"]))
        # hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f

        return false
    end
    if isnothing(match(r"^(amb|blu|brn|gry|grn|hzl|oth){1}$", p["ecl"]))
        # ecl (Eye Color) - exactly one of: amb blu brn gry grn hzl oth

        return false
    end
    if isnothing(match(r"^\d{9}$", p["pid"]))
        # pid (Passport ID) - a nine-digit number, including leading zeroes

        return false
    end
    true
end

passports = readinput("inputs/day04.txt")
count(allvldflds, passports)
# Count the number of valid passports - those that have
# all required fields and valid values


text:
function Test_AdaptiveSteppers(t0, tf)
    function func(t, y)
        [exp(t)]
    end

    y0, exact = [exp(t0)], exp(tf)
    methods = DopriStep, CarpStep, RKFStep

    for a=methods
        ts, ys = AdaptiveStepper(func, a, t0, copy(y0), tf)
        println("Actual error of ", a, ": ", abs(exact - ys[end][1]))
    end
end

#Test_AdaptiveSteppers(20, 21)

#### Systems ####

function SimpleSystem(t, y)
    [cos(t)]
end

function HarmonicOscillator(t, y)
    y1, ydot = y
    [ydot, 100*cos(20*t) - ydot/4.0 - y1]
end

function FourthOrderSystem(t, y)
    y1, ydot, yddot, ydddot = y
    [ydot, yddot, ydddot, exp(t)*cos(100*t)-7*ydddot-17*yddot-17*ydot-6*y1]
end

#t0, tf = 0, 15
#y0 = [0, 0, 0, 0]

#ts, ys = Stepper(FourthOrderSystem, BackwardEulerStep, t0, y0, tf, 0.01)
#ts, ys = MultiStepper(FourthOrderSystem, AdamsBash4, t0, y0, tf, 0.001)
#plot(ts, [a[1] for a=ys], "o-")

max_stars_repo_name: jiosue/QuadratureCandJulia

max_issues_repo_name: jiosue/QuadratureCandJulia

max_forks_repo_name: jiosue/QuadratureCandJulia


text:
"""
Convenience code for plotting CCFs

Author: Eric Ford
Date:   October 2020
"""

"""   `make_plot_ccf_vs_time(ccfs; v_grid, t_idx, output_path, save_fig )`
Overplot CCFs (2d array) at multiple times.
Inputs:
- `ccfs`:  2d array of CCFs at different times
Optional Inputs:
- `v_grid`: velocities that CCF was evaluated at
- `t_idx`:  restrict plotting to these time indicies
- `save_fig`: save figures to file if true
- `output_path`: prefix to filename
Returns:
- plot object
"""
function make_plot_ccf_vs_time(ccfs::AbstractArray{T2,2}
                ; v_grid::AbstractRange = 1:size(ccfs,1),
                t_idx = 1:size(ccfs,2),
                output_path::String = "", save_fig::Bool = false ) where
                { T1<:Real, T2<:Real, T3<:Real, T4<:Real  }

    @assert length(v_grid) == size(ccfs,1)
    @assert length(t_idx) >= 1

    plt = plot()
    plot!(plt,v_grid,ccfs[:,t_idx]./maximum(ccfs[:,t_idx],dims=1),label=:none)
    xlabel!("v (m/s)")
    ylabel!("CCF")
    if save_fig   savefig(plt,output_path * "_ccf_sum.png")   end
    return plt
end

"""   `make_heatmap_ccf_vs_time(ccfs; v_grid, t_idx, output_path, save_fig )`
Plots heatmap of CCFs at multiple times minus the mean CCF (averaged over time).
Inputs:
- `ccfs`:  2d array of CCFs at different times
Optional Inputs:
- `v_grid`: velocities that CCF was evaluated at
- `save_fig`: save figures to file if true
- `output_path`: prefix to filename
Returns:
- plot object
"""
function make_heatmap_ccf_vs_time(ccfs::AbstractArray{T2,2}
                ; v_grid::AbstractRange = 1:size(ccfs,1),
                times::AbstractVector{T4} = collect(1:length(rvs)),
                output_path::String = "", save_fig::Bool = false ) where
                { T1<:Real, T2<:Real, T3<:Real, T4<:Real  }

    @assert length(v_grid) == size(ccfs,1)

    ccfs_minus_mean = ccfs./maximum(ccfs,dims=1) .- mean( ccfs./ maximum(ccfs,dims=1),dims=2)
    zvals = ccfs_minus_mean
    colorscale = cgrad(:balance)
    plt = plot()
    heatmap!(plt, v_grid, 1:size(ccfs,2),zvals',c=colorscale, clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
    add_time_gap_lines(plt,times)
    title!("CCF(v,t)-<CCF> vs Time")
    xlabel!("v (m/s)")
    ylabel!("Observation #")
    if save_fig
        savefig(joinpath(output_path,"_ccf_sum_vs_time_heatmap.png"))
    end
    return plt
end

"""   `make_heatmap_ccf_vs_order(ccfs; v_grid, t_idx, output_path, save_fig )`
Plots heatmap of CCFs for multiple overs minus the mean CCF (averaged over time).
Inputs:
- `ccfs`:  2d array of CCFs at different times
Optional Inputs:
- `v_grid`: velocities that CCF was evaluated at
- `save_fig`: save figures to file if true
- `output_path`: prefix to filename
Returns:
- plot object
"""
function make_heatmap_ccf_vs_order(order_ccfs::AbstractArray{T2,3}
                ; v_grid::AbstractRange = 1:size(order_ccfs,1),
                #= times::AbstractVector{T4} = collect(1:length(rvs)), =#
                order_labels::AbstractVector{T5} = 1:size(order_ccfs,2),
                output_path::String = "", save_fig::Bool = false ) where
                { T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real  }

    @assert length(v_grid) == size(order_ccfs,1)

    obs = 1:size(order_ccfs,3)
    #order_labels = map(c->order_list_timeseries.chunk_list[1].data[c].λ.indices[2], 1:size(order_ccfs,2))
    orders_to_plot = findall(map(c->sum(order_ccfs[:,c,obs])>0, 1:size(order_ccfs,2)))
    zvals =  reshape(sum(order_ccfs[:,orders_to_plot,obs],dims=3)./maximum(sum(order_ccfs[:,orders_to_plot,obs],dims=3),dims=1),size(order_ccfs,1),size(order_ccfs[:,orders_to_plot,obs],2)) .-
             reshape(sum(order_ccfs[:,orders_to_plot,obs],dims=(2,3))./maximum(sum(order_ccfs[:,orders_to_plot,obs],dims=(2,3))),size(order_ccfs,1) )
    colorscale = cgrad(:balance)
    plt = plot()
    heatmap!(plt,v_grid,order_labels[orders_to_plot], zvals', clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
    #println("size(v_grid) = ", size(v_grid_order_ccfs), "  size(order_labels[plt]) = ", size(order_labels[orders_to_plot]), "  size(zvals') = ", size(zvals') )
    xlabel!("v (m/s)")
    ylabel!("Order ID")
    title!("CCF-<CCF> for order_idx = " * string(order_labels[minimum(orders_to_plot)]) * " - " * string(order_labels[maximum(orders_to_plot)]))
    if save_fig
            if save_plot(pipeline_plan,:ccf_orders)   savefig(plt,output_path * " _ccf_orders.png")   end
    end

    return plt
end


text:
for (op, i) in ((:two,:2), (:three,:3), (:four, :4))
  @eval begin
    function ($op){P}(::Type{ArbFloat{P}})
        z = ArbFloat{P}()
        ccall(@libarb(arb_set_si), Void, (Ptr{ArbFloat}, Int), &z, $i)
        return z
    end
    ($op)(::Type{ArbFloat}) = ($op)(ArbFloat{precision(ArbFloat)})
  end
end

weakcopy{P}(x::ArbFloat{P}) = WeakRef(x)

for fn in (:copy, :deepcopy)
  @eval begin
    function ($fn){P}(x::ArbFloat{P})
        z = (ArbFloat{P})()
        ccall(@libarb(arb_set), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
        return z
    end
    function ($fn){T<:ArbFloat}(x::T)
        z = T()
        ccall(@libarb(arb_set), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
        return z
    end
  end
end

function copyradius{T<:ArbFloat}(target::T, source::T)
    z = deepcopy(target)
    z.radius_exponentOf2 = source.radius_exponentOf2
    z.radius_significand = source.radius_significand
    return z
end

function copymidpoint{T<:ArbFloat}(target::T, source::T)
    z = deepcopy(target)
    z.exponentOf2  = source.exponentOf2
    z.nwords_sign  = source.nwords_sign
    z.significand1 = source.significand1
    z.significand2 = source.significand2
    return z
end

function midpoint_radius{T<:ArbFloat}(x::T)
    return midpoint(x), radius(x)
end

function midpoint_radius{T<:ArbFloat}(midpt::T, radius::T)
    z = midpoint(midpt)
    ccall(@libarb(arb_add_error), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &radius)
    return z
end

function midpoint_radius{T<:ArbFloat}(midpoint::T, radius::Float64)
    rad = convert(T, radius)
    return midpoint_radius(midpoint, rad)
end

function bounds{T<:ArbFloat}(lower::T, upper::T)
    lowerlo, lowerhi = bounds(lower)
    upperlo, upperhi = bounds(upper)
    lo = lowerlo <= upperlo ? lowerlo : upperlo
    hi = lowerhi >= upperhi ? lowerhi : upperhi
    # rad = (hi - lo) * 0.5
    mid = hi*0.5 + lo*0.5
    rad = hi - mid
    r = Float32(rad)
    dr = 0.5 * eps(r)
    z = midpoint_radius(mid, r)
    tstlo, tsthi = bounds(z)
    while (tstlo > lo) || (tsthi < hi)
        z = midpoint_radius(mid, r+dr)
        tstlo, tsthi = bounds(z)
    end
    return z
end

"""
Rounds x to a number of bits equal to the accuracy of x (as indicated by its radius), plus a few guard bits.
The resulting ball is guaranteed to contain x, but is more economical if x has less than full accuracy.
(from arb_trim documentation)
"""
function trim{T<:ArbFloat}(x::T)
    z = T()
    ccall(@libarb(arb_trim), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
    return z
end

"""
Rounds x to a clean estimate of x as a point value.
"""
function tidy{T<:ArbFloat}(x::T)
    s = smartarbstring(x)
    return (T)(s)
end

function decompose{T<:ArbFloat}(x::T)
    # decompose x as num * 2^pow / den
    # num, pow, den = decompose(x)
    P = precision(T)
    bfprec=precision(BigFloat)
    setprecision(BigFloat,P)
    bf = convert(BigFloat, x)
    n,p,d = decompose(bf)
    setprecision(BigFloat,bfprec)
    return n,p,d
end


text:
using Glob
using JSON

any_errors = false

Threads.@threads for cfg_path in glob("languages/*/exercises/concept/*/.meta/config.json")
    split_path = splitpath(cfg_path)
    exercise_path = joinpath(split_path[1:end - 2]...)
    exercise = joinpath(split_path[2], split_path[5]) # e.g. julia/lasagna

    local cfg
    try
        cfg = JSON.parsefile(cfg_path, use_mmap=false) # use_mmap=false makes it work on Windows
    catch e
        @warn "$cfg_path is not a valid JSON file" e
        continue
    end

    if !haskey(cfg, "editor")
        @warn "$exercise has no editor settings"# cfg
        continue
    end

    missing_files = [f for f in union(cfg["editor"]["test_files"], cfg["editor"]["solution_files"]) if !isfile(joinpath(exercise_path, f))]
    if length(missing_files) > 0
        @error "$exercise contains references to files that don't exist" cfg["editor"] missing_files
        global any_errors = true
    end
end

any_errors && exit(1)


text:
@testset "682.baseball-game.jl" begin
    @test calculate_points(["5", "2", "C", "D", "+"]) == 30
    @test calculate_points(["5", "-2", "4", "C", "D", "9", "+", "+"]) == 27
    @test calculate_points(["1"]) == 1
    @test calculate_points(["4", "2", "C", "C"]) == 0
end


text:
using Distributed

mutable struct WorkMaster <: Master
    tasks
    results
    slaves
end

nparalel(master::WorkMaster) = length(master.slaves)

"""
Releases one worker (who is the next one of taking a new task) from the duty calculating function values given him by the master. Returns pid.
"""
function releaseslave!(master::WorkMaster)
    @assert length(master.slaves)>0

    put!(master.tasks,nothing)

    while true
        for i in 1:length(master.slaves)
            ### I need to find first one which is free
            if isready(master.slaves[i])==true
                pid = master.slaves[i].where
                deleteat!(master.slaves,i)
                return pid
            end
        end
    end
end

"""
Releases all workers from the Master's duties.
"""
function releaseall!(master::WorkMaster)
    for s in master.slaves
        put!(master.tasks,nothing)
    end

    for s in master.slaves
        wait(s)
    end

    pids = [s.where for s in master.slaves]
    master.slaves = []

    return pids
end

"""
Gives the slave a duty to follow orders of his new Master
"""
function captureslave!(pid,f::Function,master::WorkMaster)
    tasks, results = master.tasks, master.results
    wp = @spawnat pid begin
        while true
            x = take!(tasks)
            if x==nothing
                break
            else
                y = f(x)
                put!(results,(x,y))
            end
        end
    end
    push!(master.slaves,wp)
end

WorkMaster(tasks,results) = WorkMaster(tasks,results,[])

function WorkMaster(f::Function,wpool::AbstractWorkerPool)
    tasks = RemoteChannel(()->Channel{Any}(10))
    results = RemoteChannel(()->Channel{Tuple{Any,Any}}(10))

    master = WorkMaster(tasks,results)
    for p in wpool.workers
        captureslave!(p,f,master)
    end

    return master
end

WorkMaster(f::Function) = WorkMaster(f,WorkerPool(nprocs()==1 ? [1] : workers()))


text:
module IfPartialBase

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export if_pb

function if_pb(m)

    #if (flbz_front_left_bot_z[i] == flbz_front_left_bot_z[j] + hob_height_of_box[j] &&
    #    flbz_front_left_bot_z[i] - hob_height_of_box[j] >= 0 &&
    #    flbx_front_left_bot_x[j] > flbx_front_left_bot_x[i] - box_i_x_axis_cover &&
    #    flbx_front_left_bot_x[j] < flbx_front_left_bot_x[i] + box_i_x_axis_cover &&
    #    flby_front_left_bot_y[j] > flby_front_left_bot_y[i] - box_i_y_axis_cover &&
    #    flby_front_left_bot_y[j] < flby_front_left_bot_y[i] + box_i_y_axis_cover)

    # IF 1
    @variable(m, if11_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @variable(m, if12_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @variable(m, if1_pbs[i=1:length(boxes), j=1:length(boxes); i!=j] >= 0, Bin)

    #Hz1 >= x-b+e
    #H(1-z1) >= b-x-e
    @constraint(m, c1_if11_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * m[:if11_pbs][i, j] >=
        m[:flbz_front_left_bot_z][i] - (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) + 0.001
    )
    @constraint(m, c2_if11_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * (1 - m[:if11_pbs][i, j]) >=
        (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - m[:flbz_front_left_bot_z][i] - 0.001
    )

    #Hz2 >= b-x+e
    #H(1-z2) >= x-b-e
    @constraint(m, c1_if12_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * m[:if12_pbs][i, j] >=
        (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - m[:flbz_front_left_bot_z][i] + 0.001
    )
    @constraint(m, c2_if12_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * (1 - m[:if12_pbs][i, j]) >=
        m[:flbz_front_left_bot_z][i] - (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - 0.001
    )

    #If if11, if12 = 1, then box i and j have z-bot/top value within interval -e <= z <= +e, hence, expression = 1
    @constraint(m, c_if1_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if11_pbs][i, j] + m[:if12_pbs][i, j] >=
        2 * m[:if1_pbs][i, j]
    )

    #Xj < Xi + Xcoveri
    @variable(m, if2_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if2_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * (1 - m[:if2_pbs][i, j]) >=
        m[:flbx_front_left_bot_x][j] + 0.001 - (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i])
    )

    @constraint(m, c2_if2_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * m[:if2_pbs][i, j] >=
        (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i]) - m[:flbx_front_left_bot_x][j] - 0.001
    )

    # Xj + Xcoverj > Xi
    @variable(m, if3_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if3_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * (1 - m[:if3_pbs][i, j]) >=
        m[:flbx_front_left_bot_x][i] + 0.001 - (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j])
    )

    @constraint(m, c2_if3_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * m[:if3_pbs][i, j] >=
        (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j]) - m[:flbx_front_left_bot_x][i] - 0.001
    )

    # Yj < Yi + Ycoveri
    @variable(m, if4_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if4_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * (1 - m[:if4_pbs][i, j]) >=
        m[:flby_front_left_bot_y][j] + 0.001 - (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i])
    )

    @constraint(m, c2_if4_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * m[:if4_pbs][i, j] >=
        (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i]) - m[:flby_front_left_bot_y][j] - 0.001
    )

    # Yj + Ycoverj > Yi
    @variable(m, if5_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if5_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * (1 - m[:if5_pbs][i, j]) >=
        m[:flby_front_left_bot_y][i] + 0.001 - (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j])
    )

    @constraint(m, c2_if5_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * m[:if5_pbs][i, j] >=
        (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j]) - m[:flby_front_left_bot_y][i] - 0.001
    )

    @variable(m, if_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if1_pbs][i, j] +
        m[:if2_pbs][i, j] +
        m[:if3_pbs][i, j] +
        m[:if4_pbs][i, j] +
        m[:if5_pbs][i, j] >=
        5 * m[:if_on_top][i, j]
    )

    return m
end

end


text:
@i function imul!(C::StridedVecOrMat, A::AbstractSparseMatrix{T}, B::DenseInputVecOrMat, α::Number, β::Number) where T
    @safe size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    @safe size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    @safe size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    if (β != 1, ~)
        @safe error("only β = 1 is supported, got β = $(β).")
    end
    # Here, we close the reversibility check inside the loop to increase performance
    @invcheckoff for k = 1:size(C, 2)
        @inbounds for col = 1:size(A, 2)
            @routine begin
                αxj ← zero(eltype(B))
                αxj += B[col,k] * α
            end
            for j = nzrange(A, col)
                C[A.rowval[j], k] += A.nzval[j] * αxj
            end
            ~@routine
        end
    end
end

@i function imul!(C::StridedVecOrMat, xA::Adjoint{T, <:AbstractSparseMatrix}, B::DenseInputVecOrMat, α::Number, β::Number) where T
    @safe size(xA.parent, 2) == size(C, 1) || throw(DimensionMismatch())
    @safe size(xA.parent, 1) == size(B, 1) || throw(DimensionMismatch())
    @safe size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    if (β != 1, ~)
        @safe error("only β = 1 is supported, got β = $(β).")
    end
    @invcheckoff for k in 1:size(C, 2)
        @inbounds for col in 1:size(xA.parent, 2)
                for j in nzrange(xA.parent, col)
                    @routine begin
                        @zeros T anc1 anc2
                        anc1 += (xA.parent.nzval[j])'* α
                        anc2 += B[xA.parent.rowval[j], k]
                    end
                C[col,k] += anc1 * anc2
                ~@routine
            end
        end
    end
end

@i function imul!(C::StridedVecOrMat, X::DenseMatrixUnion, A::AbstractSparseMatrix{T}, α::Number, β::Number) where T
    @safe size(X, 2) == size(A, 1) || throw(DimensionMismatch())
    @safe size(X, 1) == size(C, 1) || throw(DimensionMismatch())
    @safe size(A, 2) == size(C, 2) || throw(DimensionMismatch())
    if (β != 1, ~)
        @safe error("only β = 1 is supported, got β = $(β).")
    end
    @invcheckoff for col in 1:size(A, 2)
        @inbounds for k in nzrange(A, col)
            @simd for multivec_row in 1:size(X,1)
                @routine begin
                    anc1 ← zero(T)
                    anc2 ← zero(promote_eltype(X, A))
                    anc1 += X[multivec_row, A.rowval[k]]
                    anc2 += anc1 * A.nzval[k]
                end
                C[multivec_row, col] += anc2 * α
                ~@routine
            end
        end
    end
end

@i function imul!(C::StridedVecOrMat, X::Adjoint{T1, <:DenseMatrixUnion}, A::AbstractSparseMatrix{T2}, α::Number, β::Number) where {T1, T2}
    @safe size(X.parent, 1) == size(A, 1) || throw(DimensionMismatch())
    @safe size(X.parent, 2) == size(C, 1) || throw(DimensionMismatch())
    @safe size(A, 2) == size(C, 2) || throw(DimensionMismatch())
    if (β != 1, ~)
        @safe error("only β = 1 is supported, got β = $(β).")
    end
    @invcheckoff for multivec_row in 1:size(X.parent, 2)
        for col in 1:size(A, 2)
            @inbounds for k in nzrange(A, col)
                @routine begin
                    anc1 ← zero(T1)
                    anc2 ← zero(promote_type(T1, T2))
                    anc1 += X.parent[A.rowval[k], multivec_row]'
                    anc2 += anc1 * A.nzval[k]
                end
                C[multivec_row, col] += anc2 * α
                ~@routine
            end
        end
    end
end

@i function imul!(C::StridedVecOrMat, X::DenseMatrixUnion, xA::Adjoint{T, <:AbstractSparseMatrix}, α::Number, β::Number) where T
    @safe size(X, 2) == size(xA.parent, 2) || throw(DimensionMismatch())
    @safe size(X, 1) == size(C, 1) || throw(DimensionMismatch())
    @safe size(xA.parent, 1) == size(C, 2) || throw(DimensionMismatch())
    if (β != 1, ~)
        @safe error("only β = 1 is supported, got β = $(β).")
    end
    @invcheckoff for col in 1:size(xA.parent, 2)
        @inbounds for k in nzrange(xA.parent, col)
                        @routine begin
                            anc1 ← zero(T)
                            anc1 += (xA.parent.nzval[k])'*α
                        end
            @simd for multivec_row in 1:size(X,1)
                C[multivec_row, xA.parent.rowval[k]] += X[multivec_row, col] * anc1
            end
            ~@routine
        end
    end
end

# https://github.com/JuliaLang/julia/blob/e2aeefb60eaffbd2807089155789e19585f4e6fe/stdlib/SparseArrays/src/linalg.jl#L301
@i function idot(r::T, A::SparseMatrixCSC{T},B::SparseMatrixCSC{T}) where {T}
    @routine begin
        (m, n) ← size(A)
        branch_keeper ← zeros(Bool, 2*m)
    end
    @safe size(B) == (m,n) || throw(DimensionMismatch("matrices must have the same dimensions"))
    @invcheckoff @inbounds for j = 1:n
        @routine begin
            ia1 ← A.colptr[j]
            ib1 ← B.colptr[j]
            ia2 ← A.colptr[j+1]
            ib2 ← B.colptr[j+1]
            ia ← ia1
            ib ← ib1
        end
        if ib2>ib1 && ia2>ia1
            for i=1:ia2-ia1+ib2-ib1-1
                ra ← A.rowval[ia]
                rb ← B.rowval[ib]
                if ra == rb
                    r += A.nzval[ia]' * B.nzval[ib]
                end
                ## b move -> true, a move -> false
                branch_keeper[i] ⊻= @const ia >= ia2-1 || (ib < ib2-1 && ra > rb)
                ra → A.rowval[ia]
                rb → B.rowval[ib]
                if (branch_keeper[i], ~)
                    INC(ib)
                else
                    INC(ia)
                end
            end
            ~for i=1:ia2-ia1+ib2-ib1-1
                ## b move -> true, a move -> false
                branch_keeper[i] ⊻= @const ia >= ia2-1 || (ib < ib2-1 && A.rowval[ia] > B.rowval[ib])
                if (branch_keeper[i], ~)
                    INC(ib)
                else
                    INC(ia)
                end
            end
        end
        ~@routine
    end
    ~@routine
end

@i function idot(r, x::AbstractVector, A::AbstractSparseMatrix{T1}, y::AbstractVector{T2}) where {T1, T2}
    @safe length(x) == size(A, 1) || throw(DimensionMismatch())
    @safe length(y) == size(A, 2) || throw(DimensionMismatch())

    @invcheckoff @inbounds for col in 1:size(A, 2)
        for k in nzrange(A, col)
            @routine begin
                anc ← zero(promote_type(T1, T2))
                anc += A.nzval[k] * y[col]
            end
            r += (x[A.rowval[k]])' * anc
            ~@routine
        end
    end
end

@i function idot(r, x::SparseVector, A::AbstractSparseMatrix{T1}, y::SparseVector{T2}) where {T1, T2}
    @safe length(x) == size(A, 1) || throw(DimensionMismatch())
    @safe length(y) == size(A, 2) || throw(DimensionMismatch())

    @invcheckoff for m in 1:nnz(y)
        for j in  1:nnz(x)
            for k in nzrange(A, y.nzind[m])
                if A.rowval[k] == x.nzind[j]
                    @routine begin
                        anc ← zero(promote_type(T1, T2))
                        anc += A.nzval[k] * y.nzval[m]
                    end
                    r += x.nzval[j]' * anc
                    ~@routine
                end
            end
        end
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -8.966433788908274 <= q <= 11.565431391005125)

                     add_NL_constraint(m, :(log(1 + exp(0.8271637266959466 + -0.29461988642396975*$(x[1]) + 0.9225334019614664*$(x[2]) + -0.12844201057842541*$(x[3]) + -0.9795769523241842*$(x[4]) + 0.7579371331675153*$(x[5]))) + log(1 + exp(0.9202467558270198 + 0.816513326282978*$(x[1]) + 0.1577763582178746*$(x[2]) + 0.37608788225609757*$(x[3]) + 0.8366994089368287*$(x[4]) + 0.28442920912454284*$(x[5]))) + log(1 + exp(-0.23055889086042392 + -0.12029185584830726*$(x[1]) + 0.7385896748366556*$(x[2]) + 0.04575327278060515*$(x[3]) + 0.05400397335314544*$(x[4]) + -0.13287838316667244*$(x[5]))) + log(1 + exp(0.42724629038540085 + 0.286409417560213*$(x[1]) + -0.3599014579136197*$(x[2]) + -0.15981498421265528*$(x[3]) + -0.05370590457282587*$(x[4]) + -0.6772611532220885*$(x[5]))) + log(1 + exp(-0.6445990809995172 + -0.48590839796846375*$(x[1]) + 0.5462761102799756*$(x[2]) + -0.2471758620035076*$(x[3]) + 0.5819097771517305*$(x[4]) + 0.22143679581234954*$(x[5]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Oceananigans.BoundaryConditions: BoundaryFunction

function test_z_boundary_condition_simple(arch, FT, fldname, bctype, bc, Nx, Ny)
    Nz = 16
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), extent=(0.1, 0.2, 0.3))

    bc = BoundaryCondition(bctype, bc)
    field_bcs = TracerBoundaryConditions(grid, top=bc)
    model_bcs = NamedTuple{(fldname,)}((field_bcs,))

    model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT,
                                boundary_conditions=model_bcs)

    time_step!(model, 1e-16, euler=true)

    return model isa IncompressibleModel
end

function test_z_boundary_condition_array(arch, FT, fldname)
    Nx = Ny = Nz = 16

    bcarray = rand(FT, Nx, Ny)

    if arch == GPU()
        bcarray = CuArray(bcarray)
    end

    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), extent=(0.1, 0.2, 0.3))

    value_bc = BoundaryCondition(Value, bcarray)
    field_bcs = TracerBoundaryConditions(grid, top=value_bc)
    model_bcs = NamedTuple{(fldname,)}((field_bcs,))

    model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT,
                                boundary_conditions=model_bcs)

    time_step!(model, 1e-16, euler=true)

    field = get_model_field(fldname, model)
    bcs = field.boundary_conditions
    return bcs.top[1, 2] == bcarray[1, 2]
end

function test_flux_budget(arch, FT, fldname)
    N, κ, Lz = 16, 1, 0.7
    grid = RegularCartesianGrid(FT, size=(N, N, N), extent=(1, 1, Lz))

    bottom_flux = FT(0.3)
    flux_bc = BoundaryCondition(Flux, bottom_flux)

    if fldname == :u
        field_bcs = UVelocityBoundaryConditions(grid, bottom=flux_bc)
    elseif fldname == :v
        field_bcs = VVelocityBoundaryConditions(grid, bottom=flux_bc)
    else
        field_bcs = TracerBoundaryConditions(grid, bottom=flux_bc)
    end

    model_bcs = NamedTuple{(fldname,)}((field_bcs,))

    closure = IsotropicDiffusivity(FT, ν=κ, κ=κ)
    model = IncompressibleModel(grid=grid, closure=closure, architecture=arch, tracers=(:T, :S),
                                float_type=FT, buoyancy=nothing, boundary_conditions=model_bcs)

    field = get_model_field(fldname, model)
    @. field.data = 0

    τκ = Lz^2 / κ   # Diffusion time-scale
    Δt = 1e-6 * τκ  # Time step much less than diffusion time-scale
    Nt = 10         # Number of time steps

    for n in 1:Nt
        time_step!(model, Δt, euler= n==1)
    end

    # budget: Lz*∂<ϕ>/∂t = -Δflux = -top_flux/Lz (left) + bottom_flux/Lz (right)
    # therefore <ϕ> = bottom_flux * t / Lz
    return mean(interior(field)) ≈ bottom_flux * model.clock.time / Lz
end

function fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT)
    Lz = 1
    κ₀ = FT(exp(-3))
    bz = FT(π)
    flux = - κ₀ * bz

    grid = RegularCartesianGrid(FT, size=(16, 16, 16), extent=(1, 1, Lz))

    buoyancy_bcs = TracerBoundaryConditions(grid, bottom=BoundaryCondition(Gradient, bz))
    κₑ_bcs = DiffusivityBoundaryConditions(grid, bottom=BoundaryCondition(Value, κ₀))
    model_bcs = (b=buoyancy_bcs, κₑ=(b=κₑ_bcs,))

    model = IncompressibleModel(
        grid=grid, architecture=arch, float_type=FT, tracers=:b, buoyancy=BuoyancyTracer(),
        closure=AnisotropicMinimumDissipation(), boundary_conditions=model_bcs
    )

    b₀(x, y, z) = z * bz
    set!(model, b=b₀)

    b = model.tracers.b
    mean_b₀ = mean(interior(b))

    τκ = Lz^2 / κ₀  # Diffusion time-scale
    Δt = 1e-6 * τκ  # Time step much less than diffusion time-scale
    Nt = 10         # Number of time steps

    for n in 1:Nt
        time_step!(model, Δt, euler= n==1)
    end

    # budget: Lz*∂<ϕ>/∂t = -Δflux = -top_flux/Lz (left) + bottom_flux/Lz (right)
    # therefore <ϕ> = bottom_flux * t / Lz
    #
    # Use an atol of 1e-6 so test passes with Float32 as there's a big cancellation
    # error due to buoyancy order of magnitude.
    #
    # Float32:
    # mean_b₀ = -1.5707965f0
    # mean(interior(b)) = -1.5708286f0
    # mean(interior(b)) - mean_b₀ = -3.20673f-5
    # (flux * model.clock.time) / Lz = -3.141593f-5
    #
    # Float64
    # mean_b₀ = -1.5707963267949192
    # mean(interior(b)) = -1.57082774272148
    # mean(interior(b)) - mean_b₀ = -3.141592656086267e-5
    # (flux * model.clock.time) / Lz = -3.141592653589793e-5
    return isapprox(mean(interior(b)) - mean_b₀, flux * model.clock.time / Lz, atol=1e-6)
end

@testset "Time stepping with boundary conditions" begin
    @info "Testing stepping with boundary conditions..."

    funbc(args...) = π
    boundaryfunbc = BoundaryFunction{:z, Face, Cell}((ξ, η, t) -> exp(ξ) * cos(η) * sin(t))

    @testset "Boundary condition instatiation and time-stepping" begin
        Nx = Ny = 16
        for arch in archs
	    ArrayType = array_type(arch)
            for FT in float_types
                @info "  Testing boundary condition instantiation and time-stepping [$(typeof(arch)), $FT]..."

                for fld in (:u, :v, :T, :S)
                    for bctype in (Gradient, Flux, Value)
                        arraybc = rand(FT, Nx, Ny) |> ArrayType
                        for bc in (FT(0.6), arraybc, funbc, boundaryfunbc)
                            @test test_z_boundary_condition_simple(arch, FT, fld, bctype, bc, Nx, Ny)
                        end
                    end

                    @test test_z_boundary_condition_array(arch, FT, fld)
                    @test test_flux_budget(arch, FT, fld)
                end
            end
        end
    end

    @testset "Custom diffusivity boundary conditions" begin
        for arch in archs, FT in float_types
            @info "  Testing flux budgets with custom diffusivity boundary conditions [$(typeof(arch)), $FT]..."
            @test fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT)
        end
    end
end


text:
######### Jags line program example  ###########

using Mamba, Jags

ProjDir = joinpath(dirname(@__FILE__))
cd(ProjDir) do

  line = "
  model {
    for (i in 1:n) {
          mu[i] <- alpha + beta*(x[i] - x.bar);
          y[i]   ~ dnorm(mu[i],tau);
    }
    x.bar   <- mean(x[]);
    alpha    ~ dnorm(0.0,1.0E-4);
    beta     ~ dnorm(0.0,1.0E-4);
    tau      ~ dgamma(1.0E-3,1.0E-3);
    sigma   <- 1.0/sqrt(tau);
  }
  "

  monitors = Dict(
    "alpha" => true,
    "beta" => true,
    "tau" => true,
    "sigma" => true
    )

  jagsmodel = Jagsmodel(
    name="line2",
    model=line,
    monitor=monitors,
    ncommands=4, nchains=1,
    #deviance=true, dic=true, popt=true,
    pdir=ProjDir
    );

  println("\nJagsmodel that will be used:")
  jagsmodel |> display

  data = Dict{String, Any}(
    "x" => [1, 2, 3, 4, 5],
    "y" => [1, 3, 3, 3, 5],
    "n" => 5
  )

  inits = [
    Dict("alpha" => 0,"beta" => 0,"tau" => 1,".RNG.name" => "base::Wichmann-Hill"),
    Dict("alpha" => 1,"beta" => 2,"tau" => 1,".RNG.name" => "base::Marsaglia-Multicarry"),
    Dict("alpha" => 3,"beta" => 3,"tau" => 2,".RNG.name" => "base::Super-Duper"),
    Dict("alpha" => 5,"beta" => 2,"tau" => 5,".RNG.name" => "base::Mersenne-Twister")
  ]

  println("Input observed data dictionary:")
  data |> display
  println("\nInput initial values dictionary:")
  inits |> display
  println()

  sim = jags(jagsmodel, data, inits, ProjDir)
  describe(sim)
  println()

  ## Brooks, Gelman and Rubin Convergence Diagnostic
  try
    gelmandiag(sim1, mpsrf=true, transform=true) |> display
  catch e
    #println(e)
    gelmandiag(sim, mpsrf=false, transform=true) |> display
  end

  ## Geweke Convergence Diagnostic
  gewekediag(sim) |> display

  ## Highest Posterior Density Intervals
  hpd(sim) |> display

  ## Cross-Correlations
  cor(sim) |> display

  ## Lag-Autocorrelations
  autocor(sim) |> display

  ## Deviance Information Criterion
  #dic(sim) |> display

  ## Plotting
  p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
  draw(p, nrow=4, ncol=4, filename="$(jagsmodel.name)-summaryplot", fmt=:svg)
  draw(p, nrow=4, ncol=4, filename="$(jagsmodel.name)-summaryplot", fmt=:pdf)

  # Below will only work on OSX, please adjust for your environment.
  # JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
  @static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
          for i in 1:4
            isfile("$(jagsmodel.name)-summaryplot-$(i).svg") &&
              run(`open -a $(JULIA_SVG_BROWSER) "$(jagsmodel.name)-summaryplot-$(i).svg"`)
          end
        end : println()

end #cd