text:
"Method to print the responses, for debug purposes"
function show_responses(responses)
    local shown::Set
    function subshow(prefix, value) #Non-dict
        println("** $prefix **:")
        println("\t",value)
    end

    function subshow(prefix, dict::Dict)
        for (key, value) in dict
            value ∈ shown && continue  # Skip any fields already shown
            push!(shown, value)
            subshow("$prefix $key", value)
        end
    end

    function subshow(prefix, vec::Vector)
        for entry in vec
            entry ∈ shown && continue  # Skip any fields already shown
            push!(shown, entry)
            subshow("$prefix - ", entry)
        end
    end

    println("Swift  Responses:")
    println("-------------------------")
    for resp in responses
        shown = Set()
        push!(shown, resp)
        subshow("", resp)
        println()
    end
end


text:
abstract type AbstractReadWriteLock <: Base.AbstractLock end

const NOTLOCKED = UInt64(0)
const NREADERS_INC = UInt64(2)
const WRITELOCK_MASK = UInt64(1)

const NReadersAndWritelock = UInt64

mutable struct ReadWriteLock <: AbstractReadWriteLock
    @atomic nreaders_and_writelock::NReadersAndWritelock
    # TODO: use condition variables with lock-free notify
    @const lock::ReentrantLock
    @const cond_read::Threads.Condition
    @const cond_write::Threads.Condition
end

const OFFSET_NREADERS_AND_WRITELOCK =
    fieldoffset_by_name(ReadWriteLock, :nreaders_and_writelock)

function ReadWriteLock()
    lock = ReentrantLock()
    cond_read = Threads.Condition(lock)
    cond_write = Threads.Condition(lock)
    return ReadWriteLock(NOTLOCKED, lock, cond_read, cond_write)
end

# Not very efficient but lock-free
function ConcurrentUtils.trylock_read(rwlock::ReadWriteLock; nspins = -∞, ntries = -∞)
    local ns::Int = 0
    local nt::Int = 0
    while true
        old = @atomic :monotonic rwlock.nreaders_and_writelock
        if iszero(old & WRITELOCK_MASK)
            # Try to acquire reader lock without the responsibility to receive or send the
            # notification:
            old, success = @atomicreplace(
                :acquire_release,
                :monotonic,
                rwlock.nreaders_and_writelock,
                old => old + NREADERS_INC,
            )
            success && return true
            nt += 1
            nt < ntries || return false
        end
        ns += 1
        ns < nspins || return false
    end
end

function ConcurrentUtils.lock_read(rwlock::ReadWriteLock)

    # Using hardware FAA
    ptr = Ptr{NReadersAndWritelock}(
        pointer_from_objref(rwlock) + OFFSET_NREADERS_AND_WRITELOCK,
    )
    GC.@preserve rwlock begin
        _, n = UnsafeAtomics.modify!(ptr, +, NREADERS_INC, acq_rel)
    end
    # n = @atomic :acquire_release rwlock.nreaders_and_writelock += NREADERS_INC

    if iszero(n & WRITELOCK_MASK)
        return
    end
    lock(rwlock.lock) do
        while true
            local n = @atomic :acquire rwlock.nreaders_and_writelock
            if iszero(n & WRITELOCK_MASK)
                @assert n > 0
                return
            end
            wait(rwlock.cond_read)
        end
    end
end

function ConcurrentUtils.unlock_read(rwlock::ReadWriteLock)

    # Using hardware FAA
    ptr = Ptr{NReadersAndWritelock}(
        pointer_from_objref(rwlock) + OFFSET_NREADERS_AND_WRITELOCK,
    )
    GC.@preserve rwlock begin
        _, n = UnsafeAtomics.modify!(ptr, -, NREADERS_INC, acq_rel)
    end
    # n = @atomic :acquire_release rwlock.nreaders_and_writelock -= NREADERS_INC

    @assert iszero(n & WRITELOCK_MASK)
    if iszero(n)
        lock(rwlock.lock) do
            notify(rwlock.cond_write; all = false)
        end
    end
    return
end

function Base.trylock(rwlock::ReadWriteLock)
    _, success = @atomicreplace(
        :acquire_release,
        :monotonic,
        rwlock.nreaders_and_writelock,
        NOTLOCKED => WRITELOCK_MASK,
    )
    return success::Bool
end

function Base.lock(rwlock::ReadWriteLock)
    if trylock(rwlock)
        return
    end
    lock(rwlock.lock) do
        while true
            if trylock(rwlock)
                return
            end
            wait(rwlock.cond_write)
        end
    end
end

function Base.unlock(rwlock::ReadWriteLock)
    @assert !iszero(rwlock.nreaders_and_writelock & WRITELOCK_MASK)
    @atomic :acquire_release rwlock.nreaders_and_writelock &= ~WRITELOCK_MASK
    lock(rwlock.lock) do
        notify(rwlock.cond_read)
        notify(rwlock.cond_write; all = false)
    end
    return
end

###
### High-level APIs
###

ConcurrentUtils.lock_read(lck) = lock(lck)
ConcurrentUtils.unlock_read(lck) = unlock(lck)

function ConcurrentUtils.lock_read(f, lock)
    lock_read(lock)
    try
        return f()
    finally
        unlock_read(lock)
    end
end

ext: jl


text:
#=
priority_queue:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-10-18
=#

module PQ
    using StatsBase

    buffer = []

    # move added node to root
    function upheap(buffer, n)
        while true
            p = Int(trunc((n - 1)/2))
            if (p < -1) || (buffer[p+1] <= buffer[n+1])
                break
            end
            tmp = buffer[n+1]
            buffer[n+1] = buffer[p+1]
            buffer[p+1] = tmp
            n = p
        end
    end

    # move root node to leaf
    function downheap(buffer, n)
        len = length(buffer)
        while true
            c = 2 * n + 1
            if c >= len
                break
            end
            # exchange parent and smaller child
            if (c + 1) < len
                if buffer[c+1] > buffer[c+2]
                    c += 1
                end
            end
            if buffer[n+1] <= buffer[c+1]
                break
            end
            tmp = buffer[n+1]
            buffer[n+1] = buffer[c+1]
            buffer[c+1] = tmp
            n = c
        end
    end

    # add data to tail of buffer
    # after adding, sort as upheap
    function push(buffer, data)
        push!(buffer, data)
        upheap(buffer, length(buffer)-1)
    end

    # get minimum data
    function peek(buffer)
        return buffer[1]
    end

    # get and remove minimum data
    function pop(buffer)
        min_val = buffer[1]
        last_val = buffer[end]
        deleteat!(buffer, length(buffer))
        if length(buffer) > 0
            # create heap again
            buffer[1] = last_val
            downheap(buffer, 0)
        end
        return min_val
    end

    function main()
        input = sample(1:100, 10, replace=false)
        # create by upheap
        for n in input
            push(buffer, n)
            println("Add $(n), Min data = $(peek(buffer))")
        end
        # remove minimum data and create by downheap
        while length(buffer) > 0
            print("$(pop(buffer)) ")
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .PQ
    PQ.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
###
# Standard QR Pivot Optimizer
###

mutable struct QRPivot
	Ψ::AbstractArray #Basis matrix from SVD, RPCA, etc.
	pivots::AbstractArray #Ranked list of sensor locations
end

function QRPivot(Ψ)
	n,m = size(Ψ)
	pivots = zeros((max(n,m),1))
	return QRPivot(Ψ,pivots)
end

function QRPivot(Ψ,n_sensors::Int)
	n,m = size(Ψ)
	pivots = zeros((n_sensors,1))
	return QRPivot(Ψ,pivots)
end

function QRPivot(Ψ,pivots::AbstractArray)
	return QRPivot(Ψ,pivots)
end

function fit(qr_pivot::QRPivot;optimizer_kwargs...)
	"""
	Fits the QRPivot sensor placement.
	qrpivot::QRPivot
		QRPivot object.
	optimizer_kwargs: dictionary
		Optional settings for optimizer
	"""
	qr_pivot.pivots = sensor_placement(qr_pivot.Ψ)
end

function sensor_placement(Ψ)
	# --> Compute the QRPivot w/ column pivoting decomposition of Ψ.
	_, _, p = qr(conj(Ψ), Val(true))
	return p[1:size(Ψ, 2)]
end


text:
struct RandomModel{T} <: Persa.Model{T}
    preference::Persa.Preference{T}
    users::Int
    items::Int
end

RandomModel(dataset::Persa.Dataset) = RandomModel(dataset.preference, Persa.users(dataset), Persa.items(dataset))

@testset "Model Learn Tests" begin
    df1 = createdummydatasetone()
    dataset1 = Persa.Dataset(df1)

    df2 = createdummydatasettwo()
    dataset2 = Persa.Dataset(df2, 10, 10)

    Persa.predict(model::RandomModel, user::Int, item::Int) = rand(model.preference.possibles)

    model1 = RandomModel(dataset1)
    model2 = RandomModel(dataset2)

    @testset "Dummy Tests" begin
        @test model1[1,1] > 0

        values = model1[dataset1]
        @test length(values) == length(dataset1)

        values = model2[dataset2]
        @test length(values) == length(dataset2)

        @test ismissing(dataset1[3,6])
    end
end


text:
module TimeSeriesTools

using TimeSeries

export  is_equally_spaced, is_complete, is_clean,
        require_equally_spaced, require_complete, require_clean,
        plot, histogram, density,
        directedscatter, laghistogram, acf, pacf, ccf,
        cumulative_periodogram

include("data_quality.jl")
include("visualizations.jl")
include("stat_tests.jl")

end # module


text:
include("cell/cell.jl")

using .Model

basic_protocol = JSON.parsefile("../data/com_protocol_basic.json")

include("build_neuron.jl")

# A simulation is broken up into Spans. This allows notifications to
# be sent to the client.
# Each Span (aka span_time) equals: Duration where
# a span contains samples from all synapses.
# Duration is specified in microseconds.
#
# |----------------------------- Duration ----------------------------|
# |----- Span -----|----- Span -----|----- Span -----|----- Span -----|

# Total simulation time = Sum of all spans or Spans * Duration

function build_simulation!(model::Model.ModelData, samples::Model.Samples, streams::Streams)
    # Build neuron for simulation.

    println("Building neuron model.")
    cell = build(model, samples, streams)

    # Initialize cell just once based on Model data
    initialize!(cell)

    cell
end

# -_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_
# NOTE: simulate() is called from run.jl
# -_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_-_---_
# Once the server is started it builds the components to support the simualtion.
# Simulate does not create anything it simply resets for a simulation run.
function simulate(chan::Channel{String},
    model::Model.ModelData, samples::Model.Samples, streams::Streams,
    cell::Cell)
    println("Configuring simulation...")

    reset!(cell)

    span_time = Model.span_time(model)

    # ~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--
    # Now the simulation starts.
    # ~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--
    spans = Model.spans(model)

    # duration = Model.duration(model)

    println("### --------- Starting Simulation -------- ###")

    # Master time.
    t = 1

    # Run each span
    for span in 1:spans

        # "t" is a single time step representing a single TimeScale.
        # If the TimeScale = 100us then if t = 2 then 100us passed
        # and when t = 3 then 200us passed etc.
        for span_t in 1:span_time  # a single tick of the simulation for one span
            # This is the core of the simulation.
            integrate!(cell, span_t, t)

            # Collecting is centralized in streams.jl for consistency.
            # Collect all data for analysis by client.
            collect!(streams, samples, span_t)

            # Exercise all stimulus which means all merger streams.
            # This causes each stream to update and move its internal value to its output.
            exercise!(streams)

            t += 1
        end

        # We have finished a span, write it out to disk
        Model.write_samples!(samples, model, span)

        # Notify client that a span completed.
        notify_client_span_completed(chan, span)

        # Yield for channel tasks
        yield()

        # We "reset" the samples for the next span.
        Model.reset_samples!(samples)
    end

    println("### --------- Simulaton Complete -------- ###")
end

function reset!()
end

function notify_client_span_completed(chan::Channel{String}, span::Int64)
    data = basic_protocol

    # Populate
    data["From"] = "Simulation"
    data["To"] = "Client"
    data["Type"] = "Status"
    data["Data"] = "Span Completed"
    data["Data1"] = "poi_sample_c$span.data"
    data["Data2"] = span

    put!(chan, JSON.json(data))
end


text:
using KernelRidgeRegression
using MLKernels
using Base.Test
using StatsBase

@test GaussianKernel(3.0) == GaussianKernel(3.0)

N = 5000
x = rand(1, N) * 4π - 2π
yy = sinc.(x) # vec(sinc.(4 .* x) .+ 0.2 .* sin.(30 .* x))
y = squeeze(yy + 0.1randn(1, N), 1)
xnew = collect(-2.5π:0.01:2.5π)'

mykrr = fit(KRR, x, y, 1e-3/5000,
            GaussianKernel(100.0))
ynew = predict(mykrr, xnew);

mynystkrr = fit(NystromKRR,
                x, y, 10.0, 280,
                GaussianKernel(100.0))
ynystnew = predict(mynystkrr, xnew)

mysrkrr = fit(SubsetRegressorsKRR,
                x, y, 1.0, 280,
                GaussianKernel(100.0))
ysrnew = predict(mysrkrr, xnew)

myfastkrr = fit(FastKRR,
                x, y, 4/5000, 11,
                GaussianKernel(100.0))
yfastnew = predict(myfastkrr, xnew)

myfastkrr2 = fitPar(FastKRR, length(x),
                    (i) -> x[:, i], (i) -> y[i],
                    4/5000, 11, GaussianKernel(100.0))
yfastnew2 = predict(myfastkrr, xnew)

mytnkrr = fit(TruncatedNewtonKRR,
              x, y, 4/5000, GaussianKernel(100.0),
              0.5, 200)
ytnnew = predict(mytnkrr, xnew)

myrandkrr = fit(RandomFourierFeatures,
                x, y, 1/500.0, 500, 1.0)
yrandnew = predict(myrandkrr, xnew)

emean = sqrt(mean((vec(sinc.(xnew)) - mean(xnew)) .^ 2))
ekrr  = sqrt(mean((vec(sinc.(xnew)) - ynew)       .^ 2))
enyst = sqrt(mean((vec(sinc.(xnew)) - ynystnew)   .^ 2))
esr   = sqrt(mean((vec(sinc.(xnew)) - ysrnew)     .^ 2))
efast = sqrt(mean((vec(sinc.(xnew)) - yfastnew)   .^ 2))
erand = sqrt(mean((vec(sinc.(xnew)) - yrandnew)   .^ 2))
etn   = sqrt(mean((vec(sinc.(xnew)) - ytnnew)     .^ 2))

@test eltype(emean) == Float64
@test eltype(ekrr ) == Float64
@test eltype(enyst) == Float64
@test eltype(esr)   == Float64
@test eltype(efast) == Float64
@test eltype(erand) == Float64
@test eltype(etn  ) == Float64

@test emean > ekrr
@test emean > enyst
@test emean > esr
@test emean > efast
@test emean > erand
@test emean > etn

# from julia/test/show.jl:
replstr(x) = sprint((io, y′) -> show(IOContext(io, :limit => true), MIME("text/plain"), y′), x)

@test replstr(mykrr) == "KernelRidgeRegression.KRR{Float64}:\n    λ = 2.0e-7\n    ϕ = SquaredExponentialKernel(100.0)"
@test contains(replstr(myrandkrr), "KernelRidgeRegression.RandomFourierFeatures{Float64,Complex{Float64}}:\n    λ = 0.002\n:    σ = 1.0\n:    K = 500\n    ϕ = KernelRidgeRegression")
@test replstr(mynystkrr) == "KernelRidgeRegression.NystromKRR{Float64}:\n    λ = 10.0\n    ϕ = SquaredExponentialKernel(100.0)\n    m = 280"
@test replstr(mysrkrr) == "KernelRidgeRegression.SubsetRegressorsKRR{Float64}:\n    λ = 1.0\n    ϕ = SquaredExponentialKernel(100.0)\n    m = 280"
@test replstr(myfastkrr) == "KernelRidgeRegression.FastKRR{Float64}:\n    λ = 0.0008\n    m = 11\n    ϕ = SquaredExponentialKernel(100.0)"
@test replstr(myfastkrr2) == "KernelRidgeRegression.FastKRR{Float64}:\n    λ = 0.0008\n    m = 11\n    ϕ = SquaredExponentialKernel(100.0)"
@test replstr(mytnkrr) == "KernelRidgeRegression.TruncatedNewtonKRR{Float64}:\n    λ = 0.0008\n    ϕ = SquaredExponentialKernel(100.0)"


text:
export read_results!

"""

    read_results!(optimizer::Optimizer{T}, filepath::String, redundant_entries::Vector)
Populates `optimizer` with results in a SDPA-formatted output file specified by `filepath`. Redundant entries corresponding to linearly dependent constraints are set to 0.

"""
function read_results!(
    optimizer::Optimizer{T},
    filepath::String,
    redundant_entries::Vector
) where T

    endswith(filepath, ".dat") || error("Filename '$filepath' must end with .dat")
    getnextline(io::IO) =
        eof(io) ?
        error("The output file is possibly corrupted. Check that $filepath conforms to the SDPA output format.") :
        chomp(readline(io))

    function replace_brackets!(str::SubString)
        str = replace(str, "{" => "[")
        str = replace(str, "}" => "]")
        return str
    end
    function remove_brackets!(str)
        str = replace(str, "{" => "")
        str = replace(str, "}" => "")
        return str
    end

    phasevalue = :noINFO
    objValPrimalstring = ""
    objValDualstring = ""
    xVecstring = ""
    xMatvec = T[]
    yMatvec = T[]

    file = open(filepath, "r") do io
            line = getnextline(io)

            while !startswith(line, "phase.value")
                line = getnextline(io)
            end
            optimizer.phasevalue = Symbol(split(line)[3])

            while !startswith(line, "objValPrimal")
                line = getnextline(io)
            end
            objValPrimalstring = split(line)[3]

            while !startswith(line, "objValDual")
                line = getnextline(io)
            end
            objValDualstring = split(line)[3]

            while !startswith(line, "xVec")
                line = getnextline(io)
            end
            line = getnextline(io)
            xVecstring = line

            while !startswith(line, "xMat")
                line = getnextline(io)
            end
            line = getnextline(io)
            line = getnextline(io)
            while !startswith(line, "}")
                xMatstring = remove_brackets!(line)
                if endswith(xMatstring, ",")
                    append!(xMatvec, parse.(T, split(xMatstring[1:end-1], ",")))
                else
                    append!(xMatvec, parse.(T, split(xMatstring, ",")))
                end
                line = getnextline(io)
            end
            line = getnextline(io)
            line = getnextline(io)
            line = getnextline(io)
            while !startswith(line, "}")
                yMatstring = remove_brackets!(line)
                if endswith(yMatstring, ",")
                    append!(yMatvec, parse.(T, split(yMatstring[1:end-1], ",")))
                else
                    append!(yMatvec, parse.(T, split(yMatstring, ",")))
                end
                line = getnextline(io)
            end
        end
    xVecstring = remove_brackets!(xVecstring)
    xVecstring = split(xVecstring, ",")
    xVec = parse.(T, xVecstring)
    optimizer.primalobj = parse(T, objValPrimalstring)
    optimizer.dualobj = parse(T, objValDualstring)
    for i in redundant_entries
        splice!(xVec, i:i-1, zero(T))
    end
    optimizer.y = xVec
    structurevec = optimizer.blockdims
    yMatbm = PrimalSolution{T}(map(
        n -> zeros(T, abs(n), abs(n)),
        optimizer.blockdims
    ))
    xMatbm = VarDualSolution{T}(map(
        n -> zeros(T, abs(n), abs(n)),
        optimizer.blockdims
    ))
    for i = 1:length(structurevec)
        dim = structurevec[i]
        if dim < 0
            dim = abs(dim)
            xblock = xMatvec[1:dim]
            deleteat!(xMatvec, 1:dim)
            xblock = Diagonal(xblock)
            xMatbm.blocks[i] = xblock

            yblock = yMatvec[1:dim]
            deleteat!(yMatvec, 1:dim)
            yblock = Diagonal(yblock)
            yMatbm.blocks[i] = yblock

        elseif dim > 0
            xblock = xMatvec[1:dim^2]
            deleteat!(xMatvec, 1:dim^2)
            xblock = reshape(xblock, dim, dim)
            xMatbm.blocks[i] = xblock

            yblock = yMatvec[1:dim^2]
            deleteat!(yMatvec, 1:dim^2)
            yblock = reshape(yblock, dim, dim)
            yMatbm.blocks[i] = yblock

        end
    end
    optimizer.X = yMatbm
    optimizer.Z = xMatbm

end

"""

    inputElement(optimizer::Optimizer, constr_number::Int, blk::Int, i::Int, j::Int, value::T) where T

Stores the constraint data in `optimizer.elemdata` as a vector of tuples. Each tuple corresponds to one line in the SDPA-formatted input file.
"""
function inputElement(
    optimizer::Optimizer,
    constr_number::Int,
    blk::Int,
    i::Int,
    j::Int,
    value::T
) where T
    push!(optimizer.elemdata, (constr_number, blk, i, j, value))
end

"""

    initializeSolve(optimizer::Optimizer)

Writes problem data into an SDPA-formatted file named `input.dat-s`. `presolve.jl` routine is applied as indicated by `optimizer.presolve`.

Returns a vector of indices for redundant constraints, which are omitted from the input file.
"""
function initializeSolve(optimizer::Optimizer)
    if !optimizer.presolve
        filename = joinpath(optimizer.tempfile, "input.dat-s")
        file = open(filename, "w") do io
            nconstrs = length(optimizer.b)
            nblocks = length(optimizer.blockdims)
            println(io, nconstrs)
            println(io, nblocks)
            writedlm(io, optimizer.blockdims', " ")
            writedlm(io, optimizer.b', " ")
            for line in optimizer.elemdata
                for i in line
                    print(io, i)
                    print(io, " ")
                end
                println(io)
            end
        end
        return []
    else
        redundant_F = presolve(optimizer)
        reduced = joinpath(optimizer.tempfile, "input.dat-s")
        file =
            open(reduced, "w") do io
                nconstrs = length(optimizer.b) - length(redundant_F)
                nblocks = length(optimizer.blockdims)
                println(io, nconstrs)
                println(io, nblocks)
                writedlm(io, optimizer.blockdims', " ")
                cVec = deleteat!(copy(optimizer.b), redundant_F)
                writedlm(io, cVec', " ")
                for entry in optimizer.elemdata
                    constr_index, blk, i, j, value = entry
                    if !in(constr_index, redundant_F)
                        constr_index = constr_index -
                                       count(x -> x < constr_index, redundant_F)
                        for i in (constr_index, blk, i, j, value)
                            print(io, i)
                            print(io, " ")
                        end
                        println(io)
                    end
                end
            end
        return redundant_F
    end
end


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Tests that all fields of these two structs are the same.

# Arguments
- `struct1`: a struct
- `struct2`: a struct of the same type as struct1
"""
function test_fields_equal(
  struct1::TypeName,
  struct2::TypeName,
) where {TypeName}
  @testset "$field_name matches" for field_name in fieldnames(TypeName)
    @test getfield(struct1, field_name) == getfield(struct2, field_name)
  end
end

function test_fields_equal(
  struct1::FirstOrderLp.IterationStats,
  struct2::FirstOrderLp.IterationStats,
)
  @testset "$field_name matches" for field_name in
                                     fieldnames(FirstOrderLp.IterationStats)
    struct1_field = getfield(struct1, field_name)
    struct2_field = getfield(struct2, field_name)
    if field_name == :convergence_information ||
       field_name == :infeasibility_information
      @test length(struct1_field) == length(struct2_field)
      for (left_entry, right_entry) in
          Iterators.zip(struct1_field, struct2_field)
        test_fields_equal(left_entry, right_entry)
      end
    else
      @test struct1_field == struct2_field
    end
  end
end

"""
Tests that all fields of these two structs are approximately the same.

# Arguments
- `struct1`: a struct
- `struct2`: a struct of the same type as struct1
"""
function test_fields_approx_equal(
  struct1::TypeName,
  struct2::TypeName,
) where {TypeName}
  @testset "$field_name matches" for field_name in fieldnames(TypeName)
    @test getfield(struct1, field_name) ≈ getfield(struct2, field_name)
  end
end

"""
Tests all entries of two tuples are approximately same.
"""
function test_tuple_approx_equal(tuple1, tuple2)
  @test length(tuple1) == length(tuple2)
  @testset "tuple index $i" for i in 1:length(tuple2)
    @test tuple1[i] ≈ tuple2[i]
  end
end

"""
Terminate based on iteration limit only.
The termination tolerance are set to zero to stop any tests from
breaking. This is because when the tests were orginally designed when there
was no termination criteria.
"""
function terminate_on_iteration_limit(n::Int64)
  termination_criteria = FirstOrderLp.construct_termination_criteria(
    optimality_norm = FirstOrderLp.L_INF,
    eps_optimal_absolute = 0.0,
    eps_optimal_relative = 0.0,
    eps_primal_infeasible = 0.0,
    eps_dual_infeasible = 0.0,
    time_sec_limit = 100.0,
    iteration_limit = n,
    kkt_matrix_pass_limit = Inf,
  )
  return termination_criteria
end


text:
write_test_output = false

path = dirname(@__FILE__)

if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

###################################################################
# Test: initial_draw!()
###################################################################
include("modelsetup.jl")

m = setup_linear_model()

# Read in generated data
data = h5read("reference/test_data.h5", "data")
X = h5read("reference/test_data.h5", "X")

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

####################################################################
init_cloud = SMC.Cloud(length(m.parameters), get_setting(m,:n_particles))

@everywhere Random.seed!(42)
SMC.initial_draw!(loglik_fn, m.parameters, data, init_cloud)

if write_test_output
    JLD2.jldopen(string("reference/initial_draw_out_version=", ver, ".jld2"), "w") do file
        write(file, "cloud", init_cloud)
    end
end

saved_init_cloud = load(string("reference/initial_draw_out_version=", ver, ".jld2"), "cloud")

@testset "Initial draw" begin
    @test @test_matrix_approx_eq SMC.get_vals(init_cloud) SMC.get_vals(saved_init_cloud)
    @test @test_matrix_approx_eq SMC.get_loglh(init_cloud) SMC.get_loglh(saved_init_cloud)
end

###################################################################
# Test: one_draw()
###################################################################
draw = SMC.one_draw(loglik_fn, m.parameters, data)

if write_test_output
    JLD2.jldopen(string("reference/one_draw_out_version=", ver, ".jld2"), true, true, true, IOStream) do file
        file["draw"] = draw
    end
end

test_draw = JLD2.jldopen(string("reference/one_draw_out_version=", ver, ".jld2"), "r") do file
    file["draw"]
end

###################################################################
@testset "One draw" begin
    @test draw[1] == test_draw[1]
    @test draw[2] ≈ test_draw[2]
end

###################################################################
# Test: draw_likelihood()
###################################################################
draw_lik = SMC.draw_likelihood(loglik_fn, m.parameters, data, vec(draw[1]))

if write_test_output
    JLD2.jldopen(string("reference/draw_likelihood_out_version=", ver, ".jld2"), true, true, true, IOStream) do file
        file["draw_lik"] = draw_lik
    end
end
test_draw_lik = JLD2.jldopen(string("reference/draw_likelihood_out_version=", ver, ".jld2"), "r") do file
    file["draw_lik"]
end

###################################################################
@testset "Draw Likelihood" begin
    @test draw_lik[1][1] ≈ test_draw_lik[1][1]
    @test draw_lik[2] == test_draw_lik[2]
end

###################################################################
# Test: initialize_likelihoods!()
###################################################################
SMC.initialize_likelihoods!(loglik_fn, m.parameters, data, init_cloud)

if write_test_output
    JLD2.jldopen(string("reference/initialize_likelihood_out_version=", ver, ".jld2"), true, true, true, IOStream) do file
        file["init_lik_cloud"] = init_cloud
    end
end
test_init_cloud = JLD2.jldopen(string("reference/initialize_likelihood_out_version=", ver, ".jld2"), "r") do file
    file["init_lik_cloud"]
end

###################################################################
@testset "Initialize Likelihoods" begin
    @test @test_matrix_approx_eq SMC.get_vals(init_cloud) SMC.get_vals(test_init_cloud)
    @test @test_matrix_approx_eq SMC.get_loglh(init_cloud) SMC.get_loglh(test_init_cloud)
end

###################################################################
# Test: initialize_cloud_settings!()
###################################################################
SMC.initialize_cloud_settings!(init_cloud)

if write_test_output
    JLD2.jldopen(string("reference/initialize_cloud_settings_version=", ver, ".jld2"), true, true, true, IOStream) do file
        file["init_cloud"] = init_cloud
    end
end
test_init_cloud = JLD2.jldopen(string("reference/initialize_cloud_settings_version=", ver, ".jld2"), "r") do file
    file["init_cloud"]
end

###################################################################
@testset "Initialize Cloud" begin
    @test init_cloud.ESS         == test_init_cloud.ESS
    @test init_cloud.stage_index == test_init_cloud.stage_index
    @test init_cloud.n_Φ         == test_init_cloud.n_Φ
    @test init_cloud.resamples   == test_init_cloud.resamples
    @test init_cloud.c           == test_init_cloud.c
    @test init_cloud.total_sampling_time == test_init_cloud.total_sampling_time
    @test init_cloud.tempering_schedule  == test_init_cloud.tempering_schedule
end


text:
# Generic functions that apply to all (or most) objects

set_pref_scale!(switches :: AbstractEntrySwitches{F1},
    prefScale :: F1) where F1 =
    switches.entryPrefScale = prefScale;

"""
	$(SIGNATURES)

Scale entry probabilities by (j,l) -> c to bound away from 0.
Bound entry probs by (j,l) away from 1.
Does not guarantee particular min values (that is hard to do).
"""
function scale_entry_probs!(entryProb_jlcM :: Array{F1, 3},
    minEntryProb :: F1, maxEntryProb :: F1) where F1 <: Real

    J, nl, nc = size(entryProb_jlcM);
    for iType = 1 : J
        for l = 1 : nl
            prob_cV =
                max.(minEntryProb, entryProb_jlcM[iType, l, :]);

            pSum = sum(prob_cV);
            if pSum > maxEntryProb
                prob_cV .*= (maxEntryProb / pSum);
            end
            entryProb_jlcM[iType, l, :] .= prob_cV;
        end
        # for ic = 1 : nc
        #     # Scale so that total entry prob j -> c is at least minEntryProb
        #     prob_lV = entryProb_jlcM[iType, :, ic];
        #     pSum = sum(prob_lV);
        #     if pSum < 0.000001
        #         # Basically all zero entries. Set to a constant
        #         prob_lV = fill(minEntryProb / nl, nl);
        #     elseif pSum < minEntryProb
        #         prob_lV .*= (minEntryProb / pSum);
        #     end
        #     entryProb_jlcM[iType, :, ic] = prob_lV;
        # end

        # # Scale so that entry prob j <= maxEntryProb
        # pSum = sum(entryProb_jlcM[iType,:,:]);
		# if pSum > maxEntryProb
		# 	entryProb_jlcM[iType,:,:] .*= (maxEntryProb / pSum);
        # end
	end
	return nothing
end

## ---------------  Constructors

"""
    `init_entry_decision(objId :: ObjectId, switches :: AbstractEntrySwitches)`

Initializes an `AbstractEntryDecision` from its switches.
"""
function init_entry_decision end

# Initialize entry preference scale parameter.
function init_entry_prefscale(switches :: AbstractEntrySwitches{F1},
    st :: SymbolTable) where F1
    entryPrefScale = switches.entryPrefScale;
    pEntryPref = Param(:entryPrefScale,
    LatexLH.description(st, :pScaleEntry), latex(st, :pScaleEntry), entryPrefScale, entryPrefScale, F1(0.1), F1(3.0),
        switches.calEntryPrefScale);
end

## --------------  Solve entry decisions

"""
    $(SIGNATURES)

Entry probability for a student who is admitted to colleges in `admitV`.
Returns: Entry prob by [type, college], expected value at decision stage by type.
This function does not handle the sequential nature of admissions. It is mainly here for unified interface.
It defaults to the one step entry protocol. If that does not apply for a protocol, need to define a new method (e.g. two step entry).
"""
function entry_probs(e :: AbstractEntryDecision{F1},
    vWork_jV :: Vector{F1}, vCollege_jcM :: Matrix{F1}, admitV;
    prefShocks :: Bool = true) where F1 <: Real

    if prefShocks
        prefScale = entry_pref_scale(e);
    else
        prefScale = zero(F1);
    end
    return one_step_entry_probs(prefScale, vWork_jV, vCollege_jcM, admitV);
end

# The same for one student
function entry_probs(e :: AbstractEntryDecision{F1},
    vWork :: F1, vCollege_cV :: Vector{F1}, admitV;
    prefShocks :: Bool = true) where F1

    if prefShocks
        prefScale = entry_pref_scale(e);
    else
        prefScale = zero(F1);
    end
    return one_step_entry_probs(prefScale, vWork, vCollege_cV, admitV)
end

# The same shutting down preference shocks
# function entry_probs_no_pref_shocks(e :: AbstractEntryDecision{F1},
#     vWork :: F1, vCollege_cV :: Vector{F1}, admitV) where F1

#     return one_step_entry_probs(zero(F1), vWork, vCollege_cV, admitV)
# end

# Generic entry decision. One step. Given pref shock scale.
# Works for any entry protocol where entry works in one step (work/study and which college at the same time).
function one_step_entry_probs(entryPrefScale :: F1,
    vWork_jV :: Vector{F1}, vCollege_jcM :: Matrix{F1}, admitV) where F1 <: Real

    J, nc = size(vCollege_jcM);
    prob_jxM = zeros(F1, J, nc);
    if isempty(admitV)
        eVal_jV = copy(vWork_jV);
    else
        if entryPrefScale > 0.00001
            # Prob of work in column 1. Then admitted colleges.
            probM, eVal_jV = EconLH.extreme_value_decision(
                hcat(vWork_jV, vCollege_jcM[:, admitV]),
                entryPrefScale; demeaned = true);
            for j = 1 : J
                probV = probM[j,:];
                make_valid_probs!(probV);
                prob_jxM[j, admitV] .= probV[2 : end];
            end
        else
            eVal_jV, icV = max_choices(vWork_jV, vCollege_jcM, admitV);
            for j = 1 : J
                if icV[j] > 0
                    prob_jxM[j, icV[j]] = one(F1);
                end
            end
        end
    end
    return prob_jxM, eVal_jV
end

# The same for one type.
# Preference shocks may be zero.
function one_step_entry_probs(entryPrefScale :: F1,
    vWork :: F1, vCollege_cV :: Vector{F1}, admitV) where F1 <: Real

    nc = length(vCollege_cV);
    prob_cV = zeros(F1, nc);
    if isempty(admitV)
        eVal = vWork;
    else
        if entryPrefScale > 0.00001
            probV, eVal = EconLH.extreme_value_decision_one(
                vcat(vWork, vCollege_cV[admitV]), entryPrefScale;
                demeaned = true);
            make_valid_probs!(probV);
            prob_cV[admitV] .= probV[2 : end];
        else
            # No preference shocks - just choose the best option
            eVal, ic = max_choice(vWork, vCollege_cV, admitV);
            if ic > 0
                prob_cV[ic] = one(F1);
            end
        end
    end
    return prob_cV, eVal
end

"""
	$(SIGNATURES)

Entry probs for one student across all admissions sets. Handles the case where some colleges are full.
Always for multiple locations (matrix inputs). But one location is allowed.

# Outputs
- `entryProb_clM`: probability of entering each college
- `eVal`: expected value
- `entryProbBest_clM`: fraction of students who enter college `c, l` and for who `c` is the best available college (in any location). For the top college, this is identical to the entry rate.
"""
function entry_decisions_one_student(entryS :: AbstractEntryDecision{F1},
    admissionS :: AbstractAdmissionsRule{I1, F1},
    admProbFct :: AF1,
    vWork :: F1, vCollege_cV :: AbstractVector{F1},
    endowPct :: F1, full_clM :: AbstractMatrix{Bool}, l :: Integer;
    prefShocks :: Bool = true)  where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1}

    nl = n_locations(entryS);
    nc = n_colleges(entryS);
    @check size(full_clM) == (nc, nl)

    # Value of college is the same for all locations, except local
    vCollege_clM = repeat(vCollege_cV, outer = (1, nl));
    vCollege_clM[:, l] .+= value_local(entryS);

    entryProb_clM = zeros(F1, nc, nl);
    eVal = zero(F1);
    # Fraction of students who attend the best available college (regardless of location). By quality attended. For the best college, this equals `entryProb_cV`.
    entryProbBest_clM = zeros(F1, nc, nl);
    # Admission rule gives admission to one college type in all locations
    for (iSet, admitV) in enumerate(admissionS)
        # Prob that each person draws this college set
        probSet = prob_coll_set(admissionS, admProbFct, iSet, endowPct);
        avail_clM = available_colleges(entryS, full_clM, admitV, l);
        @check size(avail_clM) == size(vCollege_clM)

        # Entry probs for this set
        prob_clV, eValSet =
            entry_probs(entryS, vWork, vec(vCollege_clM), vec(avail_clM);
                prefShocks = prefShocks);
        # Tested separately that `reshape` undoes `vec`
        prob_clM = reshape(prob_clV, nc, nl);
        entryProb_clM .+= probSet .* prob_clM;
        eVal += probSet * eValSet;

        # Fraction attending best available.
        cBest = best_available(avail_clM);
        entryProbBest_clM[cBest,:] .+= probSet .* prob_clM[cBest, :];
    end

    make_valid_probs!(entryProb_clM);
    make_valid_probs!(entryProbBest_clM);
    # Deal with rounding errors that arise during scaling.
    bracket_array!(entryProbBest_clM, zeros(F1, nc, nl), entryProb_clM);
    return entryProb_clM, eVal, entryProbBest_clM
end

# Can only attend colleges that are not full, where student is admitted.
function available_colleges(entryS :: AbstractEntryDecision{F1},
    full_clM :: AbstractMatrix{Bool},
    admitV, l :: Integer) where F1

    # Do not use falses here. It creates a BitArray
    avail_clM = fill(false, size(full_clM));
    # Admitted counts for all locations
    avail_clM[admitV, :] .= true;
    avail_clM[full_clM] .= false;
    # If there are local-only colleges, mark those as not available
    mark_local_only_colleges!(entryS, avail_clM, l);
    return avail_clM
end

function best_available(avail_clM)
    cBest = 1;
    for ic = 2 : size(avail_clM, 1)
        if any(avail_clM[ic,:])
            cBest = ic;
        end
    end
    return cBest
end

# Mark local only colleges as unavailable when not local.
function mark_local_only_colleges!(entryS :: AbstractEntryDecision{F1},
    avail_clM, l) where F1

    idxV = local_only_colleges(entryS);
    if !isempty(idxV)
        for ic ∈ idxV
            aLocal = avail_clM[ic, l];
            avail_clM[ic, :] .= false;
            avail_clM[ic, l] = aLocal;
        end
    end
    return nothing
end

# The same with one location. Simply calls the multi-location code with one location dimension.
function entry_decisions_one_student(entryS :: AbstractEntryDecision{F1},
    admissionS :: AbstractAdmissionsRule{I1, F1},
    admProbFct :: AF1,
    vWork :: F1, vCollege_cV :: AbstractVector{F1},
    endowPct :: F1, full_cV :: AbstractVector{Bool};
    prefShocks :: Bool = true)  where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1}

    # Just solve the multi-location version for the first location.
    entryProb_clM, eVal, entryProbBest_clM =
        entry_decisions_one_student(entryS, admissionS, admProbFct,
        vWork, vCollege_cV, endowPct,  repeat(full_cV, outer = (1,1)), 1;
        prefShocks = prefShocks);

    make_valid_probs!(entryProb_clM)
    return vec(entryProb_clM), eVal, entryProbBest_clM
end

## ------------  Implied outcomes

# """
# 	$(SIGNATURES)

# Compute college enrollment from type mass and entry probabilities.
# """
# function college_enrollment(entryProb_jcM :: Matrix{F1},
#     typeMass_jV) where F1 <: Real

#     J, nc = size(entryProb_jcM);
#     enrollV = zeros(F1, nc);
#     for ic = 1 : nc
#         enrollV[ic] = sum(entryProb_jcM[:, ic] .* typeMass_jV);
#     end
#     return enrollV
# end

# function college_enrollment(e :: AbstractEntryDecision{F1},
#     entryProb_jcM :: Matrix{F1}) where F1 <: Real
#     return college_enrollment(e.switches, entryProb_jcM);
# end

# function college_enrollment(e :: AbstractEntrySwitches{F1},
#     entryProb_jcM :: Matrix{F1}) where F1 <: Real

#     @assert n_locations(e) == 1
#     J = size(entryProb_jcM, 1);
#     return college_enrollment(entryProb_jcM, type_masses(e));
# end

# function college_enrollment(e :: AbstractEntryDecision{F1},
#     entryProb_cV :: Vector{F1},  j :: Integer) where F1

#     @assert n_locations(e) == 1
#     return entryProb_cV .* type_mass(e, j);
# end

# """
#     $(SIGNATURES)

# Return `Bool` vector that indicates which colleges are full. Only matters for entry structures with capacity constraints.
# """
# function colleges_full(e :: AbstractEntryDecision{F1},
#     entryProb_jcM :: Matrix{F1}) where F1 <: Real

#     @assert n_locations(e) == 1  "Not valid for multiple locations"
#     return college_enrollment(e, entryProb_jcM) .>= capacities(e)
# end

# --------------


text:
cfs = [:s2, :l2, :c2, :surfsurf, :surfvoid]
dirs = Dict(
    1 => [:x],
    2 => [:x, :y, :xy, :yx],
    3 => [:x, :y, :z]
)

function test_cf_on_img(img, cf_map, cf_dir)
    for periodic in [true, false]
        cmap = cf_map(img; periodic)
        cdir = cf_dir(img, true; len=size(img, 1), periodic, directions=dirs[ndims(img)])

        for (direction, data) in cdir
            @test Map.dir_from_map(cmap, direction) ≈ data
        end
    end
end

function test_all()
    dims = [1, 2, 3]
    n = 10

    for dim in dims
        ns = Tuple(ones(Int, dim) * n)
        img = rand(Bool, ns)

        for cf in cfs
            @testset "random $(dim)D image, $(cf)" begin
                test_cf_on_img(img, Map.eval(cf), Directional.eval(cf))
            end
        end
    end
end

test_all()


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

mutable struct GenericIterator{N} end
Base.iterate(::GenericIterator{N}, i=1) where {N} = i > N ? nothing : (i, i + 1)
Base.IteratorSize(::Type{GenericIterator{N}}) where {N} = Base.SizeUnknown()

function generic_map_tests(mapf, inplace_mapf=nothing)
    for typ in (Float16, Float32, Float64,
                Int8, Int16, Int32, Int64, Int128,
                UInt8, UInt16, UInt32, UInt64, UInt128),
        arg_typ in (Integer,
                    Signed,
                    Unsigned)
        X = typ[1:10...]
        _typ = typeof(arg_typ(one(typ)))
        @test mapf(arg_typ, X) == _typ[1:10...]
    end

    # generic map
    f(x) = x + 1
    I = GenericIterator{10}()
    @test mapf(f, I) == Any[2:11...]

    # AbstractArray map for 2 arg case
    f(x, y) = x + y
    B = Float64[1:10...]
    C = Float64[1:10...]
    @test mapf(f, convert(Vector{Int},B), C) == Float64[ 2 * i for i in 1:10 ]
    @test mapf(f, Int[], Float64[]) == Union{}[]
    # map with different result types
    let m = mapf(x->x+1, Number[1, 2.0])
        @test isa(m, Vector{Real})
        @test m == Real[2, 3.0]
    end

    # AbstractArray map for N-arg case
    A = Array{Int}(undef, 10)
    f(x, y, z) = x + y + z
    D = Float64[1:10...]

    @test map!(f, A, B, C, D) == Int[ 3 * i for i in 1:10 ]
    @test mapf(f, B, C, D) == Float64[ 3 * i for i in 1:10 ]
    @test mapf(f, Int[], Int[], Complex{Int}[]) == Union{}[]

    # In-place map
    if inplace_mapf != nothing
        A = Float64[1:10...]
        inplace_mapf(x -> x*x, A, A)
        @test A == map(x -> x*x, Float64[1:10...])

        # Map to destination collection
        B = inplace_mapf((x,y,z)->x*y*z, A, Float64[1:10...], Float64[1:10...], Float64[1:10...])
        @test A == map(x->x*x*x, Float64[1:10...])
        @test A === B
    end
end

function testmap_equivalence(mapf, f, c...)
    x1 = mapf(f,c...)
    x2 = map(f,c...)

    if Base.IteratorSize isa Base.HasShape
        @test size(x1) == size(x2)
    else
        @test length(x1) == length(x2)
    end

    @test eltype(x1) == eltype(x2)

    for (v1,v2) in zip(x1,x2)
        @test v1==v2
    end
end

function run_map_equivalence_tests(mapf)
    testmap_equivalence(mapf, identity, (1,2,3,4))
    testmap_equivalence(mapf, (x,y,z)->x+y+z, 1,2,3)
    testmap_equivalence(mapf, x->x ? false : true, BitMatrix(undef, 10,10))
    testmap_equivalence(mapf, x->"foobar", BitMatrix(undef, 10,10))
    testmap_equivalence(mapf, (x,y,z)->string(x,y,z), BitVector(undef, 10), fill(1.0, 10), "1234567890")
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
function implement(line::LineNumberNode,
                   constant_lookup :: Bool,
                   trait::Type{<:Trait},
                   @nospecialize(impls),
                   sup_impls   :: Vector{Expr},
                   type_params :: AbstractArray,
                   freshvars :: AbstractArray)
    mnames    = fieldnames(CanonicalTraits.instance(trait)) |> collect
    not_impls = mnames |> Set
    methods = @when Expr(:block, stmts...) = impls begin
        stmts
    @otherwise
        [impls]
    end
    defs = Pair{Symbol, Any}[]

    for i in 1:length(methods)
        push!(
            defs,
            @when (name, impl) = extract_impl(methods[i]) begin
                delete!(not_impls, name)
                name => impl
            @otherwise
                :_ => methods[i]
            end
        )
    end
    # generate default methods not implemented by user
    for mname in not_impls
        maker = default_method_maker(trait, Val(mname))
        bd    = :($maker($(type_params...)))
        push!(defs, mname => bd)
    end

    name_maps = Dict{Symbol, Symbol}()
    basename  = gensym(string(trait))

    cg_defs = map(defs) do (a, b)
        if a === :_
            b
        else
            sym = Symbol(basename, "#", a)
            name_maps[a] = sym
            :($sym = let; $b end)
        end
    end

    argtypes = [:(::Type{$tp}) for tp in type_params]
    if constant_lookup
        @q begin
            $line
            Base.@generated $line function (::Type{$trait})($(argtypes...)) where {$(freshvars...)}
                $line
                $check_inheritance($trait, $(type_params...))
                $(sup_impls...)
                $(cg_defs...)
                $(instance(trait))($([name_maps[mname] for mname in mnames]...))
            end
        end
    else
        @q begin
            $line
            function (::Type{$trait})($(argtypes...)) where {$(freshvars...)}
                $line
                $check_inheritance($trait, $(type_params...))
                $(sup_impls...)
                $(cg_defs...)
                $(instance(trait))($([name_maps[mname] for mname in mnames]...))
            end
        end
    end
end

function implement(source::LineNumberNode, constant_lookup, @nospecialize(sig), @nospecialize(block), mod::Module)
    freshvars = Any[]
    sups, sig = extract_trait_mk(sig)
    @when :($trait where {$(freshvars_...)}) = sig begin
        sig = trait
        freshvars = freshvars_
    end
    @when :($trait{$(type_params...)}) = sig begin
        implement(source, constant_lookup, mod.eval(trait), block, sups, type_params, freshvars)
    @otherwise
        error("Instance should be in form of 'Trait{A, B, C, D}' instead of '$sig'.")

    end
end

"""
If the implementation is pure static, try `@implement!` which
guarantees the constant-time instance lookup.
"""
macro implement(@nospecialize(sig), @nospecialize(block))
    implement(__source__, false, sig, block, __module__) |> esc
end

macro implement(@nospecialize(sig))
    implement(__source__, false, sig, Expr(:block), __module__) |> esc
end

"""
Implementing traits with a guarantee of constant-time instance lookup.
It's made by generated function, but once you call any methods who
belong to not yet implemented instances, you can never implement
it in this runtime.
"""
macro implement!(@nospecialize(sig), @nospecialize(block))
    implement(__source__, true, sig, block, __module__) |> esc
end

macro implement!(@nospecialize(sig))
    implement(__source__, true, sig, Expr(:block), __module__) |> esc
end


text:
# Test that an example fixed point problem converges to the exact solution. Based on the pyblp code.
@testset "Iteration converges $it" for it in (SimpleFixedPointIteration(), SQUAREMIteration())
    function contraction!(x_out::AbstractVector, x::AbstractVector)
        c1 = [10, 12]
        c2 = [3, 5]

        @. x_out = sqrt(c1 / (x + c2))
        return x_out
    end

    initial = ones(2)
    result = fixed_point_iteration(it, initial, contraction!)

    @test result.status == ITERATION_CONVERGED
    @test result.delta ≈ [1.4920333, 1.37228132]
end

# Test that the solution to the fixed point problem from Hasselblad (1969) is reasonably close to the exact
# solution. This same problem is used in an original SQUAREM unit test.
@testset "SQUAREM $scheme Hasselblad" for scheme in (SQUAREM1, SQUAREM2, SQUAREM3)
    function contraction!(x_out::AbstractVector, x::AbstractVector)
        y = [162, 267, 271, 185, 111, 61, 27, 8, 3, 1]
        i = 0:length(y)-1
        z = @. (x[1] * exp(-x[2]) * x[2]^i) / (x[1] * exp(-x[2]) * x[2]^i + (1 - x[1]) * exp(-x[3]) * x[3]^i)

        x_out[:] = @. [
            $sum(y * z) / $sum(y),
            $sum(y * i * z) / $sum(y .* z),
            $sum(y * i * (1-z)) / $sum(y * (1-z))
        ]
        return x_out
    end

    it = SQUAREMIteration(scheme=scheme)

    initial = [0.2, 2.5, 1.5]

    result = fixed_point_iteration(it, initial, contraction!)
    @test result.status == ITERATION_CONVERGED
    @test result.delta ≈ [0.6401146029910, 2.6634043566619, 1.2560951012662]
end


text:
@testset "Analytic" begin
    seed!(42)
    D = 10
    N = 20
    b = 5
    x = [rand(D) for i in 1:D]
    y = rand(N)
    i = Analytic()
    @test Analytic(ϵ = 0.0001f0) isa Analytic{Float32}
    @test repr(i) == "Analytic Inference"
    xview = view(x, :)
    yview = view(y, :)
    i = AGP.init_inference(i, N, xview, yview)

    @test AGP.xview(i) == view(x, :)
    @test AGP.yview(i) == view(y, :)

    @test AGP.nMinibatch(i) == N
    @test AGP.getρ(i) == 1.0
    @test AGP.isStochastic(i) == false
    @test AGP.MBIndices(i) == 1:N
end


text:
# Store network information
# TODO: echenolize model

struct MetNet

    ## LP (original)
    S::Matrix{Float64}
    b::Vector{Float64}
    lb::Vector{Float64}
    ub::Vector{Float64}
    c::Vector{Float64}
    rxns::Vector{String}
    mets::Vector{String}

end

MetNet(;S, b, lb, ub, c, rxns, mets) = MetNet(S, b, lb, ub, c, rxns, mets)

## -------------------------------------------------------------------
# Utils
_IDER_TYPE = Union{String, Symbol, Int}
_IDER_TYPE_STR = Union{String, Symbol}
_IDER_TYPE_INT = Int

rxnindex(net::MetNet, rxn::_IDER_TYPE_STR) = findfirst(isequal(string(rxn)), net.rxns)
rxnindex(net::MetNet, rxn::_IDER_TYPE_INT) = findfirst(isequal(rxn), eachindex(net.rxns))
rxnindex(net::MetNet, rxns::Vector) = rxnindex.([net], rxns)
rxnindex(net::MetNet, rxns) = rxnindex(net, collect(rxns))

metindex(net::MetNet, met::_IDER_TYPE_STR) = findfirst(isequal(string(met)), net.mets)
metindex(net::MetNet, met::_IDER_TYPE_INT) = findfirst(isequal(met), eachindex(net.mets))
metindex(net::MetNet, mets::Vector) = metindex.([net], mets)
metindex(net::MetNet, mets) = metindex(net, collect(mets))

rxns(net::MetNet, rxns) = net.rxns[rxnindex(net, rxns)]

## -------------------------------------------------------------------
function reducebox!(net::MetNet)
    lb, ub = fva(net)
    net.lb .= lb
    net.ub .= ub
    return net
end
_reducebox!(net, reduce) = reduce ? reducebox!(net) : net

function fix!(net::MetNet, rxns, val; reducebox = false)
    reducebox && reducebox!(net)
    idxs = rxnindex(net, rxns)
    idxs = (idxs isa Vector) ? idxs : [idxs]
    net.lb[idxs] .= val
    net.ub[idxs] .= val
    net
end

function fixxing(f::Function, net::MetNet, rxns, val; tol = 0.0)

    idxs = rxnindex(net, rxns)
    bk_lb = net.lb[idxs]
    bk_ub = net.ub[idxs]
    net.lb[idxs] = (val .- tol)
    net.ub[idxs] = (val .+ tol)
    try; return f()
    finally;
        net.lb[idxs] = bk_lb
        net.ub[idxs] = bk_ub
    end
end

rxnid(net::MetNet, rxn) = net.rxn[rxnindex(net, rxn)]
lb(net::MetNet, rxn) = net.lb[rxnindex(net, rxn)]
lb!(net::MetNet, rxn, val; reducebox = false) =
    (net.lb[rxnindex(net, rxn)] = val; _reducebox!(net, reducebox))
ub(net::MetNet, rxn) = net.ub[rxnindex(net, rxn)]
ub!(net::MetNet, rxn, val; reducebox = false) =
    (net.ub[rxnindex(net, rxn)] = val; _reducebox!(net, reducebox))
bounds!(net::MetNet, rxn, lb, ub; reducebox = false) =
    (idx = rxnindex(net, rxn); net.lb[idx] = lb; net.ub[idx] = ub; _reducebox!(net, reducebox))
bounds(net::MetNet, rxn) =
    (idx = rxnindex(net, rxn); (net.lb[idx],  net.ub[idx]))

## -------------------------------------------------------------------
_get_elstr(v, i) = i > lastindex(v) ? "" : string(v[i])
function pretty(;kwargs...)
    kwargs = collect(kwargs)
    ks = first.(kwargs)
    vs = last.(kwargs)
    # maximum length
    padl = map(zip(ks, vs)) do (k, v)
        maxvl = maximum(length.(string.(v)))
        max(maxvl, length(string(k)))
    end

    tab = "   "
    println(join(rpad.(first.(kwargs), padl), tab))
    println(join( ["-"].^padl , tab))
    maxl = maximum(length.(last.(kwargs)))
    for i in 1:maxl
        println(join(rpad.(_get_elstr.(vs, i), padl), tab))
    end

end

## -------------------------------------------------------------------
_round(n) = round(n; sigdigits = 3)
function prettyrxn(net::MetNet, v::Vector, rxns = eachindex(net.rxns))
    rxns = rxnindex(net, rxns)
    pretty(;
        n = collect(eachindex(net.rxns))[rxns],
        rxn = net.rxns[rxns],
        v = _round.(v)[rxns],
        lb = _round.(net.lb)[rxns],
        ub = _round.(net.ub)[rxns],
        eq = rxn_str.([net], net.rxns)[rxns]
    )
end

## -------------------------------------------------------------------
# LP
fba(net::MetNet; kwargs...) = fba(net.S, net.b, net.lb, net.ub, net.c; kwargs...)
fba(net::MetNet, obj, sense = MAX_SENSE; kwargs...) =
    fba(net.S, net.b, net.lb, net.ub, net.c, rxnindex(net, obj), sense; kwargs...)
fva(net::MetNet, rxns = eachindex(net.rxns); kwargs...) =
    fva(net.S, net.b, net.lb, net.ub, net.c, rxnindex(net, rxns); kwargs...)
Δv(net, rxns = eachindex(net.rxns)) =
    ((lb, ub) = fva(net, rxnindex(net, rxns)); ub - lb)
U(net, rxns = eachindex(net.rxns)) =
    ((lb, ub) = fva(net, rxnindex(net, rxns)); ub)
L(net, rxns = eachindex(net.rxns)) =
    ((lb, ub) = fva(net, rxnindex(net, rxns)); lb)

## -------------------------------------------------------------------
Base.hash(net::MetNet) = hash((:MetNet, net.S, net.b, net.lb, net.ub))
Base.size(net::MetNet) = size(net.S)

## -------------------------------------------------------------------
function rxn_str(net::MetNet, rxn)
    ri = rxnindex(net, rxn)
    rS = net.S[:, ri]

    reacts = String[]
    prods = String[]
    for (mi, ms) in enumerate(rS)
        mid = net.mets[mi]
        ms = round(ms; sigdigits = 3)
        (ms < 0.0) && push!(reacts, string("(", ms, ")", mid))
        (ms > 0.0) && push!(prods, string("(", ms, ")", mid))
    end

    reacts_str = join(reacts, " + ")
    prods_str = join(prods, " + ")

    lb, ub = net.lb[ri], net.ub[ri]
    arrow = (lb == ub == 0) ?
        ">-<" : (lb < 0.0 && ub > 0.0) ?
            "<-->" : (lb >= 0.0 && ub > 0.0) ?
                "-->" : (lb < 0.0 && ub <= 0.0) ?
                    "<--" : "error"

    string(reacts_str, " ", arrow, " ", prods_str)

end

## -------------------------------------------------------------------
function Base.show(io::IO, net::MetNet)

    println(io, "MetNet: ", size(net))

    M, N = size(net)

    ri_strs = String["ri"; string.(1:N)]
    rxnid_strs = String["rxnid"]
    bounds_strs = String["bounds"]
    rxni_strs = String["rxn str"]
    for ri in 1:N
        lb, ub = net.lb[ri], net.ub[ri]
        lb = round(lb; sigdigits = 3)
        ub = round(ub; sigdigits = 3)
        bounds_str = string("[", lb, ", " ,ub, "]")
        push!(bounds_strs, bounds_str)

        rid = net.rxns[ri]
        rxnid_str = string(rid)
        push!(rxnid_strs, rxnid_str)

        rxni_str = rxn_str(net, ri)
        push!(rxni_strs, rxni_str)

    end

    ri_pad = maximum(length.(ri_strs))
    rxnid_pad = maximum(length.(rxnid_strs))
    rxni_pad = maximum(length.(rxni_strs))
    bounds_pad = maximum(length.(bounds_strs))
    tab = "   "
    for ri in eachindex(rxnid_strs)
        println(io,
            join([
                rpad(ri_strs[ri], ri_pad),
                rpad(rxnid_strs[ri], rxnid_pad),
                rpad(bounds_strs[ri], bounds_pad),
                rpad(rxni_strs[ri], rxni_pad)
            ], tab)
        )
    end
    nothing
end

Base.show(net::MetNet) = show(stdout, net)

## -------------------------------------------------------------------
function is_feasible(net::MetNet, freeis, v; tol::Vector = Δv(net, freeis) .* 0.01)
    fullv = fixxing(net, freeis, v; tol) do
        fba(net)
    end
    any(fullv .!= 0.0)
end

## -------------------------------------------------------------------
function BoxGrid(net::MetNet, frees::Vector, δs::Vector)
    idxs = rxnindex.([net], frees)
    dimdat = []
    for (idx, free, δ) in zip(idxs, frees, δs)
        isnothing(idx) && error(free, " not found")
        id = net.rxns[idx]
        lb, ub = net.lb[idx], net.ub[idx]
        push!(dimdat, (Symbol(id), lb, ub, δ))
    end
    BoxGrid(dimdat)
end
BoxGrid(net::MetNet, frees::Vector, δ::Int) = BoxGrid(net, frees, fill(δ, length(frees)))
BoxGrid(net::MetNet, δ::Int) = BoxGrid(net, net.rxns, fill(δ, length(net.rxns)))


text:
#=
setindex.jl

Provide setindex!() capabilities like A[i,j] = s for LinearMapAA objects.

The setindex! function here is (necessarily) quite inefficient.
It is provided solely so that a LinearMapAA conforms to the AbstractMatrix
requirement that a setindex! method exists.
Use of this function is strongly discouraged, other than for testing.

2018-01-19, Jeff Fessler, University of Michigan
=#

using Test: @test
using LinearMaps #: LinearMaps.FunctionMap

Indexer = AbstractVector{Int}

"""
`A[i,j] = X`

Mathematically, if B = copy(A) and then we say `B[i,j] = s`, for scalar `s`,
then `B = A + (s - A[i,j]) e_i e_j'`
where `e_i` and `e_j` are the appropriate unit vectors.
Using this basic math, we can perform `B*x` using `A*x` and `A[i,j]`.

This idea generalizes to `B[ii,jj] = X` where `ii` and `jj` are vectors.

This method works because LinearMapAA is a *mutable* struct.
"""
function Base.setindex!(A::LinearMapAA, X::AbstractMatrix,
		ii::Indexer, jj::Indexer)
	# todo: handle WrappedMap differently
	L = A._lmap
	Aij = A[ii,jj] # must extract this outside of forw() to avoid recursion!
	Aij = reshape(Aij, size(X)) # needed for [i,:] case
	if true
		forw = (x) ->
			begin
				tmp = L * x
				tmp[ii] += (X - Aij) * x[jj]
				return tmp
			end
	end

	has_adj = !(typeof(L) <: LinearMaps.FunctionMap) || (L.fc !== nothing)
	if has_adj
		back = (y) ->
			begin
				tmp = L' * y
				tmp[jj] += (X - Aij)' * y[ii]
				return tmp
			end
	end

#=
	# for future reference, some of this could be implemented
	A.issymmetric = false # could check if i==j
	A.ishermitian = false # could check if i==j and v=v'
	A.posdef = false # could check if i==j and v > A[i,j]
=#

	if has_adj
		A._lmap = LinearMap(forw, back, size(L)...)
	else
		A._lmap = LinearMap(forw, size(L)...)
	end

	nothing
end

# [i,j] = s
Base.setindex!(A::LinearMapAA, s::Number, i::Int, j::Int) =
	setindex!(A, fill(s,1,1), [i], [j])

# [ii,jj] = X
Base.setindex!(A::LinearMapAA, X::AbstractMatrix,
	ii::AbstractVector{Bool}, jj::AbstractVector{Bool}) =
	setindex!(A, X, findall(ii), findall(jj))

# [:,jj] = X
Base.setindex!(A::LinearMapAA, X::AbstractMatrix, ::Colon, jj::AbstractVector{Bool}) =
	setindex!(A, X, :, findall(jj))
Base.setindex!(A::LinearMapAA, X::AbstractMatrix, ::Colon, jj::Indexer) =
	setindex!(A, X, 1:size(A,1), jj)

# [ii,:] = X
Base.setindex!(A::LinearMapAA, X::AbstractMatrix, ii::AbstractVector{Bool}, ::Colon) =
	setindex!(A, X, findall(ii), :)
Base.setindex!(A::LinearMapAA, X::AbstractMatrix, ii::Indexer, ::Colon) =
	setindex!(A, X, ii, 1:size(A,2))

# [:,j] = v
Base.setindex!(A::LinearMapAA, v::AbstractVector, ::Colon, j::Int) =
	setindex!(A, reshape(v,:,1), :, [j])

# [ii,j] = v
Base.setindex!(A::LinearMapAA, v::AbstractVector, ii::Indexer, j::Int) =
	setindex!(A, reshape(v,:,1), ii, [j])

# [i,:] = v
Base.setindex!(A::LinearMapAA, v::AbstractVector, i::Int, ::Colon) =
	setindex!(A, reshape(v,1,:), [i], :)

# [i,jj] = v
Base.setindex!(A::LinearMapAA, v::AbstractVector, i::Int, jj::Indexer) =
	setindex!(A, reshape(v,1,:), [i], jj)

# [ii,jj] = s
Base.setindex!(A::LinearMapAA, s::Number, ii::AbstractVector{Bool}, jj::Indexer) =
	setindex!(A, s, findall(ii), jj)
Base.setindex!(A::LinearMapAA, s::Number, ii::Indexer, jj::AbstractVector{Bool}) =
	setindex!(A, s, ii, findall(jj))
Base.setindex!(A::LinearMapAA, s::Number,
	ii::AbstractVector{Bool}, jj::AbstractVector{Bool}) =
	setindex!(A, s, findall(ii), findall(jj))
Base.setindex!(A::LinearMapAA, s::Number, ii::Indexer, jj::Indexer) =
	setindex!(A, fill(s,length(ii),length(jj)), ii, jj)

# [:,j] = s
Base.setindex!(A::LinearMapAA, s::Number, ::Colon, j::Int) =
	setindex!(A, s, 1:size(A,1), [j])

# [ii,:] = s
Base.setindex!(A::LinearMapAA, s::Number, ii::AbstractVector{Bool}, ::Colon) =
	setindex!(A, s, findall(ii), :)
Base.setindex!(A::LinearMapAA, s::Number, ii::Indexer, ::Colon) =
	setindex!(A, s, ii, 1:size(A,2))

# [ii,j] = s
Base.setindex!(A::LinearMapAA, s::Number, ii::AbstractVector{Bool}, j::Int) =
	setindex!(A, s, findall(ii), [j])
Base.setindex!(A::LinearMapAA, s::Number, ii::Indexer, j::Int) =
	setindex!(A, fill(s, length(ii), 1), ii, [j])

# [i,:] = s
Base.setindex!(A::LinearMapAA, s::Number, i::Int, ::Colon) =
	setindex!(A, fill(s, 1, size(A,2)), [i], :)

# [i,jj] = s
Base.setindex!(A::LinearMapAA, s::Number, i::Int, jj::AbstractVector{Bool}) =
	setindex!(A, s, i, findall(jj))
Base.setindex!(A::LinearMapAA, s::Number, i::Int, jj::Indexer) =
	setindex!(A, fill(s, 1, length(jj)), [i], jj)

# [:,jj] = s
Base.setindex!(A::LinearMapAA, s::Number, ::Colon, jj::AbstractVector{Bool}) =
	setindex!(A, s, :, findall(jj))
Base.setindex!(A::LinearMapAA, s::Number, ::Colon, jj::Indexer) =
	setindex!(A, s, 1:size(A,1), jj)

# [kk]
#= too much work, so unsupported
Base.setindex!(A::LinearMapAA, v::AbstractVector,
		kk::Indexer) =
=#

# [:] = v
Base.setindex!(A::LinearMapAA, v::AbstractVector, ::Colon) =
	setindex!(A, reshape(v,size(A)), :, :)

# [:] = s
Base.setindex!(A::LinearMapAA, s::Number, ::Colon) =
	begin
		(M,N) = size(A)
		forw = x -> fill(s, M) * (ones(1,N) * x)[1]
		back = y -> fill(conj(s), N) * (ones(1, M) * y)[1]
		A._lmap = LinearMap(forw, back, M, N)
	end

# [:,:] = s
Base.setindex!(A::LinearMapAA, s::Number, ::Colon, ::Colon) =
	setindex!(A, s, :)

# [:,:] = X
Base.setindex!(A::LinearMapAA, X::AbstractMatrix, ::Colon, ::Colon) =
	begin
		A._lmap = LinearMap(X)
	end

# [k] = s
"""
`setindex!(A::LinearMapAA, s::Number, k::Int)`

Provide the single index version to meet the `AbstractArray` spec:
https://docs.julialang.org/en/latest/manual/interfaces/#Indexing-1
"""
Base.setindex!(A::LinearMapAA, s::Number, k::Int) =
	setindex!(A, s, Tuple(CartesianIndices(size(A))[k])...) # any better way?

#=
"""
`setindex!(A::LinearMapAA, v::Number, k::Int)`

Provide the single index version to meet the `AbstractArray` spec:
https://docs.julialang.org/en/latest/manual/interfaces/#Indexing-1
"""
function Base.setindex!(A::LinearMapAA, v::Number, k::Int)
	c = CartesianIndices(size(A))[k] # is there a more elegant way?
	setindex!(A, v::Number, c[1], c[2])
end
=#

"""
`LinearMapAA_test_setindex(A::LinearMapAA)`
"""
function LinearMapAA_test_setindex(A::LinearMapAA)

    @test all(size(A) .>= (4,4)) # required by tests

	tf1 = [false; trues(size(A,1)-1)]
	tf2 = [false; trues(size(A,2)-2); false]
	ii1 = (3, 2:4, [2,4], :, tf1)
	ii2 = (2, 3:4, [1,4], :, tf2)

	# test all [?,?] combinations with array "X"
	for i2 in ii2
		for i1 in ii1
			B = copy(A)
			X = 2 .+ A[i1,i2].^2 # values must differ from A[ii,jj]
			B[i1,i2] = X
			Am = Matrix(A)
			Bm = Matrix(B)
			Am[i1,i2] = X
			@test isapprox(Am, Bm)
		end
	end

	# test all [?,?] combinations with scalar "s"
	for i2 in ii2
		for i1 in ii1
			B = copy(A)
			s = maximum(abs.(A[:])) + 2
			B[i1,i2] = s
			Am = Matrix(A)
			Bm = Matrix(B)
			if (i1 == Colon() || ndims(i1) > 0) || (i2 == Colon() || ndims(i2) > 0)
				Am[i1,i2] .= s
			else
				Am[i1,i2] = s
			end
			@test isapprox(Am, Bm)
		end
	end

	# others not supported for now
	set1 = (3, ) # [3], [2,4], 2:4, (1:length(A)) .== 2), end # todo
	for s1 in set1
		B = copy(A)
		X = 2 .+ A[s1].^2 # values must differ from A[ii,jj]
		B[s1] = X
		Am = Matrix(A)
		Bm = Matrix(B)
		Am[s1] = X
		@test isapprox(Am, Bm)
	end

	# insanity below here

	# A[:] = s
	B = copy(A)
	B[:] = 5
	Am = Matrix(A)
	Am[:] .= 5
	Bm = Matrix(B)
	@test Bm == Am

	# A[:] = v
	B = copy(A)
	v = 1:length(A)
	B[:] = v
	Am = Matrix(A)
	Am[:] .= v
	Bm = Matrix(B)
	@test Bm == Am

	# A[:,:]
	B = copy(A)
	B[:,:] = 6
	Am = Matrix(A)
	Am[:,:] .= 6
	Bm = Matrix(B)
	@test Bm == Am

	true
end


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

@doc raw"""SubjectAccessReviewStatus

    IoK8sApiAuthorizationV1SubjectAccessReviewStatus(;
        allowed=nothing,
        denied=nothing,
        evaluationError=nothing,
        reason=nothing,
    )

    - allowed::Bool : Allowed is required. True if the action would be allowed, false otherwise.
    - denied::Bool : Denied is optional. True if the action would be denied, otherwise false. If both allowed is false and denied is false, then the authorizer has no opinion on whether to authorize the action. Denied may not be true if Allowed is true.
    - evaluationError::String : EvaluationError is an indication that some error occurred during the authorization check. It is entirely possible to get an error and be able to continue determine authorization status in spite of it. For instance, RBAC can be missing a role, but enough roles are still present and bound to reason about the request.
    - reason::String : Reason is optional.  It indicates why a request was allowed or denied.
"""
mutable struct IoK8sApiAuthorizationV1SubjectAccessReviewStatus <: SwaggerModel
    allowed::Any # spec type: Union{ Nothing, Bool } # spec name: allowed
    denied::Any # spec type: Union{ Nothing, Bool } # spec name: denied
    evaluationError::Any # spec type: Union{ Nothing, String } # spec name: evaluationError
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason

    function IoK8sApiAuthorizationV1SubjectAccessReviewStatus(;allowed=nothing, denied=nothing, evaluationError=nothing, reason=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewStatus, Symbol("allowed"), allowed)
        setfield!(o, Symbol("allowed"), allowed)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewStatus, Symbol("denied"), denied)
        setfield!(o, Symbol("denied"), denied)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewStatus, Symbol("evaluationError"), evaluationError)
        setfield!(o, Symbol("evaluationError"), evaluationError)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewStatus, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        o
    end
end # type IoK8sApiAuthorizationV1SubjectAccessReviewStatus

const _property_map_IoK8sApiAuthorizationV1SubjectAccessReviewStatus = Dict{Symbol,Symbol}(Symbol("allowed")=>Symbol("allowed"), Symbol("denied")=>Symbol("denied"), Symbol("evaluationError")=>Symbol("evaluationError"), Symbol("reason")=>Symbol("reason"))
const _property_types_IoK8sApiAuthorizationV1SubjectAccessReviewStatus = Dict{Symbol,String}(Symbol("allowed")=>"Bool", Symbol("denied")=>"Bool", Symbol("evaluationError")=>"String", Symbol("reason")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1SubjectAccessReviewStatus }) = collect(keys(_property_map_IoK8sApiAuthorizationV1SubjectAccessReviewStatus))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1SubjectAccessReviewStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthorizationV1SubjectAccessReviewStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1SubjectAccessReviewStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1SubjectAccessReviewStatus[property_name]

function check_required(o::IoK8sApiAuthorizationV1SubjectAccessReviewStatus)
    (getproperty(o, Symbol("allowed")) === nothing) && (return false)
    true
end

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


text:
module PrescriptionUtils

using Compat

import GLMNet
import MLDataUtils
import NearestNeighbors
import PyCall
import ScikitLearn
import StatsBase

import Compat.Statistics: mean

export
  evaluatebaseline,
  evaluateoracle,
  evaluateprescriptions,
  getcounterfactuals,
  getoutcomes,
  makeprescriptions,
  predictoutcomes

# @sk_import creates a global pointer that doesn't work with precompilation
# copy directly: https://github.com/cstjean/ScikitLearn.jl/issues/50
const ensemble = PyCall.PyNULL()

function __init__()
  ScikitLearn.Skcore.import_sklearn()
  copy!(ensemble, PyCall.pyimport("sklearn.ensemble"))
end

include("data.jl")
include("prescriptions.jl")
include("evaluation.jl")

include("knn.jl")
include("lasso.jl")
include("randomforest.jl")

end # module


text:
# Syntax highlighting and other enhancements for the Julia REPL with Monokai theme by default
using OhMyREPL
@async begin
    # reinstall keybindings to work around https://github.com/KristofferC/OhMyREPL.jl/issues/166
    sleep(1)
    OhMyREPL.Prompt.insert_keybindings()
end
colorscheme!("Monokai24bit")
# Load Debugger with Monokai theme by default
using Highlights
using Debugger
Debugger.set_theme(Highlights.Themes.MonokaiTheme)
# Automatically update function definitions in a running Julia session
# Use includet() instead include() to track the files
atreplinit() do repl
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
    catch
    end
end


text:
#=
Set up the infrastructure, and trigger a quick failure if appropriate.
=#
if Sys.isapple()
    using Pkg
    Pkg.add("Conda")
    using Conda
    Conda.add("pyqt")
end
Pkg.add("Plots")
using Plots
pyplot()
p=scatter(rand(10),rand(10))
png(p,"/tmp/fubar.png")


text:
export @esc, isexpr, isline, rmlines, unblock, block, inexpr, namify, isdef,
  longdef, shortdef, @expand, makeif, prettify

assoc!(d, k, v) = (d[k] = v; d)

macro esc(xs...)
  :($([:($x = esc($x)) for x in map(esc, xs)]...);)
end

"""
    isexpr(x, ts...)

Convenient way to test the type of a Julia expression.
Expression heads and types are supported, so for example
you can call

    isexpr(expr, String, :string)

to pick up on all string-like expressions.
"""
isexpr(x::Expr) = true
isexpr(x) = false
isexpr(x::Expr, ts...) = x.head in ts
isexpr(x, ts...) = any(T->isa(T, Type) && isa(x, T), ts)

isline(ex) = isexpr(ex, :line) || isa(ex, LineNumberNode)

"""
    rmlines(x)

Remove the line nodes from a block or array of expressions.
"""
rmlines(x) = x
rmlines(x::Expr) = Expr(x.head, filter(x->!isline(x), x.args)...)

striplines(ex) = prewalk(rmlines, ex)

"""
    unblock(expr)

Remove outer `begin` blocks from an expression, if the block is
redundant (i.e. contains only a single expression).
"""
function unblock(ex)
  isexpr(ex, :block) || return ex
  exs = rmlines(ex).args
  length(exs) == 1 || return ex
  return unblock(exs[1])
end

block(ex) = isexpr(ex, :block) ? ex : :($ex;)

"""
An easy way to get pull the (function/type) name out of
expressions like `foo{T}` or `Bar{T} <: Vector{T}`.
"""
namify(s::Symbol) = s
namify(ex::Expr) = namify(ex.args[1])

Base.macroexpand(m::Module, ex) =
  eval(m, :(macroexpand($(Expr(:quote, ex)))))

walk(x, inner, outer) = outer(x)
walk(x::Expr, inner, outer) = outer(Expr(x.head, map(inner, x.args)...))

postwalk(f, x) = walk(x, x -> postwalk(f, x), f)
prewalk(f, x)  = walk(f(x), x -> prewalk(f, x), identity)

replace(ex, s, s′) = prewalk(x -> x == s ? s′ : x, ex)

function inexpr(ex, x)
  result = false
  MacroTools.postwalk(ex) do y
    y == x && (result = true)
  end
  return result
end

global const animals = split(readstring(joinpath(dirname(@__FILE__), "..", "animals.txt")))

isgensym(s::Symbol) = contains(string(s), "#")
isgensym(s) = false

function alias_gensyms(ex)
  syms = Dict{Symbol, Symbol}()
  s(x) = get!(syms, x, lowercase(rand(filter(s->!(s in values(syms)), animals))))
  prewalk(ex) do x
    isgensym(x) ? s(x) : x
  end
end

"""
More convenient macro expansion, e.g.

    @expand @time foo()
"""
macro expand(ex)
  ex = postwalk(ex) do ex
    isexpr(ex, :$) ? Expr(:quote, ex) : ex
  end
  :(alias_gensyms(macroexpand($(Expr(:quote, ex)))))
end

"Test for function definition expressions."
isdef(ex) = ismatch(or_(:(function _(__) _ end),
                        :(f_(__) = _)),
                    ex)

function longdef(ex)
  prewalk(ex) do ex
    @match ex begin
      (f_(args__) = body_) => :(function $f($(args...)) $body end)
      ((args__,) -> body_) => :(function ($(args...),) $body end)
      (arg_ -> body_) => :(function ($arg,) $body end)
      _ => ex
    end
  end
end

function shortdef(ex)
  prewalk(ex) do ex
    @match ex begin
      function f_(args__) body_ end => :($f($(args...)) = $body)
      function (args__,) body_ end => :(($(args...),) -> $body)
      (arg_ -> body_) => :(($arg,) -> $body)
      _ => ex
    end
  end
end

function flatten1(ex)
  isexpr(ex, :block) || return ex
  ex′ = quote end
  for x in ex.args
    isexpr(x, :block) ? append!(ex′.args, x.args) : push!(ex′.args, x)
  end
  return ex′
end

flatten(ex) = postwalk(flatten1, ex)

function makeif(clauses, els = nothing)
  foldr((c, ex)->:($(c[1]) ? $(c[2]) : $ex), els, clauses)
end

unresolve1(x) = x
unresolve1(f::Function) = methods(f).mt.name

unresolve(ex::Expr) = prewalk(unresolve1, ex)

function resyntax(ex)
  prewalk(ex) do x
    @match x begin
      setfield!(x_, :f_, x_.f_ + v_) => :($x.$f += $v)
      setfield!(x_, :f_, v_) => :($x.$f = $v)
      getindex(x_, i__) => :($x[$(i...)])
      tuple(xs__) => :($(xs...),)
      ctranspose(x_) => :($x')
      transpose(x_) => :($x.')
      _ => x
    end
  end
end

prettify(ex) = ex |> flatten |> unresolve |> resyntax |> alias_gensyms |> striplines


text:
using Test, MyTestRepo

x = [0.0, 2.0, 4.0]
y = [0.0, 3.0, 5.0]
α = lsq_fit(x, y)
@test isapprox(α, [1/6, 1.25])

plot_fit(x, y, "solution.png")


text:
#=

This file is the entry-point for the RadiativeTransfer module.

It includes this module's source files and exports the relevant keywords.

=#

module RadiativeTransfer
using Pkg.Artifacts
using LinearAlgebra
using Distributions
using CUDA
using Parameters
using DocStringExtensions

# Export Architecture functions
export CPU, GPU, default_architecture, array_type

# Export the artifact convenience function
export artifact

# GPU/CPU Architecture (from Oceanigans)
include("Architectures.jl")
using .Architectures

# Artifacts
include("Artifacts/artifact_helper.jl")

# Absorption Cross Section module:
include("Absorption/Absorption.jl")

# Mie Phase Function module:
include("Scattering/Scattering.jl")

# vSmartMOM module:
include("vSmartMOM/vSmartMOM.jl")
using .vSmartMOM

# SolarModel module:
include("SolarModel/SolarModel.jl")

# Export some vSmartMOM functions
export default_parameters, parameters_from_yaml, model_from_parameters, rt_run

using .Architectures
using .Absorption

# Perform some GPU setup when the module is loaded
function __init__()
    @hascuda begin
        @info "CUDA-enabled GPU(s) detected:"
        for (gpu, dev) in enumerate(CUDA.devices())
            @info "$dev: $(CUDA.name(dev))"
        end
	CUDA.allowscalar(false)
    end
end

end


text:
using Mixins
using Test

@testset "Mixins.jl" begin
  include("Mixins.jl")
end


text:
export rcd, ∥ₛ, ∥
"""
    `rcd(x, θ)`

Random Conditional Distribution of `x` given `Θ`.

Random conditional distributions provide a mechanism to condition distributional properties.
Given two random variables ``X`` and ``Θ``, the random conditional distribution of ``X`` given ``Θ`` -- which we denote ``X ∥ Θ`` -- is a a random distribution:
a random variable whose `elemtype` is a random variables.
In particular, each realization of ``X ∥ Θ`` is the random variable ``X`` conditioned on ``Θ = θ`` where ``θ ∼ Θ`` is a realization of ``Θ``:

The random conditional distribution (rcd) of a random variable ``X: Ω → τ_1`` given ``Θ: Ω → τ_2`` is a random variable ``X ∥ Θ:  Ω → ( Ω → τ_1)``, defined as:

```
(x ∥ θ) = ω -> x | (Θ = Θ(ω))
```
"""

rcd(x, θ, eq = ==ₚ) =  ω -> cond(x, eq(θ, θ(ω)))
# rcd(x, θs::Tuple, eq = ==ₚ) = rcd(x, randtuple(θs), eq)

"`rcd`, x ∥ θ"
x ∥ θ = rcd(x, θ)

"`rcd` with soft equality"
# rcdₛ(x, θ) = rcd(x, θ, ==ₛ)

# "`rcd` with soft equality"
# x ∥ₛ θ = rcd(x, θ, ==ₛ)


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script collecting all functions used in the repository
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Functions for making raw data operational
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# initial cleaning of the data
function initial_cleaning(year::String, month::String)

    path = dir_dropbox * "rawdata/comext/" * "full" * year * month * ".dat"
    df = CSV.read(path, DataFrame)

    # column types
    # Notes:
    #   - treat PERIOD as numeric to allow sorting
    transform!(df, [:DECLARANT, :PARTNER, :PERIOD] .=> ByRow(Int64), renamecols=false)

    # additional PRODUCT_BEC5 column from 2017 onwards
    if year in string.(2001:2016) # why does shorthand not work?
        df.PRODUCT_BEC5 .= missing
    else
    end

    # column names different: PRODUCT_cpa2002, PRODUCT_cpa2008
    if year in string.(2001:2001)
        rename!(df, :PRODUCT_cpa2002 => :PRODUCT_CPA2002, :PRODUCT_cpa2008 => :PRODUCT_CPA2008)
    elseif year in string.(2002:2007)
        rename!(df, :PRODUCT_cpa2008 => :PRODUCT_CPA2008)
    end

    cols_string = ["DECLARANT_ISO", "PARTNER_ISO", "TRADE_TYPE", "PRODUCT_NC", "PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1",
                   "PRODUCT_BEC", "PRODUCT_BEC5", "PRODUCT_SECTION", "FLOW", "STAT_REGIME"]
    transform!(df,  cols_string .=> ByRow(string), renamecols=false)

    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :SUP_QUANTITY] .=> ByRow(Float64), renamecols=false)

    # formatting
    df[:, :TRADE_TYPE] .= ifelse.(df[:, :TRADE_TYPE] .== "I", "intra", "extra")
    df[:, :FLOW] .= ifelse.(df[:, :FLOW] .== "1", "imports", "exports")
    df[:, :SUPP_UNIT] .= ifelse.(ismissing.(df[:, :SUPP_UNIT]), missing, string.(df[:, :SUPP_UNIT])) # to make sure type is of Union{Missing, String}

    # missing data
    # Notes:
    #   - if value or quanity is reported as 0 use missing instead
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :SUP_QUANTITY] .=> ByRow(x->ifelse(x == 0.0, missing, x)), renamecols=false)

    # aggregate over STAT_REGIME, SUPP_UNIT
    # Notes:
    #   - only lose few rows ~1%
    #   - disregard SUPP_UNIT and SUP_QUANTITY for the time being
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "TRADE_TYPE", "PRODUCT_NC", "PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1",
                     "PRODUCT_BEC", "PRODUCT_BEC5", "PRODUCT_SECTION", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)

    return df
end

# function to create WORLD as a seperate partner
# Notes:
#   - need to create WORLD before EU, otherwise EU countries counted twice
function append_WORLD(df::DataFrame)

    # sum over PARTNER
    df_WORLD = subset(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x)))
    subset!(df_WORLD, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> x != 0))
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1",
                    "PRODUCT_BEC", "PRODUCT_BEC5", "PRODUCT_SECTION", "FLOW", "PERIOD"]
    gdf = groupby(df_WORLD, cols_grouping)
    df_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)

    df_WORLD.TRADE_TYPE .= "total"
    df_WORLD.PARTNER_ISO .= "WORLD"

    df = vcat(df, df_WORLD)

    return df
end

# function to create EU as a seperate partner
# Notes:
#   - with the vector "EU" one can specify the exact countries
function append_EU(df::DataFrame, EU::Vector{String})

    # subset
    df_EU = subset(df, :TRADE_TYPE => ByRow(x-> x == "intra"), :PARTNER_ISO => ByRow(x -> x in EU))
    subset!(df_EU, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x)))
    subset!(df_EU, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> x != 0))
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1",
                    "PRODUCT_BEC", "PRODUCT_BEC5", "PRODUCT_SECTION", "FLOW", "PERIOD"]
    gdf = groupby(df_EU, cols_grouping)
    df_EU = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)
    df_EU.TRADE_TYPE .= "intra"
    df_EU.PARTNER_ISO .= "EU"

    df = vcat(df, df_EU)

    return df
end

# function to import data for single country
function import_ctry(year::String, month::String, ctry::String)

    path = dir_dropbox * "rawdata/comext/" * ctry * "/" * "full" * year * month * "_" * ctry * ".csv"
    df = CSV.read(path, DataFrame)

    cols_string = ["DECLARANT_ISO", "PARTNER_ISO", "TRADE_TYPE", "PRODUCT_NC", "PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1",
                "PRODUCT_BEC", "PRODUCT_BEC5", "PRODUCT_SECTION", "FLOW"]
    transform!(df,  cols_string .=> ByRow(string), renamecols=false)

    # take out rows with missing values and quantities (already taken out when computing EU/WORLD)
    subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x)))
    subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> x != 0))
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(Float64), renamecols=false)

    return df
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Functions for computing statistics
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# computes pct_change between successive periods
# Notes:
#   - percentages?
#   - not every product imported/exported every period
#       + calculate price changes between periods when products are bought, i.e. price change between Feb2019 and March2019,
#         but could also be Feb2019 and Dec2020!
function pct_change(input::AbstractVector)
    [i == 1 ? missing : (input[i]-input[i-1])/input[i-1]*100 for i in eachindex(input)]
end

# computes pct_change between successive months
# Notes:
#   - functions computes UNIT_PRICE_CHANGE only for successive periods (months)
#       + notice the double ifelse statement in shorthand, find a way to use ifelse instead to make it easier to read!
function mom_change(period::AbstractVector, input::AbstractVector)
    [i == 1 ? missing : period[i]-period[i-1] != 1 ? missing : (input[i]-input[i-1])/input[i-1]*100 for i in eachindex(input)]
end

# moving average
# Notes:
#   - specify range with 'n'
movingaverage(input::AbstractArray, n::Int64) = [i < n ? missing : mean(input[i-n+1:i]) for i in eachindex(input)]

# computes standard deviation for rolling window
# requires PERIOD to be sorted when grouped
function rolling_std(input::AbstractArray, n::Int64)
    @assert 1 <= n <= length(input)
    output = missings(Float64, length(input)) # initialize vector (keep first n as missing)
    for i in eachindex(output)[n:end]
        output[i] = std(input[i-n+1:i])
    end
    return output
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Median absolute deviation
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# Wikipedia: https://en.wikipedia.org/wiki/Median_absolute_deviation
#   The median absolute deviation is a measure of statistical dispersion. Moreover, the MAD is a robust statistic,
#   being more resilient to outliers in a data set than the standard deviation. In the standard deviation, the distances from the mean are squared,
#   so large deviations are weighted more heavily, and thus outliers can heavily influence it.
#   In the MAD, the deviations of a small number of outliers are irrelevant.

function MAD_method(x)
    med = median(x)
    values = [abs(x[i] - med) for i in eachindex(x)]
    mad = median(values)
    outlyingness = [abs(x[i] - med)/mad for i in eachindex(x)]
    return outlyingness
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# HP filter (Hodrick–Prescott)
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# Wikipedia: https://en.wikipedia.org/wiki/Hodrick%E2%80%93Prescott_filter
#   The first term of the equation is the sum of the squared deviations, which penalizes the cyclical component.
#   The second term is a multiple λ of the sum of the squares of the trend component's second differences.
#   This second term penalizes variations in the growth rate of the trend component. The larger the value of λ, the higher is the penalty.
#   Hodrick and Prescott suggest 1600 as a value for λ for quarterly data.
#   Ravn and Uhlig (2002) state that λ should vary by the fourth power of the frequency observation ratio.
#   Thus, λ should equal 6.25 (1600/4^4) for annual data and 129,600 (1600*3^4) for monthly data.
#   In practice, λ = 100 for yearly data and λ = 14400 for monthly data are commonly used.

function HP(x::AbstractArray, λ::Int64)
    n = length(x)
    m = 2
    @assert n > m
    I = Diagonal(ones(n))

    # use diagm instead of spdiagm otherwise error with grouped dataframe
    D = diagm(0 => fill(1, n-m),
        -1 => fill(-2, n-m),
        -2 => fill(1, n-m) )
    @inbounds D = D[1:n,1:n-m]

    return (I + λ * D * D') \ x
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Functions for creating data for figures
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# function to create data for figure 1
function data_fig1(df::DataFrame, declarants::Vector{String}, partners::Vector{String})

    # subset dataframe
    df = subset(df, :DECLARANT_ISO => ByRow(x -> x in declarants), :PARTNER_ISO => ByRow(x -> x in partners))
    subset!(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL")) # take out total

    # compute UNIT_PRICE
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)

    # compute weights for price index
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_total = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL)
    df_join = leftjoin(df, df_total, on=cols_grouping)
    transform!(df_join, [:VALUE_IN_EUROS, :VALUE_TOTAL] => ByRow((v,s) -> v/s) => :WEIGHTS_VALUE)
    transform!(df_join, [:UNIT_PRICE, :WEIGHTS_VALUE] => ByRow((p, w) -> p * w) => :PRICE_INDEX)

    # aggregate over products
    gdf = groupby(df_join, cols_grouping)
    df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX] .=> sum, renamecols=false)

    return df
end

# function to create data for figure 2
# Notes:
#   - essentially the same function as for figure 1 but with additional grouping at CN digit level
function data_fig2(df::DataFrame, declarants::Vector{String}, partners::Vector{String}, digits::Int64)

    # subset dataframe
    df = subset(df, :DECLARANT_ISO => ByRow(x -> x in declarants), :PARTNER_ISO => ByRow(x -> x in partners))
    subset!(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL")) # take out total

    # compute UNIT_PRICE
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)

    # add 2 digit CN classification
    # Notes:
    #   - is simply the first digits
    g(x) = x[1:digits]
    df.PRODUCT_NC_digits = g.(df.PRODUCT_NC)

    # compute weights for price index
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "PRODUCT_NC_digits"]
    gdf = groupby(df, cols_grouping)
    df_total = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL)
    df_join = leftjoin(df, df_total, on=cols_grouping)
    transform!(df_join, [:VALUE_IN_EUROS, :VALUE_TOTAL] => ByRow((v,s) -> v/s) => :WEIGHTS_VALUE)
    transform!(df_join, [:UNIT_PRICE, :WEIGHTS_VALUE] => ByRow((p, w) -> p * w) => :PRICE_INDEX)

    # aggregate over products
    gdf = groupby(df_join, cols_grouping)
    df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX] .=> sum, renamecols=false)

    return df
end

# function to create data for figure 3
function data_fig3(df::DataFrame, declarants::Vector{String}, partners::Vector{String})

    # clean df
    subset!(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL")) # take out TOTAL

    # compute total exports/imports (i.e. to WORLD)
    df_total = subset(df, :DECLARANT_ISO => ByRow(x -> x in declarants), :PARTNER_ISO => ByRow(x -> x == "WORLD"))
    cols_grouping = ["DECLARANT_ISO", "FLOW", "PERIOD"]
    gdf = groupby(df_total, cols_grouping)
    df_total = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL, :QUANTITY_IN_KG => sum => :QUANTITY_TOTAL)

    # compute total exports/imports for each partner
    df_total_partner = subset(df, :DECLARANT_ISO => ByRow(x -> x in declarants), :PARTNER_ISO => ByRow(x -> x in partners))
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
    gdf = groupby(df_total_partner, cols_grouping)
    df_total_partner = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL_PARTNER, :QUANTITY_IN_KG => sum => :QUANTITY_TOTAL_PARTNER)

    # compute export/import shares
    # Notes:
    #   - shares are in percentages
    df_join = leftjoin(df_total_partner, df_total, on=["DECLARANT_ISO", "FLOW", "PERIOD"])
    transform!(df_join, [:VALUE_TOTAL_PARTNER, :VALUE_TOTAL] => ByRow((v,s) -> v/s*100) => :VALUE_SHARE)
    transform!(df_join, [:QUANTITY_TOTAL_PARTNER, :QUANTITY_TOTAL] => ByRow((v,s) -> v/s*100) => :QUANTITY_SHARE)

    df = df_join[:,Not([:VALUE_TOTAL_PARTNER, :VALUE_TOTAL, :QUANTITY_TOTAL_PARTNER, :QUANTITY_TOTAL])]

    return df
end

# function to create data for table 1 and 2
function tab1(df::DataFrame, declarants::Vector{String}, partners::Vector{String}, digits::Int64)

    # clean df
    subset!(df, :DECLARANT_ISO => ByRow(x -> x in declarants))
    subset!(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL")) # take out TOTAL

    # add 2 digit CN classification
    # Notes:
    #   - is simply the first digits
    g(x) = x[1:digits]
    df.PRODUCT_NC_digits = g.(df.PRODUCT_NC)

    # aggregate over new CN classification and PARNTERS => WORLD (to compute shares)
    subset!(df, :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)

    return df
end

# function to create data for figure 5
function fig5(df::DataFrame, declarants::Vector{String}, partners::Vector{String})

    # clean df
    subset!(df, :DECLARANT_ISO => ByRow(x -> x in declarants), :PARTNER_ISO => ByRow(x -> x in partners))
    subset!(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL")) # take out TOTAL

    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, n)

    return df
end


text:
function convert_it(read_only_numpy_arr::PyObject)
    # See https://github.com/JuliaPy/PyCall.jl/blob/master/src/pyarray.jl#L14-L24

    # instead of PyBUF_ND_STRIDED =  Cint(PyBUF_WRITABLE | PyBUF_FORMAT | PyBUF_ND | PyBUF_STRIDES)
    # See https://github.com/JuliaPy/PyCall.jl/blob/master/src/pybuffer.jl#L113
    pybuf = PyBuffer(read_only_numpy_arr, PyCall.PyBUF_FORMAT | PyCall.PyBUF_ND | PyCall.PyBUF_STRIDES)

    T, native_byteorder = PyCall.array_format(pybuf)
    sz = size(pybuf)
    strd = PyCall.strides(pybuf)
    length(strd) == 0 && (sz = ())
    N = length(sz)
    isreadonly = pybuf.buf.readonly==1
    info = PyCall.PyArray_Info{T,N}(native_byteorder, sz, strd, pybuf.buf.buf, isreadonly, pybuf)

    # See https://github.com/JuliaPy/PyCall.jl/blob/master/src/pyarray.jl#L123-L126
    PyCall.PyArray{T,N}(read_only_numpy_arr, info)
end

function get_units(unit_str)
    # workaround stuff like kg.m/s
    unit_str = replace(unit_str, "."=>"*")
    # workaround stuff like 1/m^3
    if !occursin(r"1\/([a-z]*\^[0-9]?)", unit_str)
        units = @eval @u_str $unit_str
    else
        u = replace(unit_str, r"1\/([a-z]*)\^([0-9]?)" => s"\g<1>^-\g<2>")
        units = @eval @u_str $u
    end
end


text:
if Sys.isunix() && !success(pipeline(`which curl`, devnull))
    error("curl was not found. Please install curl through your system's package manager.")
end


text:
# add
pkgs = [
    "GridInterpolations",
    "ODE",
    "JLD"
]

for pkg in pkgs
    try
        Pkg.add(pkg)
    catch e
        println("Exception when adding $(pkg): $(e)")
    end
end


text:
@testset "Clustering" begin
    g10 = CompleteGraph(10)
    for g in testgraphs(g10)
      @test @inferred(local_clustering_coefficient(g)) == ones(10)
      @test @inferred(global_clustering_coefficient(g)) == 1
      @test @inferred(local_clustering(g)) == (fill(36, 10), fill(36, 10))
      @test @inferred(triangles(g)) == fill(36, 10)
      @test @inferred(triangles(g, 1)) == 36
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
"""
    WeightedBagNode{T <: Union{AbstractMillNode, Missing}, B <: AbstractBags, W, C} <: AbstractBagNode

Structure like [`BagNode`](@ref) but allows to specify weights of type `W` of each instance.

See also: [`BagNode`](@ref), [`AbstractBagNode`](@ref), [`AbstractMillNode`](@ref), [`BagModel`](@ref).
"""
struct WeightedBagNode{T <: Maybe{AbstractMillNode}, B <: AbstractBags, W, C} <: AbstractBagNode
    data::T
    bags::B
    weights::Vector{W}
    metadata::C

    function WeightedBagNode(d::T, b::B, w::Vector{W}, m::C=nothing) where {T <: Maybe{AbstractMillNode}, B <: AbstractBags, W, C}
        ismissing(d) && any(length.(b.bags) .!= 0) && error("WeightedBagNode with nothing in data cannot have a non-empty bag")
        new{T, B, W, C}(d, b, w, m)
    end
end

"""
    WeightedBagNode(d::Union{AbstractMillNode, Missing}, b::AbstractBags, w::Vector, m=nothing)
    WeightedBagNode(d::Union{AbstractMillNode, Missing}, b::AbstractVector, w::Vector, m=nothing)

Construct a new [`WeightedBagNode`](@ref) with data `d`, bags `b`, weights `w` and metadata `m`. If `b` is an `AbstractVector`, [`Mill.bags`](@ref) is applied first.

# Examples
```jldoctest
julia> BagNode(ArrayNode(NGramMatrix(["s1", "s2"])), bags([1:2, 0:-1]), [0.2, 0.8])
BagNode 	# 2 obs, 184 bytes
  └── ArrayNode(2053×2 NGramMatrix with Int64 elements) 	# 2 obs, 140 bytes

julia> BagNode(ArrayNode(zeros(2, 2)), [1, 2], [1, 2])
BagNode 	# 2 obs, 160 bytes
  └── ArrayNode(2×2 Array with Float64 elements) 	# 2 obs, 80 bytes
```

See also: [`BagNode`](@ref), [`AbstractBagNode`](@ref), [`AbstractMillNode`](@ref), [`BagModel`](@ref).
"""
WeightedBagNode(d::Maybe{AbstractMillNode}, b::AbstractVector, weights::Vector, metadata=nothing) =
    WeightedBagNode(d, bags(b), weights, metadata)

mapdata(f, x::WeightedBagNode) = WeightedBagNode(mapdata(f, x.data), x.bags, x.weights, x.metadata)

dropmeta(x::WeightedBagNode) = WeightedBagNode(dropmeta(x.data), x.bags, x.weights)

function Base.getindex(x::WeightedBagNode{T, B, W}, i::VecOrRange{<:Int}) where {T, B, W}
    nb, ii = remapbags(x.bags, i)
    emptyismissing() && isempty(ii) && return(WeightedBagNode(missing, nb, W[], nothing))
    WeightedBagNode(subset(x.data,ii), nb, subset(x.weights, ii), subset(x.metadata, ii))
end

function reduce(::typeof(catobs), as::Vector{<:WeightedBagNode})
    d = filter(!ismissing, data.(as))
    md = filter(!isnothing, metadata.(as))
    bags = reduce(vcat, [n.bags for n in as])
    WeightedBagNode(reduce(catobs, d), bags, reduce(vcat, [a.weights for a in as]), reduce(catobs, md))
end

function removeinstances(a::WeightedBagNode, mask)
    WeightedBagNode(subset(a.data, findall(mask)),
        adjustbags(a.bags, mask),
        subset(a.weights, findall(mask)), a.metadata)
end

Base.hash(e::WeightedBagNode, h::UInt) = hash((e.data, e.bags, e.weights, e.metadata), h)
(e1::WeightedBagNode == e2::WeightedBagNode) =
    e1.data == e2.data && e1.bags == e2.bags && e1.weights == e2.weights && e1.metadata == e2.metadata
Base.isequal(e1::WeightedBagNode, e2::WeightedBagNode) =
    isequal(e1.data, e2.data) && isequal(e1.bags, e2.bags) && isequal(e1.weights, e2.weights) && isequal(e1.metadata, e2.metadata)


text:
module TestMerCounting

using MerCounting, BioSequences, Test

include("MerCount.jl")

end # module