text:
immutable DBNModel
    BN            :: BayesNet
    statsvec      :: Vector{Matrix{Float64}} # each matrix is [r×q], nvar_instantiations × nparent_instantiations
    features      :: Vector{AbstractFeature}
    discretizers  :: Vector{AbstractDiscretizer}
    istarget      :: BitVector # whether a feature is a target feature
end

function Base.print(io::IO, model::DBNModel)

    target_lat = get_target_lat(model)
    target_lon = get_target_lon(model)

    println(io, "DBNModel")
    println(io, "\t", symbol(target_lat), " <- ", map(f->symbol(f), get_indicators_for_target(model, target_lat)))
    println(io, "\t", symbol(target_lon), " <- ", map(f->symbol(f), get_indicators_for_target(model, target_lon)))
end

function dbnmodel{R<:Real, D<:AbstractDiscretizer}(
    BN:: BayesNet,
    statsvec::Vector{Matrix{R}},
    discretizerdict::Dict{Symbol, D},
    targets::Vector{AbstractFeature},
    indicators::Vector{AbstractFeature}
    )

    features = AbstractFeature[symbol2feature(sym) for sym in names(BN)]
    discretizers = AbstractDiscretizer[discretizerdict[sym] for sym in names(BN)]
    istarget = falses(length(features))
    for (i,f) in enumerate(features)
        if in(f, targets)
            istarget[i] = true
        elseif !in(f, indicators)
            error("Feature not in targets or indicators")
        end
    end

    statsvec_float = convert(Vector{Matrix{Float64}}, statsvec)
    DBNModel(BN, statsvec_float, features, discretizers, istarget)
end
function dbnmodel{R<:Real, D<:AbstractDiscretizer}(
    BN::BayesNet,
    statsvec::Vector{Matrix{R}},
    discretizers::Vector{D},
    targets::Vector{AbstractFeature},
    indicators::Vector{AbstractFeature},
    )

    features = AbstractFeature[symbol2feature(sym) for sym in names(BN)]
    istarget = falses(length(features))
    for (i,f) in enumerate(features)
        if in(f, targets)
            istarget[i] = true
        elseif !in(f, indicators)
            error("Feature $(symbol(f)) neither in targets nor indicators")
        end
    end

    statsvec_float = convert(Vector{Matrix{Float64}}, statsvec)
    DBNModel(BN, statsvec_float, features, discretizers, istarget)
end
function dbnmodel{R<:Real, D<:AbstractDiscretizer}(
    adj::BitMatrix, # adj[i,j] = true means i → j
    statsvec::Vector{Matrix{R}},
    discretizerdict::Dict{Symbol, D},
    targets::Vector{AbstractFeature},
    indicators::Vector{AbstractFeature},
    )

    BN = build_bn(statsvec, targets, indicators, adj)
    dbnmodel(BN, statsvec, discretizerdict, targets, indicators)
end
function dbnmodel(modelstats::Dict{AbstractString, Any})

    discretizerdict = modelstats["binmaps"]
    targets    = modelstats["targets"]
    indicators = modelstats["indicators"]
    stats      = modelstats["statistics"]
    adj        = modelstats["adjacency"] # adj[i,j] = true means i → j

    dbnmodel(adj, stats, discretizerdict, convert(Vector{AbstractFeature}, targets),
                        convert(Vector{AbstractFeature}, indicators))
end
function dbnmodel(modelpstats_file::AbstractString)
    emstats = load(modelpstats_file)
    dbnmodel(emstats)
end

function build_bn{R<:Real}(
    statsvec   :: Vector{Matrix{R}},
    targets    :: Vector{AbstractFeature},
    indicators :: Vector{AbstractFeature}, # technically all of the non-target indicators
    adj        :: BitMatrix
    )

    @assert(size(adj,1) == size(adj,2)) # is square

    n_targets = length(targets)
    n_indicators = length(indicators)

    bnnames = Array(Symbol, n_targets + n_indicators)
    for (i,f) in enumerate(targets)
        bnnames[i] = symbol(f)
    end
    for (i,f) in enumerate(indicators)
        bnnames[i+n_targets] = symbol(f)
    end

    @assert(length(unique(bnnames)) == length(bnnames))
    n_nodes = length(bnnames)
    @assert(n_nodes == size(adj, 1) == size(adj, 2))

    BN = BayesNet(bnnames)

    r_arr = Array(Int, n_nodes)
    for (i,node_sym) in enumerate(bnnames)
        stats = statsvec[i]
        r, q = size(stats) # r = num node instantiations, q = num parental instantiations
        states = collect(1:r)
        BN.nodes[i].domain = BayesNets.CPDs.DiscreteDomain(states)
        r_arr[i] = r
    end

    for (node,node_sym) in enumerate(bnnames)

        stats = statsvec[node]
        r, q = size(stats)
        states = collect(1:r)

        stats .+= 1 # NOTE(tim): adding uniform prior
        probabilities = stats ./ sum(stats,1)

        # set any parents & populate probability table
        n_parents = sum(adj[:,node])
        if n_parents > 0
            bnparents = bnnames[adj[:,node]]
            for pa in bnparents
                add_edge!(BN, pa, node_sym)
            end

            # populate probability table
            assignments = BayesNets.assignment_dicts(BN, bnparents)
            # parameterFunction = BayesNets.discrete_parameter_function(assignments, vec(probabilities), r)
            # setCPD!(BN, node_sym, CPDs.Discrete(states, parameterFunction))

            parameterlookup = BayesNets.discrete_parameter_dict(assignments, vec(probabilities), r)
            set_CPD!(BN, node_sym, CPDs.DiscreteDict(states, parameterlookup))
        else
            # no parents
            # setCPD!(BN, node_sym, CPDs.Discrete(states, vec(probabilities)))

            set_CPD!(BN, node_sym, CPDs.DiscreteStatic(states, vec(probabilities)))
        end

    end

    return BN
end
# function build_bn{R<:Real}(
#     statsvec   :: Vector{Matrix{R}},
#     targets    :: Vector{AbstractFeature},
#     indicators :: Vector{AbstractFeature}, # technically all of the non-target indicators
#     adj        :: BitMatrix
#     )

#     @assert(size(adj,1) == size(adj,2)) # is square

#     n_targets = length(targets)
#     n_indicators = length(indicators)

#     bnnames = Array(Symbol, n_targets + n_indicators)
#     for (i,f) in enumerate(targets)
#         bnnames[i] = symbol(f)
#     end
#     for (i,f) in enumerate(indicators)
#         bnnames[i+n_targets] = symbol(f)
#     end

#     @assert(length(unique(bnnames)) == length(bnnames))
#     n_nodes = length(bnnames)
#     @assert(n_nodes == size(adj, 1) == size(adj, 2))

#     BN = BayesNet(bnnames)

#     r_arr = Array(Int, n_nodes)
#     for (node,node_sym) in enumerate(bnnames)
#         stats = statsvec[node]
#         r, q = size(stats) # r = num node instantiations, q = num parental instantiations
#         states = collect(1:r)
#         set_domain!(BN, node_sym, DiscreteDomain(states))
#         r_arr[node] = r
#     end

#     for (node,node_sym) in enumerate(bnnames)

#         stats = statsvec[node]
#         r, q = size(stats)
#         states = collect(1:r)

#         stats .+= 1 # NOTE(tim): adding uniform prior
#         probabilities = stats ./ sum(stats,1)

#         # set any parents & populate probability table
#         n_parents = sum(adj[:,node])
#         if n_parents > 0
#             bnparents = bnnames[adj[:,node]]
#             for pa in bnparents
#                 add_edge!(BN, pa, node_sym)
#             end

#             # populate probability table
#             assignments = BayesNets.assignment_dicts(BN, bnparents)
#             # parameterFunction = BayesNets.discrete_parameter_function(assignments, vec(probabilities), r)
#             # setCPD!(BN, node_sym, CPDs.Discrete(states, parameterFunction))

#             parameterlookup = BayesNets.discrete_parameter_dict(assignments, vec(probabilities), r)
#             setCPD!(BN, node_sym, CPDs.DiscreteDictCPD(states, parameterlookup))
#         else
#             # no parents
#             # setCPD!(BN, node_sym, CPDs.Discrete(states, vec(probabilities)))

#             setCPD!(BN, node_sym, CPDs.DiscreteStaticCPD(states, vec(probabilities)))
#         end

#     end

#     return BN
# end

indexof(f::Symbol, model::DBNModel) = model.BN.name_to_index[f]
indexof(f::AbstractFeature, model::DBNModel) = model.BN.name_to_index[symbol(f)]
is_parent(model::DBNModel, parent::Int, child::Int) = in(parent, in_neighbors(model.BN.dag, child))
is_parent(model::DBNModel, parent::Symbol, child::Symbol) = is_parent(model, model.BN.name_to_index[parent], model.BN.name_to_index[child])
function parent_indeces(varindex::Int, model::DBNModel)
    parent_names = BayesNets.parents(model.BN, names(model.BN)[varindex])
    retval = Array(Int, length(parent_names))
    for (i, name) in enumerate(parent_names)
        retval[i] = model.BN.name_to_index[name]
    end
    retval
end

get_targets(model::DBNModel) = model.features[model.istarget]
function get_target_lat(model::DBNModel, targets::Vector{AbstractFeature}=get_targets(model))
    ind = findfirst(f->isa(f, FeatureTargetLat), targets)
    targets[ind]
end
function get_target_lon(model::DBNModel, targets::Vector{AbstractFeature}=get_targets(model))
    ind = findfirst(f->isa(f, FeatureTargetLon), targets)
    targets[ind]
end

AutomotiveDrivingModels.get_indicators(model::DBNModel) = model.features[!model.istarget]
function get_indicators_for_target(model::DBNModel, i::Int)
    indicator_indeces = parent_indeces(i, model)
    model.features[indicator_indeces]
end
function get_indicators_for_target(model::DBNModel, target::AbstractFeature)
    i = indexof(target, model)
    get_indicators_for_target(model, i)
end
get_indicators_lat(model::DBNModel) = get_indicators_for_target(model, get_target_lat(model))
get_indicators_lon(model::DBNModel) = get_indicators_for_target(model, get_target_lon(model))

function find_target_indeces{F<:AbstractFeature}(targetset::ModelTargets, features::Vector{F})
    ind_lat = findfirst(f->f==targetset.lat, features)
    ind_lon = findfirst(f->f==targetset.lon, features)
    @assert(ind_lat  != 0)
    @assert(ind_lon  != 0)
    (ind_lat, ind_lon)
end

get_num_vertices(model::DBNModel) = nv(model.BN.dag)
get_num_edges(model::DBNModel) = ne(model.BN.dag)

get_total_sample_count(model::DBNModel) = sum(model.statsvec[1])
function get_bin_counts(model::DBNModel)
    n_nodes = get_num_vertices(model)
    r_arr = Array(Int, n_nodes)
    for i = 1 : n_nodes
        r_arr[i] = size(model.statsvec[i], 1)
    end
    r_arr
end
function get_marginal_probability(model::DBNModel, varindex::Int)
    binprobs = vec(sum(model.statsvec[varindex], 2))
    binprobs ./= sum(binprobs)
end
function get_counts_for_assignment(
    model::DBNModel,
    targetind::Int,
    parentindeces::Vector{Int},
    parentassignments::Vector{Int},
    bincounts::Vector{Int}
    )

    dims = tuple(bincounts[parentindeces]...)
    subs = tuple(parentassignments...)
    j = sub2ind(dims, subs...)
    copy(model.statsvec[targetind][:,j])
end
function get_counts_for_assignment(
    model::DBNModel,
    targetind::Int,
    assignments::Dict{Symbol, Int},
    bincounts::Vector{Int}
    )

    parentindeces = parent_indeces(targetind, model)
    nparents = length(parentindeces)
    parentassignments = Array(Int, nparents)
    for (i,ind) in enumerate(parentindeces)
        parentassignments[i] = assignments[model.BN.nodes[ind].name]
    end
    get_counts_for_assignment(model, targetind, parentindeces, parentassignments, bincounts)
end

function encode!(assignment::Dict{Symbol,Int}, model::DBNModel, observations::Dict{Symbol,Float64})
    # take each observation and bin it appropriately
    # returns a Dict{Symbol,Int}

    # TODO(tim): ensure order is correct
    for (i,istarget) in enumerate(model.istarget)
        if !istarget
            sym = model.BN.nodes[i].name
            val = observations[sym]
            try
                encode(model.discretizers[i], val)
            catch
                println("failed to encode ", sym, " of value ", val, "  ", model.discretizers[i])
            end
            assignment[sym] = encode(model.discretizers[i], val)
        end
    end
    assignment
end
function sample!(model::DBNModel, assignment::Dict{Symbol, Int}, ordering::Vector{Int}=topological_sort_by_dfs(model.BN.dag))
    #=
    Run through nodes in topological order, building the instantiation vector as we go
    We use nodes we already know to condition on the distribution for nodes we do not
    Modifies assignment to include newly sampled symbols
    =#

    # for name in model.BN.names[ordering]
    #     if !haskey(assignment, name)
    #         assignment[name] = BayesNets.rand(BayesNets.cpd(model.BN, name), assignment)
    #     end
    # end

    for node in model.BN.nodes[ordering]
        name = node.name
        cpd = BayesNets.cpd(model.BN, name)

        p = BayesNets.probvec(cpd, assignment)
        r = rand()
        i = 1
        p_tot = 0.0
        while p_tot + p[i] < r && i < length(p)
            p_tot += p[i]
            i += 1
        end
        assignment[name] = cpd.domain[i]
    end

    assignment
end
function sample_unset!(model::DBNModel, assignment::Dict{Symbol, Int}, ordering::Vector{Int}=topological_sort_by_dfs(model.BN.dag); missing_value::Int=-1)
    #=
    Run through nodes in topological order, building the instantiation vector as we go
    We use nodes we already know to condition on the distribution for nodes we do not
    Modifies assignment to include newly sampled symbols
    Only sample values set to missing_value
    =#

    for node in model.BN.nodes[ordering]

        name = node.name

        if get(assignment, name, missing_value) == missing_value

            cpd = BayesNets.cpd(model.BN, name)

            p = BayesNets.probvec(cpd, assignment)
            r = rand()
            i = 1
            p_tot = 0.0
            while p_tot + p[i] < r && i < length(p)
                p_tot += p[i]
                i += 1
            end
            assignment[name] = cpd.domain[i]
        end
    end

    assignment
end
function sample_and_logP!(
    model::DBNModel,
    assignment::Dict{Symbol, Int},
    logPs::Dict{Symbol, Float64},

    ordering::Vector{Int}=topological_sort_by_dfs(model.BN.dag),
    )

    for node in model.BN.nodes[ordering]
        name = node.name

        cpd = BayesNets.cpd(model.BN, name)

        p = BayesNets.probvec(cpd, assignment)

        r = rand()
        i = 1
        p_tot = 0.0
        while p_tot + p[i] < r && i < length(p)
            p_tot += p[i]
            i += 1
        end

        assignment[name] = cpd.domain[i]
        logPs[name] = log(p[i])
    end

    logPs
end
function calc_log_probability_of_assignment(model::DBNModel, assignment::Dict{Symbol, Int}, sym::Symbol)
    # returns the discrete log probability of the given bin assignment
    cpd = BayesNets.cpd(model.BN, sym)
    b = assignment[sym]
    p = cpd.parameterFunction(assignment)
    log(p[b])
end
function calc_probability_distribution_over_assignments!(
    dest::Vector{Float64},
    model::DBNModel,
    assignment::Dict{Symbol, Int},
    target::Symbol
    )

    # NOTE (tim): cpd.parameterFunction(assignment) returns the actual probability vector, not a copy
    cpd = BayesNets.cpd(model.BN, target)
    copy!(dest, BayesNets.probvec(cpd, assignment))
    @assert(abs(sum(dest)-1.0) < 0.0001) # probability sums to 1.0
    dest
end
function calc_probability_for_uniform_sample_from_bin(
    discrete_prob::Float64,
    ::CategoricalDiscretizer, # disc
    ::Int # bin
    )

    discrete_prob
end
function calc_probability_for_uniform_sample_from_bin(
    discrete_prob::Float64,
    disc::LinearDiscretizer,
    bin::Int
    )

    # widthtot = totalwidth(disc)
    widthbin = binwidth(disc, bin)

    discrete_prob / widthbin # * widthtot
end
function calc_probability_for_uniform_sample_from_bin(
    discrete_prob::Float64,
    disc::HybridDiscretizer,
    bin::Int
    )

    if bin ≤ disc.lin.nbins
        calc_probability_for_uniform_sample_from_bin(discrete_prob, disc.lin, bin)
    else
        discrete_prob
    end
end

function export_to_text(model::DBNModel, filename::AbstractString)
    # export a bayesian network to an encounter definition file

    n = model.n # the number of nodes

    open(filename, "w") do fout

        # Labels: the node names
        ###############################################
        println(fout, "# labels")
        for i = 1 : n
            @printf(fout, "\"%s\"", string(model.labels[i]))
            if i < n
                @printf(fout, ", ")
            else
                @printf(fout, "\n")
            end
        end

        # G: the graphical structure
        ##################################################
        println(fout, "# G: graphical structure")
        for i = 1 : n
            for j = 1 : n
                @printf(fout, "%c ", model.G[i,j] ? '1' : '0')
            end
            @printf(fout, "\n")
        end

        # r: number of bins for each variable
        ##################################################
        println(fout, "# r: number of bins")
        for i = 1 : n
            @printf(fout, "%d ", model.r[i])
        end
        @printf(fout, "\n")

        # N: the sufficient statistics, in integer form
        ##################################################
        println(fout, "# N: sufficient statistics")
        for i = 1 : n
            stats = model.N[i]
            r, q = size(stats)
            for b = 1:q, a = 1:r
                @printf(fout, "%d ", stats[a,b])
            end
        end
        @printf(fout, "\n")

        # discretizers, '*' if discrete
        ##################################################
        @printf(fout, "# discretizers\n")
        for i = 1 : n
            bmap = model.discretizers[i]

            if isa(bmap, CategoricalDiscretizer)
                @printf(fout, "*")
            elseif isa(bmap, LinearDiscretizer)
                for edge in bmap.binedges
                    @printf(fout, "%f ", edge)
                end
            elseif isa(bmap, HybridDiscretizer)
                for edge in bmap.lin.binedges
                    @printf(fout, "%f ", edge)
                end
            else
                error("invalid bmap type $(typeof(bmap))")
            end
            @printf(fout, "\n")
        end
    end
end
function print_structure(model::DBNModel)
    target = get_target_lat(model)
    println("target lat: ", symbol(target))
    for indicator in get_indicators_lat(model, target)
        println("\t", symbol(indicator))
    end

    target = get_target_lon(model)
    println("target lon: ", symbol(target))
    for indicator in get_indicators_lon(model, target)
        println("\t", symbol(indicator))
    end
end

function construct_model_adjacency_matrix(
    net_ind_lat :: Int,
    net_ind_lon :: Int,
    net_parent_indices_lat :: Vector{Int},
    net_parent_indices_lon :: Vector{Int},
    num_net_features :: Int
    )

    # [parent_index, child_index] = true for parent -> child edges

    adj = falses(num_net_features, num_net_features)

    adj[net_parent_indices_lat, net_ind_lat] = true
    adj[net_parent_indices_lon, net_ind_lon] = true

    adj
end


text:
function rel_err(sol, ref)
    maximum(abs.((sol.d[:,end] .- ref) ./ ref))
end

ext: jl


text:
#
# Vadd
#
# Element wise addition of three vectors at a time (R=A+B+C)
# Asks the user to select a device at runtime
#
# History: Initial version based on vadd.c, written by Tim Mattson, June 2011
#          Ported to C++ Wrapper API by Benedict Gaster, September 2011
#          Updated to C++ Wrapper API v1.2 by Tom Deakin and Simon McIntosh-Smith, October 2012
#          Ported to Python by Tom Deakin, July 2013
#          Ported to Julia by Jake Bolewski, Nov 2013

using OpenCL

# tolerance used in floating point comparisons
TOL = 1e-3

# length of vectors a, b, c
LENGTH = 1024

# Kernel: vadd
#
# To compute the elementwise sum r = a + b + c
#
# Input: a, b and c float vectors of length count
# Output r float vector of length count holding the sum a + b + cs

kernelsource = "
__kernel void vadd(
    __global float* a,
    __global float* b,
    __global float* c,
    __global float* r,
    const unsigned int count)
{
    unsigned int i = get_global_id(0);
    if (i < count)
        r[i] = a[i] + b[i] + c[i];
}"

# create a compute context
device, ctx, queue = cl.create_compute_context()

# create the compute program and build it
program = cl.Program(ctx, source=kernelsource) |> cl.build!

# create a, b and c vectors and fill with random float values
# (the result array will be created when reading back from the device)
h_a = rand(Float32, LENGTH)
h_b = rand(Float32, LENGTH)
h_c = rand(Float32, LENGTH)

d_a = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_a)
d_b = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_b)
d_c = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_c)

# create the output (r) buffer in device memory
d_r = cl.Buffer(Float32, ctx, :w, LENGTH)

# create the kernel
vadd = cl.Kernel(program, "vadd")

# execute the kernel over the entire range of the input
vadd[queue, size(h_a)](d_a, d_b, d_c, d_r, UInt32(LENGTH))

# read the results back from the compute device
# by convention..
# cl.(action) calls are blocking
# cl.enqueue_(action) calll are async/non-blocking
h_r = cl.read(queue, d_r)

# test the results
correct = 0
for i in 1:LENGTH
    tmp = h_a[i] + h_b[i] + h_c[i]
    # compute the deviation of expected and output result
    tmp -= h_r[i]
    if tmp^2 < TOL^2
        global correct += 1
    else
        println("tmp $tmp h_a $(h_a[i]) h_b $(h_b[i]) h_c $(h_c[i]) h_r $(h_r[i])")
    end
end

# summarize results
println("3 vector adds to find F=A+B+C: $correct out of $LENGTH results were correct")


text:
const USE_GPU = false
using LinearAlgebra
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using GLMakie, Makie
using Printf, ImplicitGlobalGrid
import MPI

# macros to avoid array allocation
macro qx(ix,iy,iz) esc(:( -D_dx*(H[$ix+1,$iy+1,$iz+1] - H[$ix,$iy+1,$iz+1]) )) end
macro qy(ix,iy,iz) esc(:( -D_dy*(H[$ix+1,$iy+1,$iz+1] - H[$ix+1,$iy,$iz+1]) )) end
macro qz(ix,iy,iz) esc(:( -D_dz*(H[$ix+1,$iy+1,$iz+1] - H[$ix+1,$iy+1,$iz]) )) end

"""
    norm_g(A)

Computes the l2 norm of `A`, an MPI array.
"""
norm_g(A) = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))

"""
    compute_dual_time!(H, H2, Hold, dHdt, _dt, dtau, _dx, _dy, _dz, D_dx, D_dy,
    D_dz, dmp, size_H1_2, size_H2_2, size_H3_2)

ImplicitGlobalGrid kernel that computes a virtual time step of the diffusion 3D
process.
"""
@parallel_indices (ix,iy,iz) function compute_dual_time!(H, H2, Hold, dHdt, _dt, dtau, _dx, _dy, _dz, D_dx, D_dy, D_dz, dmp, size_H1_2, size_H2_2, size_H3_2)

    if (ix<=size_H1_2 && iy<=size_H2_2 && iz<=size_H3_2)
        dHdt[ix,iy,iz] = dmp * dHdt[ix,iy,iz] -
            (H[ix+1,iy+1,iz+1] - Hold[ix+1,iy+1,iz+1]) * _dt -
                                        ((@qx(ix+1,iy,iz) - @qx(ix,iy,iz))*_dx +
                                         (@qy(ix,iy+1,iz) - @qy(ix,iy,iz))*_dy +
                                         (@qz(ix,iy,iz+1) - @qz(ix,iy,iz))*_dz)

        H2[ix+1,iy+1,iz+1] = H[ix+1,iy+1,iz+1] + dtau * dHdt[ix,iy,iz]
    end
    return
end

"""
    compute_residual!(Rh, H, Hold, _dt, _dx,
    _dy, _dz, D_dx, D_dy, D_dz, size_H1_2, size_H2_2, size_H3_2)

ImplicitGlobalGrid kernel that computes the residual of the dual-time diffusion
equation, used to know when to stop the virtual time iterations.
"""
@parallel_indices (ix,iy,iz) function compute_residual!(Rh, H, Hold, _dt, _dx, _dy, _dz, D_dx, D_dy, D_dz, size_H1_2, size_H2_2, size_H3_2)

    if (ix<=size_H1_2 && iy<=size_H2_2 && iz<=size_H3_2)
        Rh[ix,iy,iz] = -(H[ix+1,iy+1,iz+1] - Hold[ix+1,iy+1,iz+1]) * _dt -
                                        ((@qx(ix+1,iy,iz) - @qx(ix,iy,iz))*_dx +
                                         (@qy(ix,iy+1,iz) - @qy(ix,iy,iz))*_dy +
                                         (@qz(ix,iy,iz+1) - @qz(ix,iy,iz))*_dz)
    end
    return
end

"""
    diffusion_3D(; ttot, n, nout, do_visu=false)

Creates the 3D steady-state diffusion problem and solves, using a dual-time
approach. All the physics and numerics are defined in this function.

# Arguments
    - `ttot`: the total physical time.
    - `n`: the number of discretised cells per dimension.
    - `nout`: the number of virtual time iterations before we check the residual.
    - `do_visu`: if true, each physical time step will be ploted and saved as
        an image. In order to produce a video from those images (stored in
        `plots/part-1`), please execute the script `create_anim.sh` located in
        the same folder.
# Return values
    - `H`: The solution arrray (Water surface height in m).
    - `xc`: The x-coord vector.
"""
@views function diffusion_3D(;
        # Physics
        ttot = 1.0,
        # Numerics
        n = 32,
        nout = 10,
        tol = 1e-8,
        # Visualisation
        do_visu=false)
    # Physics
    Lx, Ly, Lz  = 10.0, 10.0, 10.0
    D       = 1.0
    # Numerics
    nx, ny, nz  = n, n, n # number of grid points
    ϵ       = tol
    maxIter = 1e5
    # Derived numerics
    me, dims = init_global_grid(nx, ny, nz)  # Initialization of MPI and more...
    @static if USE_GPU select_device() end  # select one GPU per MPI local rank (if >1 GPU per node)
    dx, dy, dz = Lx/nx_g(), Ly/ny_g(), Lz/nz_g()
    dt      = 0.2
    _dt     = 1.0 / dt
    dtau    = (1.0/(min(dx, dy, dz)^2/D/6.1) + 1.0/dt)^-1
    nt      = cld(ttot, dt)
    xc, yc, zc  = LinRange(dx/2, Lx-dx/2, nx), LinRange(dy/2, Ly-dy/2, ny), LinRange(dz/2, Lz-dz/2, nz)
    D_dx    = D/dx
    D_dy    = D/dy
    D_dz    = D/dz
    _dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz
    # Array initialisation
    H0    = zeros(nx,ny,nz)
    H0    = Data.Array([2.0 .* exp(-(x_g(ix,dx,H0)+dx/2 -Lx/2)^2 -(y_g(iy,dy,H0)+dy/2 -Ly/2)^2 -(z_g(iz,dz,H0)+dz/2 -Lz/2)^2) for ix=1:size(H0,1), iy=1:size(H0,2), iz=1:size(H0,3)])
    H     = copy(H0)
    H2    = copy(H)
    Hold  = copy(H)
    Rh    = @zeros(nx-2,ny-2,nz-2)
    dHdt  = @zeros(nx-2,ny-2,nz-2)
    size_H1_2, size_H2_2, size_H3_2 = size(H,1)-2, size(H,2)-2, size(H,3)-2
    len_Rh_g = ((nx-2-2)*dims[1]+2)*((ny-2-2)*dims[2]+2)*((nz-2-2)*dims[3]+2)
    # Preparation of visualisation
    if do_visu
        GLMakie.activate!()

        # Create directory if it does not exist
        if (~isdir("plots/part-1")) mkpath("plots/part-1") end

        if (me==0) ENV["GKSwstype"]="nul"; if isdir("viz2D_mxpu_out")==false mkdir("viz2D_mxpu_out") end; end

        nx_v, ny_v, nz_v = (nx-2)*dims[1], (ny-2)*dims[2], (nz-2)*dims[3]
        if (nx_v*ny_v*nz_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        H_v   = zeros(nx_v, ny_v, nz_v) # global array for visu
        H_inn = zeros(nx-2, ny-2, nz-2) # no halo local array for visu
        Xi_g, Yi_g, Zi_g = LinRange(dx+dx/2, Lx-dx-dx/2, nx_v), LinRange(dy+dy/2, Ly-dy-dy/2, ny_v), LinRange(dz+dz/2, Lz-dz-dz/2, nz_v)
    end

    t_tic = 0.0; niter_1 = 0
    # Time loop
    t = 0.0; it = 0; ittot = 0;
    dmp = 1-18/nx_g()
    while t<ttot
        if (it==1) t_tic = Base.time(); niter_1 = ittot end
        iter = 0; err = 2*ϵ

        while err > ϵ && iter < maxIter
            # @hide_communication (8, 2) begin
                @parallel compute_dual_time!(H, H2, Hold, dHdt, _dt, dtau, _dx, _dy, _dz, D_dx, D_dy, D_dz, dmp, size_H1_2, size_H2_2, size_H3_2)
                H, H2 = H2, H
                update_halo!(H)
            # end
            iter += 1
            if (iter % nout == 0)
                @parallel compute_residual!(Rh, H, Hold, _dt, _dx, _dy, _dz, D_dx, D_dy, D_dz, size_H1_2, size_H2_2, size_H3_2)
                err = norm_g(Rh)/sqrt(len_Rh_g)
            end
        end
        ittot += iter; it += 1; t += dt
        Hold .= H
        if isnan(err) error("Error is NaN!") end

        # Visualize
        if do_visu
            H_inn .= H[2:end-1,2:end-1,2:end-1]; gather!(H_inn, H_v)
            if (me==0)
                fig = Figure(resolution = (1200,1200))
                ax1 = Axis3(fig[1,1:2], perspectiveness = 0.5, azimuth = (it/nout)*0.2,
                    elevation = 0.57, aspect = (1,1,1), title="Diffusion 3D domain (t=$(round(t*100)/100)s)")
                ax2 = Axis(fig[2,1], title="z=16 slice")
                ax3 = Axis(fig[2,2], title="x=16 slice")
                volume!(ax1, Xi_g, Yi_g, Zi_g, H_v; colormap=:viridis, transparency=true)
                heatmap!(fig[2,1], Xi_g, Yi_g, H_v[:,:,16]')
                heatmap!(fig[2,2], Xi_g, Yi_g, H_v[16,:,:]')
                # Add information
                image_num = lpad(it, 3, "0")
                save("plots/part-1/plot_$(image_num).png", fig, px_per_unit = 2)
            end
        end
    end
    t_toc = Base.time() - t_tic
    niter = ittot - niter_1
    A_eff = 5/1e9*nx_g()*ny_g()*nz_g()*sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it  = t_toc/niter                          # Execution time per iteration [s]
    T_eff = A_eff/t_it                           # Effective memory throughput [GB/s]
    if (me==0) @printf("Time = %1.5f sec, T_eff = %1.6f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits=3), niter) end

    # Retrieve the global domain
    finalize_global_grid()
    return H, xc
end

# diffusion_3D(; do_visu=false)


text:
while true
	if read(ExplorerHat.analog[1]) < 0.5
		on(ExplorerHat.output[1])
	else
		off(ExplorerHat.output[1])
        end
	sleep(1)
end


text:
module FATEStools

# Include functions pertaining to fates and hlm inputdata (e.g. paramfiles, metdata, etc)
include("inputdata.jl")

export ncmeta_addglobalattr, ncmeta_globaltemplate, ncmeta_vartemplate, ncmeta_checkvar, ncmeta_addvarattr

end


text:
#!/usr/bin/env julia

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using JuliaFormatter

headbranch = get(ARGS, 1, "master")

for filename in
    readlines(`git diff --name-only --diff-filter=AM $headbranch...`)
    endswith(filename, ".jl") || continue

    format(filename, verbose = true, indent = 4, margin = 80)
end


text:
using Documenter
using DataStructures

makedocs(
    format = :html,
    sitename = "DataStructures.jl",
    pages = [
        "index.md",
        "deque.md",
        "circ_buffer.md",
        "circ_deque.md",
        "stack_and_queue.md",
        "priority-queue.md",
        "accumulators.md",
        "disjoint_sets.md",
        "heaps.md",
        "ordered_containers.md",
        "default_dict.md",
        "trie.md",
        "linked_list.md",
        "intset.md",
        "sorted_containers.md",
    ],
    modules = [DataStructures]
)

deploydocs(
    repo = "github.com/JuliaCollections/DataStructures.jl.git",
    julia  = "0.6",
    latest = "master",
    target = "build",
    deps = nothing,  # we use the `format = :html`, without `mkdocs`
    make = nothing,  # we use the `format = :html`, without `mkdocs`
)


text:
using DiffinDiffsBase: RefArray, validpool, scaledlabel

@testset "ScaledArray" begin
    refs = repeat(1:5, outer=2)
    pool = Date(1):Year(1):Date(5)
    invpool = Dict(p=>i for (i, p) in enumerate(pool))
    sa = ScaledArray(RefArray(refs), pool, invpool)
    ssa = view(sa, 3:5)
    @test sa.refs === refs

    @test scale(sa) == Year(1)
    @test scale(ssa) == Year(1)
    @test size(sa) == (10,)
    @test IndexStyle(typeof(sa)) == IndexLinear()

    x = 1:10
    @test validpool(x, Int, 10, -1, nothing, true) == 10:-1:1
    @test_throws ArgumentError validpool(x, Int, 1, -1, 10, true)
    p = PooledArray(x)
    p[1] = 2
    p[10] = 9
    @test validpool(p, Int, nothing, 1, nothing, true) == 1:10
    @test_throws ArgumentError validpool(p, Int, 2, 1, 9, true)
    @test validpool(p, Int, 2, 1, nothing, false) == 2:9
    @test_throws ArgumentError validpool(p, Int, 2, 1, 8, false)
    s = ScaledArray(1.0:2.0:9.0, 1.0, 2, 10)
    ss = validpool(s, Float64, 1, 2, 9, true)
    @test ss == 1.0:2.0:9.0
    @test eltype(ss) == Float64
    @test_throws ArgumentError validpool(1:5, Int, 1, nothing, 1, true)
    @test_throws ArgumentError validpool(1:5, Int, 1, Year(1), 1, true)
    @test_throws ArgumentError validpool(1:5, Int, 1, 0, 1, true)

    sa1 = ScaledArray(p, 2, usepool=false)
    @test sa1.refs == [1, repeat(1:4, inner=2)..., 4]
    @test sa1.pool == 2:2:9
    @test sa1 != sa

    sa2 = ScaledArray(sa1, 2, 1, 10)
    @test sa2.refs == [1, repeat(1:2:7, inner=2)..., 7]
    @test sa2.pool == 2:10
    @test sa2 == sa1

    sa3 = ScaledArray(sa2, reftype=Int16, start=0, stop=8, usepool=false)
    @test sa3.refs == [3, repeat(3:2:9, inner=2)..., 9]
    @test eltype(sa3.refs) == Int16
    @test sa3.pool == 0:8
    @test sa3 == sa2

    sa4 = ScaledArray(sa3, start=2, stop=8, usepool=false)
    @test sa4.refs == [1, repeat(1:2:7, inner=2)..., 7]
    @test sa4.pool == 2:8
    @test sa4 == sa3

    sa5 = ScaledArray(sa4, usepool=false)
    @test sa5.refs == sa4.refs
    @test sa5.refs !== sa4.refs

    sa6 = ScaledArray([missing, 1, 2], 1)
    @test eltype(sa6) == Union{Int,Missing}
    @test sa6.refs == 0:2
    @test eltype(sa6.refs) == Int32
    @test ismissing(sa6[1])

    @test_throws ArgumentError ScaledArray(sa5, stop=7, usepool=false)
    @test_throws ArgumentError ScaledArray(RefArray(1:3), 1.0:0.5:3.0, Dict{Int,Int}())

    @test similar(sa, 4) == ScaledArray(RefArray(ones(Int, 4)), sa.pool, Dict{Date,Int}())
    @test similar(sa, (4,)) == similar(view(sa, 1:5), 4) == similar(sa, 4)

    xs = Date(5):Year(-1):Date(1)
    sa1 = align(xs, sa)
    @test sa1.pool == pool
    xs = Date(0):Year(1):Date(4)
    @test_throws ArgumentError align(xs, sa)
    sa1 = ScaledArray(sa, Year(-1))
    xs = Date(2):Year(1):Date(6)
    @test_throws ArgumentError align(xs, sa1)

    @test refarray(sa) === sa.refs
    @test refvalue(sa, 1) == Date(1)
    @test refpool(sa) === sa.pool
    @test invrefpool(sa) === sa.invpool

    @test refarray(ssa) == view(sa.refs, 3:5)
    @test refvalue(ssa, 1) == Date(1)
    @test refpool(ssa) === sa.pool
    @test invrefpool(ssa) === sa.invpool
    @test ssa == sa[3:5]

    @test sa[1] == Date(1)
    @test sa[1:2] == sa[[1,2]] == sa[(1:10).<3] == Date.(1:2)
    @test sa[1:2] isa ScaledArray
    @test ssa[1:2] isa ScaledArray

    sa[1] = Date(10)
    @test sa[1] == Date(10)
    sa[1:2] .= Date(100)
    @test sa[1:2] == [Date(100), Date(100)]
    @test last(sa.pool) == Date(100)
    @test_throws MethodError sa[1] = missing
    sa = allowmissing(sa)
    sa[1] = missing
    @test ismissing(sa[1])
    @test_throws ArgumentError disallowmissing(sa)
    sa[1] = Date(1)
    sa = disallowmissing(sa)
    @test eltype(sa) == Date
end

@testset "scaledlabel" begin
    x = Date.(10:-2:0)
    pl = Date(-1):Year(1):Date(11)
    refs, pool, invpool = scaledlabel(x, Year(1), start=Date(-1), stop=Date(11))
    @test refs == 12:-2:2
    @test eltype(refs) == Int32
    @test pool == pl

    pl = Date(0):Year(2):Date(10)
    ip = Dict(p=>i for (i, p) in enumerate(pl))
    refs, pool, invpool = scaledlabel(x, Year(2), Int16)
    @test refs == 6:-1:1
    @test eltype(refs) == Int16
    @test pool == pl
    @test invpool == ip
    @test valtype(invpool) == Int16

    x = ScaledArray(RefArray(refs), Date(0):Year(2):Date(10), invpool)
    refs1, pool1, invpool1 = scaledlabel(x, Year(2))
    @test refs1 == x.refs && refs1 !== x.refs
    @test eltype(refs1) == Int32
    @test pool1 == x.pool
    @test invpool1 == x.invpool && invpool1 !== x.invpool

    refs1, pool1, invpool1 = scaledlabel(x, Year(2), Int16)
    @test refs1 == 6:-1:1
    @test eltype(refs1) == Int16
    @test pool1 == pl
    @test invpool1 == ip
    @test valtype(invpool1) == Int16

    refs1, pool1, invpool1 = scaledlabel(x, Year(1), start=Date(-1), stop=Date(11))
    @test refs1 == 12:-2:2
    @test pool1 == Date(-1):Year(1):Date(11)

    refs, pool, invpool = scaledlabel(1.0:200.0, 1, Int8, Int)
    @test refs == 1:200
    @test eltype(refs) == Int16
    @test pool == 1:200
    @test eltype(pool) == Int
    @test invpool == Dict(i=>Int(i) for i in 1.0:200.0)

    refs, pool, invpool = scaledlabel(1:typemax(Int16), 1, Int8)
    @test refs == 1:typemax(Int16)
    @test eltype(refs) == Int16
    @test valtype(invpool) == Int16

    refs, pool, invpool = scaledlabel([missing, 1, 2], 1)
    @test refs == 0:2
    @test eltype(refs) == Int32
end


text:
export resizable,
       draggable,
       sortable

@api resizable => (Resizable <: Behavior) begin
    arg(tile::Tile)
    kwarg(name::Symbol=:_resizer)
end

render(r::Resizable) =
    render(r.tile)

immutable Size
    proportions::(Float64, Float64)
end

@api draggable => (Draggable <: Behavior) begin
    arg(container::Tile)
    curry(tile::Tile)
    kwarg(name::Symbol=:_draggable)
end

immutable Drag
    fraction::(Float64, Float64)
end

@api sortable => (Sortable <: Behavior) begin
    arg(tile::TileList)
    kwarg(name::Symbol=:_sortable)
end

immutable PositionSwap
    swap::(Int, Int)
end


text:
import brainflow

# enable all possible logs from all three libs
brainflow.enable_dev_brainflow_logger(Integer(brainflow.BOARD_CONTROLLER))
brainflow.enable_dev_brainflow_logger(Integer(brainflow.DATA_HANDLER))
brainflow.enable_dev_brainflow_logger(Integer(brainflow.ML_MODULE))

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)
sampling_rate = brainflow.get_sampling_rate(Integer(brainflow.SYNTHETIC_BOARD))
nfft = brainflow.get_nearest_power_of_two(sampling_rate)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_board_data(board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(Integer(brainflow.SYNTHETIC_BOARD))

bands = brainflow.get_avg_band_powers(data, eeg_channels, sampling_rate, true)
feature_vector = vcat(bands[1], bands[2])

# calc concentration
model_params = brainflow.BrainFlowModelParams(Integer(brainflow.CONCENTRATION), Integer(brainflow.KNN))
concentration = brainflow.MLModel(model_params)
brainflow.prepare(concentration)
print(brainflow.predict(feature_vector, concentration))
brainflow.release(concentration)

# calc relaxation
model_params = brainflow.BrainFlowModelParams(Integer(brainflow.RELAXATION), Integer(brainflow.REGRESSION))
relaxation = brainflow.MLModel(model_params)
brainflow.prepare(relaxation)
print(brainflow.predict(feature_vector, relaxation))
brainflow.release(relaxation)

max_stars_repo_name: Cognixion-inc/brainflow

max_issues_repo_name: Cognixion-inc/brainflow

max_forks_repo_name: Cognixion-inc/brainflow


text:
export sparsePCE

# Compute a sparse basis for the provided pce model and parameters
# Parameters:
#   * op: The full candidate orthogonal basis
function sparsePCE(mop::MultiOrthoPoly, modelFun::Function; Q²tgt=.999, pMax=10, jMax=5)
    # Parameter specification
    pMax = min(mop.deg, pMax)   # pMax can be at most size of given full basis
    dim = length(mop.uni)       # basis dimension
    jMax = min(dim, jMax)       # Maximum interactions is basis dimension at max
    COND = 1e4                  # Maximum allowed matrix condition number (see Blatman2010)
    α = 0.001                   # Tuning parameter (see Blatman2010)
    ϵ = α * (1 - Q²tgt)         # Error threshold of coefficients

    # Initialize basis and determination coefficients
    R² = 0.0    # coefficient of determination
    Q² = 0.0    # leave-one out coefficient determination
    ind = mop.ind # index of complete basis
    pce = []    # pce coefficients, array of arrays
    p = 0       # max degree
    Ap = []

    # 1. Build initial ED, compute Y
    sampleSize = pMax * 20 # TODO: How to determine?
    X = sampleMeasure(sampleSize, mop)
    Y = [modelFun(x...) for x in eachrow(X)] # This is the most expensive part

    restart = true
    # Outer loop: Iterate on experimental design
    while restart
        restart = false

        # 2. Initialize loop
        p = 0
        Ap = [ind[1,:]]    # Set of basis indices for current degree, ind[1] is zero element (array of arrays)
        Φ = [ evaluate(j, X, mop) for j in Ap ] # FIX create matrix directly
        Φ = reduce(hcat, Φ)
        pce = leastSquares(Φ, Y)
        R² = 1 - empError(Y, Φ, pce)
        Q² = 1 - looError(Y, Φ, pce)
        println("** Initialization:")
        println("R²: $R²")
        println("Q²: $Q²")

        # Main loop: Iterate max degree p
        while Q² ≤ Q²tgt && p ≤ pMax && !restart
            p += 1  # current max degree
            candidates = [el for el in eachrow(ind) if sum(el) == p] # only inices with total degree p

            # Iterate number of interactions j
            jM = min(p, jMax) # p limits #interactions
            j = 0   # number of allowed interactions
            while Q² ≤ Q²tgt && j < jM && !restart
                j += 1
                J = [] # Temporary store potential new basis elements
                cands = filter(el -> count(!iszero, el) == j, candidates) # only indices with interaction order j
                println("Candidates for p = $p, j = $j: ", cands)

                # Forward step: compute R² for all candidate terms and keep the relevant ones
                for a in cands
                    A = Ap ∪ [a]
                    Φ = [ evaluate(j, X, mop) for j in A ] # FIX create matrix directly
                    Φ = reduce(hcat, Φ)
                    pce = leastSquares(Φ, Y)
                    R²new = 1 - empError(Y, Φ, pce) # Only need R² error here "due to more efficiency" (Blatman 2010)
                    # println("R²new (p=$p, j=$j): ", R²new)

                    # Compute accuracy gain. If high enough, add cadidate polynomial
                    ΔR² = R²new - R²
                    # println("ΔR²: $ΔR² \t ϵ: $ϵ")
                    if ΔR² > ϵ
                        J = J ∪ [a]
                    end
                end

                # FUTURE: Sort Jp and ΦJ according to ΔR² <- why?
                # J = sort(J)
                # R = []

                # Conditioning Check: If resulting enriched basis does not yield a well-conditioned moments matrix, we have to restart
                Ap_new = Ap ∪ J
                println("Candidates after forward step (p=$p, j=$j): $Ap_new")

                Φ = [ evaluate(j, X, mop) for j in Ap_new]
                Φ = reduce(hcat, Φ)
                # Check Φ
                if cond(Φ) > COND
                    # Increase experimental design and restart computations
                    restart = true
                    k = 3 # rescale factor according to Blatman
                    println("Moments matrix is ill-conditioned. Restart computation with new sample size: $(k * sampleSize)) (Old size: $sampleSize)")
                    sampleSize *= k   # TODO: build properly
                    X = sampleMeasure(sampleSize, mop)
                    Y = [modelFun(x...) for x in eachrow(X)] # TODO: Reuse old ED data, this part is very expensive!
                else
                #     # If conditioning is okay, update accuracy R² for backward step
                    # Φ = [ evaluate(j, X[i], op) for i = 1:sampleSize, j in Ap_new]
                    pce = leastSquares(Φ, Y)
                    R² = 1 - empError(Y, Φ, pce)
                end

                # Backward step: Remove candidate polynomials one by one and compute the effect on Q²
                if !restart
                    println("\n ** Backward Step - Ap_new: ", Ap_new)
                    Del = [] # polynomials to be discarded
                    for a in Ap_new
                        A = filter(e->e≠a,Ap_new)
                        # New candiadte basis -> compute new determination coefficients
                        Φ = [ evaluate(j, X, mop) for j in A]
                        Φ = reduce(hcat, Φ)
                        pce = leastSquares(Φ, Y)
                        R²new = 1 - empError(Y, Φ, pce)

                        # If decrease in accuracy is too small, throw polynomial away
                        ΔR² = R² - R²new
                        if ΔR² ≤ ϵ
                            Del = Del ∪ [a]
                            println("throw away a = ", a)
                        end
                    end

                    # Update basis and compute errors for next iteration
                    Ap = filter(e->e∉Del, Ap_new)
                    Φ = [ evaluate(j, X, mop) for j in Ap]
                    Φ = reduce(hcat, Φ)
                    pce = leastSquares(Φ, Y)
                    println("pce: ", pce)
                    R² = 1 - empError(Y, Φ, pce)
                    Q² = 1 - looError(Y, Φ, pce)
                    println("Q² (p=$p, j=$j): ", Q²)
                end

            end
            println()
        end

    end

    if Q² < Q²tgt
        println("Computation reached max degree $pMax. However, accuracy is below target with Q² = $(Q²). \n")
    else
        println("Computation reached target accuracy with Q² = $(Q²) and degree $p \n")
    end

    return pce, Ap, p, R², Q²
end

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos


text:
function show{T<:ArbFloat}(io::IO, x::T)
    if isexact(x)
      s = string(midpoint(x))
    else
      s = string(x)
    end
    print(io, s)
end
show{T<:ArbFloat}(x::T) = show(STDOUT, x)

function showsmall{T<:ArbFloat}(io::IO, x::T)
    s = stringsmall(x)
    print(io, s)
end
showsmall{T<:ArbFloat}(x::T) = showsmall(STDOUT, x)

function showsmall_pm{T<:ArbFloat}(io::IO, x::T)
    s = stringsmall_pm(x)
    print(io, s)
end
showsmall_pm{T<:ArbFloat}(x::T) = showsmall_pm(STDOUT, x)

function showcompact{T<:ArbFloat}(io::IO, x::T)
    s = stringcompact(x)
    print(io, s)
end
showcompact{T<:ArbFloat}(x::T) = showcompact(STDOUT, x)

function showcompact_pm{T<:ArbFloat}(io::IO, x::T)
    s = stringcompact_pm(x)
    print(io, s)
end
showcompact_pm{T<:ArbFloat}(x::T) = showcompact_pm(STDOUT, x)

function showmedium{T<:ArbFloat}(io::IO, x::T)
    s = stringmedium(x)
    print(io, s)
end
showmedium{T<:ArbFloat}(x::T) = showmedium(STDOUT, x)

function show_pm{T<:ArbFloat}(io::IO, x::T)
    s = stringmedium_pm(x)
    print(io, s)
end
show_pm{T<:ArbFloat}(x::T) = show_pm(STDOUT, x)

function showlarge{T<:ArbFloat}(io::IO, x::T)
    s = stringlarge(x)
    print(io, s)
end
showlarge{T<:ArbFloat}(x::T) = showlarge(STDOUT, x)

function showlarge_pm{T<:ArbFloat}(io::IO, x::T)
    s = stringlarge_pm(x)
    print(io, s)
end
showlarge_pm{T<:ArbFloat}(x::T) = showlarge_pm(STDOUT, x)

function showall{T<:ArbFloat}(io::IO, x::T)
    s = stringall(x)
    print(io, s)
end
showall{T<:ArbFloat}(x::T) = showall(STDOUT, x)

function showall_pm{T<:ArbFloat}(io::IO, x::T)
    s = stringall_pm(x)
    print(io, s)
end
showall_pm{T<:ArbFloat}(x::T) = showall_pm(STDOUT, x)

function showallcompact{T<:ArbFloat}(io::IO, x::T)
    s = stringallcompact(x)
    print(io, s)
end
showallcompact{T<:ArbFloat}(x::T) = showallcompact(STDOUT, x)

function showsmart{T<:ArbFloat}(io::IO, x::T)
    s = smartstring(x)
    print(io, s)
end
# showsmart is not a Base show function, it needs explict version without io parameter
showsmart{T<:ArbFloat}(x::T) = showsmart(STDOUT, x)

function showmany{T<:ArbFloat,N}(io::IO, x::NTuple{N,T}, stringformer::Function)
    if N==0
       print(io,"()")
       return nothing
    elseif N==1
       print(io,string("( ",stringformer(x[1]),", )"))
       return nothing
    end

    ss = Vector{String}(N)
    for i in 1:N
      ss[i] = stringformer(x[i])
    end

    println(io,string("( ", ss[1], ","));
    for s in ss[2:end-1]
      println(io, string("  ", s, ","))
    end
    println(io,string("  ", ss[end], " )"))
end

showmany{T<:ArbFloat,N}(x::NTuple{N,T}, stringformer::Function) =
   showmany(STDOUT,x,stringformer)

function showmany{T<:ArbFloat}(io::IO, x::Vector{T}, stringformer::Function)
    n = length(x)

    if n==0
       print(io,"[]")
       return nothing
    elseif n==1
       print(io,string("[ ",stringformer(x[1])," ]"))
       return nothing
    end

    ss = Vector{String}(n)
    for i in 1:n
      ss[i] = stringformer(x[i])
    end

    println(io,string("[ ", ss[1], ","));
    for s in ss[2:end-1]
      println(io, string("  ", s, ","))
    end
    println(io,string("  ", ss[end], " ]"))
end

showmany{T<:ArbFloat}(x::Vector{T}, stringformer::Function) =
    showmany(STDOUT,x,stringformer)

for (F,S) in [(:showsmall, :stringsmall), (:showcompact, :stringcompact),
              (:show, :stringmedium), (:showlarge, :stringlarge),
              (:showall, :stringall), (:showsmart, :smartstring)]
  @eval begin
     ($F){P,N}(io::IO, x::NTuple{N,ArbFloat{P}}) = showmany(io, x, $S)
     ($F){P,N}(x::NTuple{N,ArbFloat{P}}) = showmany(STDOUT, x, $S)
     ($F){P}(io::IO, x::Vector{ArbFloat{P}}) = showmany(io, x, $S)
     ($F){P}(x::Vector{ArbFloat{P}}) = showmany(STDOUT, x, $S)
     ($F){P,N}(io::IO, x::Vararg{ArbFloat{P},N}) = showmany(io, x, $S)
     ($F){P,N}(x::Vararg{ArbFloat{P},N}) = showmany(STDOUT, x, $S)
  end
end

for (F,S) in [(:showsmall_pm, :stringsmall_pm),
              (:showcompact_pm, :stringcompact_pm),
              (:show_pm, :stringmedium_pm),
              (:showlarge_pm, :stringlarge_pm),
              (:showall_pm, :stringall_pm)]
  @eval begin
     ($F){P,N}(io::IO, x::NTuple{N,ArbFloat{P}}) = showmany(io, x, $S)
     ($F){P,N}(x::NTuple{N,ArbFloat{P}}) = showmany(STDOUT, x, $S)
     ($F){P}(io::IO, x::Vector{ArbFloat{P}}) = showmany(io, x, $S)
     ($F){P}(x::Vector{ArbFloat{P}}) = showmany(STDOUT, x, $S)
     ($F){P,N}(io::IO, x::Vararg{ArbFloat{P},N}) = showmany(io, x, $S)
     ($F){P,N}(x::Vararg{ArbFloat{P},N}) = showmany(STDOUT, x, $S)
  end
end


text:
function custom_martingale(stake::Float64, res::Bool)
    if res == false
        return stake+(stake*0.15)
    else
        return stake
end

max_stars_repo_name: jkotra/LuckSimulations

max_issues_repo_name: jkotra/LuckSimulations

max_forks_repo_name: jkotra/LuckSimulations


text:
function runfig5()
    p = Params(Ln = 2000, W = 2000, Ls = 20, Ws = 20, scale = 40, λ = 5, α = 0,
    Δ = 1, d = 0, τ = 1)

    p = reconstruct(p, μN = 0.000001)
    sp = spectrumsweepb(p, 0:.25:5, true)
    savespectrum("fig5", p, sp[1], sp[2])

    p = reconstruct(p, μN = 0.6)
    sp = spectrumsweepb(p, 0:.25:5, true)
    savespectrum("fig5", p, sp[1], sp[2])

    p = reconstruct(p, μN = 1.192)
    sp = spectrumsweepb(p, 0:.25:5, true)
    savespectrum("fig5", p, sp[1], sp[2])

    p = reconstruct(p, EZ = SA[0, 2, 0], μN = 0.6)
    sp = spectrumsweepα(p, 0:.25:4, true)
    savespectrum("fig5", p, sp[1], sp[2])

    p = reconstruct(p, EZ = SA[0, 4, 0], μN = 1.192)
    sp = spectrumsweepα(p, 0:.25:4, true)
    savespectrum("fig5", p, sp[1], sp[2])
end

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
module CommunalHelperCassetteMoveBlock
using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/CassetteMoveBlock" CassetteMoveBlock(x::Integer,
                                                              y::Integer,
                                                              width::Integer=Maple.defaultBlockWidth,
                                                              height::Integer=Maple.defaultBlockHeight,
                                                              direction::String="Right",
                                                              fast::Bool=false,
                                                              index::Integer=0,
                                                              tempo::Number=1.0)

const colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Bright Sun" => 2,
    "Malachite" => 3
)

const colors = Dict{Int, Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
	2 => (252, 220, 58, 255) ./ 255,
	3 => (56, 224, 78, 255) ./ 255
)
const defaultColor = (73, 170, 240, 255) ./ 255

const placements = Ahorn.PlacementDict(
    "Cassette Move Block ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteMoveBlock,
        "rectangle",
        Dict{String, Any}(
            "index" => index,
        )
    ) for (color, index) in colorNames
)

Ahorn.editingOptions(entity::CassetteMoveBlock) = Dict{String, Any}(
    "direction" => Maple.move_block_directions
)
Ahorn.minimumSize(entity::CassetteMoveBlock) = 16, 16
Ahorn.resizable(entity::CassetteMoveBlock) = true, true

Ahorn.selection(entity::CassetteMoveBlock) = Ahorn.getEntityRectangle(entity)

const arrows = Dict{String, String}(
    "up" => "objects/CommunalHelper/cassetteMoveBlock/arrow02",
    "left" => "objects/CommunalHelper/cassetteMoveBlock/arrow04",
    "right" => "objects/CommunalHelper/cassetteMoveBlock/arrow00",
    "down" => "objects/CommunalHelper/cassetteMoveBlock/arrow06",
)

const block = "objects/cassetteblock/solid"

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CassetteMoveBlock, room::Maple.Room)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    index = Int(get(entity.data, "index", 0))
    color = get(colors, index, defaultColor)

    for i in 1:tileWidth, j in 1:tileHeight
        tx = (i == 1) ? 0 : ((i == tileWidth) ? 16 : 8)
        ty = (j == 1) ? 0 : ((j == tileHeight) ? 16 : 8)

        Ahorn.drawImage(ctx, block, (i - 1) * 8, (j - 1) * 8, tx, ty, 8, 8, tint=color)
    end

    direction = lowercase(get(entity.data, "direction", "up"))
    arrowSprite = Ahorn.getSprite(arrows[lowercase(direction)], "Gameplay")

    Ahorn.drawImage(ctx, arrowSprite, div(width - arrowSprite.width, 2), div(height - arrowSprite.height, 2), tint=color)
end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
function _make_test_grid(args...)
    return reshape(Iterators.flatten(
                            reduce(vcat, Iterators.product(args...))
                            ) |> collect,
                    length(args), *(length.(args)...)) |> permutedims
end

function _make_test_df(args...; vars)
    DataFrame(_make_test_grid(args...), vars)
end

function _make_test_rows(args...)
    RowVecs(_make_test_grid(args...))
end

function getrank(v; flattenzeros=true)
    r = invperm(sortperm(v))
    if flattenzeros
        z = findall(iszero, v)
        r[z] .= 1
    end
    return r
end

function invnormaltransform(v; μ=0, σ=1, c=3/8, flattenzeros=true)
    rank = getrank(v; flattenzeros)
    return [norminvcdf(μ, σ, (x - c) / (length(v) - 2c + 1)) for x in rank]
end

function _df_output(df, args)
    if !isnothing(args["output"])
        @info "Writing to $(args["output"])"
        out = args["output"]
        delim = endswith(out, "csv") ? ',' :
              endswith(out, "tsv") ? '\t' : error("--output arg must be '.tsv' or '.csv'")
        CSV.write(expanduser(args["output"]), df; delim)
    else
        @show df
    end
end

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
using RvSpectraKitLearn

# Likely need to update this to point to directory containing the spectra you wish to analyze
include(joinpath(Pkg.dir("RvSpectraKitLearn"),"path_to_spectra.jl"))
cd(path_to_spectra)
(datafiles_planet, phases_planet) = get_filenames();
(lambda_planet, obs_planet_noisefree) = read_filelist(datafiles_planet);
planet_amplitude = 10.0
rvs_planet = -planet_amplitude*cos(phases_planet*2pi)

path_to_spectra = joinpath(Pkg.dir(),"..","..","Code","XavierSpotData")
cd(path_to_spectra)

function get_filenames_multi_lats(; num_max_files::Integer = 0)
  files = readdir()
  datafiles = String[]
  strengths = Float64[]
  inclinations = Float64[]
  latitudes = Float64[]
  resolutions = Float64[]
  samplings = Float64[]
  phases = Float64[]
  for f in files
    m = match(r"(\w+)_(\d+)pc_incl(\d+)_lat(\d+)_reso(\d+)k_(-?\d+\.?\d*)pxsampl_(-?\d+\.?\d*)\.csv",f)
    if m==nothing continue end
    effect_name = m.captures[1]
    resolution = parse(m.captures[5])
    if resolution> 149
    push!(strengths,parse(m.captures[2]))
	  push!(inclinations,parse(m.captures[3]))
    push!(latitudes,parse(m.captures[4]))
    push!(resolutions,parse(m.captures[5]))
    push!(samplings,parse(m.captures[6]))
    push!(phases,parse(m.captures[7]))
    push!(datafiles,f)
    end
    if num_max_files>0 && length(datafiles)==num_max_files break end
  end
  @assert(length(datafiles)>0)
  return (datafiles,phases)
  #idx = sortperm(phases)
  #return (datafiles[idx],phases[idx])
end

(datafiles_spot, phases_spot) = get_filenames_multi_lats(); # num_max_files=4);
(lambda_spot, obs_spot_noisefree) = read_filelist(datafiles_spot);
rvs_spot = zeros(phases_spot)

rvs_all = [rvs_spot, rvs_planet]
@assert(all(lambda_spot .== lambda_planet))
obs_noisefree = Array(Float64,(size(obs_planet_noisefree,1),size(obs_planet_noisefree,2)+size(obs_spot_noisefree,2)))
obs_noisefree[:,1:size(obs_spot_noisefree,2)] = obs_spot_noisefree
obs_noisefree[:,1+size(obs_spot_noisefree,2):end] = obs_planet_noisefree
obs_noisefree

function find_important_lamba{T}(lambda::Array{T,1}, obs::Array{T,2}, rvs::Array{T,1};
                                 chunk_size::Integer=1500, frac_to_keep = 0.1 )
  Y = reshape(rvs,(length(rvs),1))
 beta = ones(length(lambda)) # coefficients
 new_beta = copy(beta)
  idx_active = find(beta)
  num_chunks = ceil(Int64,length(idx_active)/chunk_size)
  for c in 1:num_chunks
    idx_min = 1+(c-1)*chunk_size
    idx_max = min(c*chunk_size,length(idx_active))
    lambda_range = idx_active[idx_min:idx_max]
    #println("idx=",idx_min,"-",idx_max,"  active=",idx_active[idx_min],"-",idx_active[idx_max]," lambda=",lambda[idx_active[idx_min]],"-",lambda[idx_active[idx_max]])
    X = obs[lambda_range,:]'
    (gamma_list, Z_list)  = fit_arp_slr_2d(X,Y,stepsize=1.01,gamma_init=0.001, max_itterations=1000, min_active=floor(Int64,size(X,2)*size(Y,2)*frac_to_keep) )
    new_beta[lambda_range] = vec(Z_list[size(Z_list,1),:,1])
    #println("size(Z_list)=", size(Z_list,1), " nactive=",length(find(new_beta[lambda_range])))
  end
  beta[:] = new_beta
 return new_beta
end

beta = find_important_lamba(lambda_spot,obs_noisefree,rvs_all, chunk_size=200, frac_to_keep = [0.1] )

#Pkg.add("JLD")  # If you don't already have JLD installed
using JLD
save("important_lambda_0.1.jld", "frac_to_keep", [0.1], "beta", beta)

#=
# Plot which wavelengths are being used from the last itteration of the algorithm regularization path algorithm
=#
using PyPlot
lambda = lambda_planet
obs = obs_noisefree
lambda_range = 1:length(lambda_planet)
 zmax = maximum(abs(beta[lambda_range]))
  ymax = maximum(vec(obs[lambda_range,1]))
  ymin = minimum(vec(obs[lambda_range,1]))
  yrange = ymax-ymin
  #plot(lambda,vec(obs[:,1]),"r.")
  min_lambda = lambda[lambda_range][1]
  max_lambda = lambda[lambda_range][end]
  idx_around = find(x->min_lambda<=x<=max_lambda,lambda)
  plot(lambda[idx_around],vec(obs[idx_around,1]),"r-")
  plot(lambda[lambda_range],ymin+yrange/2*(1+1/zmax*(beta[lambda_range])),"b+")
  plot(lambda[lambda_range],ymin+yrange/2*(ones(length(lambda_range))),"b-")

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
# TensorOperations.jl interface.
#

contract!(α,
          A::DualTensor{N, T}, conjA::Symbol,
          B::DualTensor{N, T}, conjB::Symbol,
          β,
          C::DualTensor{N, T},
          oindA::IndexTuple, cindA::IndexTuple,
          oindB::IndexTuple, cindB::IndexTuple,
          indleft::IndexTuple,
          indright::IndexTuple, syms = nothing) where {N, T} =
    contract!(α,
              A, conjA, B, conjB,
              β, C,
              oindA, cindA,
              oindB, cindB,
              (indleft..., indright...), syms)

contract!(α,
          A::DualTensor{N, T}, conjA::Symbol,
          B::DualTensor{N, T}, conjB::Symbol,
          β,
          C::DualTensor{N, T},
          oindA::IndexTuple, cindA::IndexTuple,
          oindB::IndexTuple, cindB::IndexTuple,
          tindC::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing) where {N, T} = begin

    # Value part.
    contract!(α,
              A.A, conjA, B.A, conjB,
              β, C.A,
              oindA, cindA,
              oindB, cindB,
              tindC, syms)

    # Loop over derivatives.
    # ∂A part.
    for i=1:N
        contract!(α,
                  A.∂[i], conjA, B.A, conjB,
                  β, C.∂[i],
                  oindA, cindA,
                  oindB, cindB,
                  tindC, syms)
    end
    # ∂B part.
    for i=1:N
        contract!(α,
                  A.A, conjA, B.∂[i], conjB,
                  true, C.∂[i],
                  oindA, cindA,
                  oindB, cindB,
                  tindC, syms)
    end

    C
end


text:
using CSV
using DataFrames
using Ipopt
using JuMP

n_steps = 101 # Setting number of ODE discretisation steps

# Data
println("Reading measurement data...")
data_path = "/media/sf_DPhil_Project/Project07_Parameter Fitting/df_software/sbml2julia/tests/fixtures/0015_objectivePrior/_measurements.tsv"
df = CSV.read(data_path)
insert!(df, 1, (1:length(df[:,1])), :id)

df_by_o = groupby(df, :observableId)
df_by_c = groupby(df, :simulationConditionId)

# Setting variables
t_sim = range(0, stop=maximum(df[!, :time]), length=n_steps)
t_exp = Dict()
m_exp = Dict()
t_sim_to_exp = Dict()
deduplicated_time_idx = Dict()

cond2idx = Dict()
for i_cond in 1:length(keys(df_by_c))
    cond_name = values(keys(df_by_c)[i_cond])[1]
    cond2idx[cond_name] = i_cond
end

for i_obs in 1:length(keys(df_by_o))
    obs_name = values(keys(df_by_o)[i_obs])[1]
    t_exp[obs_name] = Dict()
    m_exp[obs_name] = Dict()
    t_sim_to_exp[obs_name] = Dict()
    deduplicated_time_idx[obs_name] = Dict()
    df_by_o_c = groupby(DataFrame(df_by_o[i_obs]), :simulationConditionId)
    for i_cond in 1:length(keys(df_by_o_c))
        cond_name = values(keys(df_by_o_c)[i_cond])[1]
        cond_idx = cond2idx[cond_name]
        t_exp[obs_name][cond_idx] = df_by_o_c[i_cond][!, :time]
        m_exp[obs_name][cond_idx] = df_by_o_c[i_cond][!, :measurement]

        t = df_by_o_c[i_cond][!, :time]
        tmp = []
        for i in 1:length(t)
            idx = argmin(abs.(t[i] .- t_sim))
            append!(tmp, idx)
        end
        t_sim_to_exp[obs_name][cond_idx] = tmp

        deduplicated_time_idx[obs_name][cond_idx] = []
        idx = 0
        prev = "a"
        for t in df_by_o_c[i_cond][!, :time]
            if t != prev
                idx = idx + 1
            end
            push!(deduplicated_time_idx[obs_name][cond_idx], idx)
            prev = t
        end
    end
end

obs2conds = Dict()
for obs_dict in t_exp
    obs_name = obs_dict[1]
    obs2conds[obs_name] = []
    for cond_dict in obs_dict[2]
        cond_idx = cond_dict[1]
        push!(obs2conds[obs_name], cond_idx)
    end
    sort!(obs2conds[obs_name])
end

results = Dict()
results["objective_value"] = Dict()
results["parameters"] = Dict()
results["species"] = Dict()
results["observables"] = Dict()

j_to_cond_par = [1, 2]
cond_without_preequ = [1]
cond_with_preequ = [2]
preequ = [3]

i_start = 1

    species_dict = Dict()
    obs_dict = Dict()

    m = Model(with_optimizer(Ipopt.Optimizer))

    # Define global parameters
    println("Defining global parameters...")
    @variable(m, 0.0 <= a0 <= 10.0, start=0.0+(10.0-(0.0))*rand(Float64))
    @variable(m, 0.0 <= b0 <= 10.0, start=0.0+(10.0-(0.0))*rand(Float64))
    @variable(m, 0.0 <= k1_free <= 10.0, start=0.0+(10.0-(0.0))*rand(Float64))
    @variable(m, 0.0 <= k2 <= 10.0, start=0.0+(10.0-(0.0))*rand(Float64))
    @variable(m, 0.005 <= noise_A1 <= 0.1, start=0.005+(0.1-(0.005))*rand(Float64))
    @variable(m, 0.01 <= noise_A2 <= 0.2, start=0.01+(0.2-(0.01))*rand(Float64))
    @variable(m, 0.01 <= noise_B <= 0.2, start=0.01+(0.2-(0.01))*rand(Float64))
    @variable(m, 0.0 <= offset_B <= 5.0, start=0.0+(5.0-(0.0))*rand(Float64))
    @variable(m, 0.0 <= scaling_B <= 10.0, start=0.0+(10.0-(0.0))*rand(Float64))

    # Define condition-specific parameters
    println("Defining condition-specific parameters...")
    @variable(m, k1[1:3])
    @constraint(m, k1[1] == k1_free)
    @constraint(m, k1[2] == k1_free)
    @constraint(m, k1[3] == 0.1)

    # Define overrides
    # Define observable overrides
    println("Defining observableParameter overrides...")

    # Define noise overrides
    println("Defining noiseParameter overrides...")

    # Model compartments
    println("Defining compartments...")
    @variable(m, compartment == 1.0, start=1.0)

    # Model species
    println("Defining species...")
    species_dict["A"] = @variable(m, 0.0 <= A[j in 1:2, k in 1:(length(t_sim)+1)] <= 20.038222743)
    species_dict["B"] = @variable(m, 0.0 <= B[j in 1:2, k in 1:(length(t_sim)+1)] <= 20.038222743)

    # Model initial assignments
    println("Defining initial assignments...")
    @constraint(m, [j in 1:1], A[cond_without_preequ[j],1] == a0)
    @constraint(m, [j in 1:1], B[cond_without_preequ[j],1] == b0)

    # Model ODEs
    println("Defining ODEs...")
    println("A")
    @NLconstraint(m, [j in 1:2, k in 1:length(t_sim)-1],
        A[j, k+1] == A[j, k] + ( -1.0*( compartment * k1[j_to_cond_par[j]] * A[j, k+1] ) +1.0*( compartment * k2 * B[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("B")
    @NLconstraint(m, [j in 1:2, k in 1:length(t_sim)-1],
        B[j, k+1] == B[j, k] + ( +1.0*( compartment * k1[j_to_cond_par[j]] * A[j, k+1] ) -1.0*( compartment * k2 * B[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )

    # Pre-equilibration constraints
    println("Defining pre-equilibration constraints...")
    @constraint(m, [j in cond_without_preequ], A[j, length(t_sim)+1] == 0) # Dummy preequilibration for these conditions
    println("A")
    @NLconstraint(m, [j in 1:1],
        -1.0*( compartment * k1[preequ[j]] * A[cond_with_preequ[j], length(t_sim)+1] ) +1.0*( compartment * k2 * B[cond_with_preequ[j], length(t_sim)+1] )  == 0 )
    @constraint(m, [j in 1:1], A[cond_with_preequ[j], length(t_sim)+1] == A[cond_with_preequ[j], 1])
    @constraint(m, [j in cond_without_preequ], B[j, length(t_sim)+1] == 0) # Dummy preequilibration for these conditions
    println("B")
    @NLconstraint(m, [j in 1:1],
        +1.0*( compartment * k1[preequ[j]] * A[cond_with_preequ[j], length(t_sim)+1] ) -1.0*( compartment * k2 * B[cond_with_preequ[j], length(t_sim)+1] )  == 0 )
    @constraint(m, [j in 1:1], B[cond_with_preequ[j], length(t_sim)+1] == B[cond_with_preequ[j], 1])

    # Define observables
    println("Defining observables...")
    obs_dict["obs_a"] = @variable(m, -0.4096061150000001 <= obs_a[j in obs2conds["obs_a"], k in 1:length(t_exp["obs_a"][j])] <= 1.8169054180000002, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_a"], k in 1:length(t_exp["obs_a"][j])], obs_a[j, k] == A[j, t_sim_to_exp["obs_a"][j][k]])
    obs_dict["obs_b"] = @variable(m, -0.6546633489999999 <= obs_b[j in obs2conds["obs_b"], k in 1:length(t_exp["obs_b"][j])] <= 3.059816594, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_b"], k in 1:length(t_exp["obs_b"][j])], obs_b[j, k] == offset_B+scaling_B*B[j, t_sim_to_exp["obs_b"][j][k]])

    # Defining objectivePriors
    println("Defining objectivePriors")
    @variable(m, prior_mean[l in 1:4], start=1.)
    @constraint(m, prior_mean[1] == 1)
    @constraint(m, prior_mean[2] == 0)
    @constraint(m, prior_mean[3] == -1.22314355131)
    @constraint(m, prior_mean[4] == -0.51082562376)

    @variable(m, 0 <= prior_std[l in 1:4], start=1.)
    @constraint(m, prior_std[1] == 1)
    @constraint(m, prior_std[2] == 1)
    @constraint(m, prior_std[3] == 1)
    @constraint(m, prior_std[4] == 2)

    @variable(m, par_est[l in 1:4], start=1.)
    @constraint(m, par_est[1] == a0)
    @constraint(m, par_est[2] == b0)
    @constraint(m, par_est[3] == k1_free)
    @constraint(m, par_est[4] == k2)

    normal_priors = [1]
    laplace_priors = [2]
    logNormal_priors = [3]
    logLaplace_priors = [4]
    @variable(m, delta_par[l in [2]])
    @constraint(m, [l in [2]], delta_par[l] >= par_est[l] - prior_mean[l])
    @constraint(m, [l in [2]], delta_par[l] >= prior_mean[l] - par_est[l])

    @variable(m, delta_log_par[l in [4]])
    @NLconstraint(m, [l in [4]], delta_log_par[l] >= log(par_est[l]) - prior_mean[l])
    @NLconstraint(m, [l in [4]], delta_log_par[l] >= prior_mean[l] - log(par_est[l]))

    @variable(m, nlp_normal)
    @NLconstraint(m, nlp_normal == sum( 0.5 * log(2*pi*prior_std[l]^2) + 0.5*((par_est[l]-prior_mean[l])/prior_std[l])^2 for l in normal_priors))
    @variable(m, nlp_laplace)
    @NLconstraint(m, nlp_laplace == sum( log(2*prior_std[l]) + delta_par[l]/prior_std[l] for l in laplace_priors))
    @variable(m, nlp_logNormal)
    @NLconstraint(m, nlp_logNormal == sum( log(par_est[l]*prior_std[l]*sqrt(2*pi)) + 0.5*((log(par_est[l])-prior_mean[l])/prior_std[l])^2 for l in logNormal_priors))
    @variable(m, nlp_logLaplace)
    @NLconstraint(m, nlp_logLaplace == sum( log(2*prior_std[l]*par_est[l]) + delta_log_par[l]/prior_std[l] for l in logLaplace_priors))

    @variable(m, neg_log_priors)
    @constraint(m, neg_log_priors == nlp_laplace + nlp_logLaplace + nlp_logNormal + nlp_normal)

    # Define objective
    println("Defining objective...")
    @NLobjective(m, Min, sum(0.5 * log(2*pi*(( noise_A1 + noise_A2 * obs_a[j, k] ))^2) + 0.5*((obs_a[j, deduplicated_time_idx["obs_a"][j][k]]-m_exp["obs_a"][j][k])/(( noise_A1 + noise_A2 * obs_a[j, k] )))^2 for j in obs2conds["obs_a"] for k in 1:length(t_exp["obs_a"][j]))
        + sum(0.5 * log(2*pi*(( noise_B ))^2) + 0.5*((obs_b[j, deduplicated_time_idx["obs_b"][j][k]]-m_exp["obs_b"][j][k])/(( noise_B )))^2 for j in obs2conds["obs_b"] for k in 1:length(t_exp["obs_b"][j]))
        + neg_log_priors)

    println("Optimizing:")
    optimize!(m)

    println("Transfering results to Python...")
    parameter_names = [a0, b0, k1_free, k2, noise_A1, noise_A2, noise_B, offset_B, scaling_B]
    parameter_values = Dict()
    for p in parameter_names
        if occursin("[", string(p))
            parameter_values[split(string(p[1]), "[")[1]] = JuMP.value.(p)
        else
            parameter_values[string(p)] = JuMP.value.(p)
        end
    end

    species_names = ["A", "B", ]
    species_values = Dict(string(spec)=>Dict(cond_idx=>[value(species_dict[spec][cond_idx, val_idx]) for val_idx=1:length(t_sim)]
        for cond_idx in 1:length(keys(df_by_c))) for spec in species_names)

    observable_values = Dict(obs=>Dict(cond_idx=>[value(obs_dict[obs][cond_idx, val_idx]) for val_idx=1:length(list)]
        for (cond_idx, list) in dict) for (obs, dict) in m_exp)

    objective_val = objective_value(m)

    results["objective_value"][string(i_start)] = objective_val
    results["parameters"][string(i_start)] = parameter_values
    results["species"][string(i_start)] = species_values
    results["observables"][string(i_start)] = observable_values

results

max_stars_repo_name: paulflang/DisFit

max_issues_repo_name: paulflang/DisFit

max_forks_repo_name: paulflang/DisFit


text:
using JLBoost: _fit_tree!, LogitLogLoss, depth_wise, get_leaf_nodes, keeprow_vec, max_depth_stopping_criterion

using RDatasets, DataFrames

loss = LogitLogLoss()

tbl = dataset("datasets", "iris")

tbl[:is_setosa] = tbl.Species .== "setosa"

target = :is_setosa
features = setdiff(Symbol.(names(tbl)), [:is_setosa, :Species])

warm_start = fill(0.0, nrow(tbl))
lambda = 0
gamma = 0
verbose = false

jlt = JLBoostTree(0.0)
tree_growth = depth_wise
stopping_criterion = max_depth_stopping_criterion(6)

# _fit_tree!(loss, tbl, target, features, warm_start,
#     jlt::AbstractJLBoostTree = JLBoostTree(0.0),
#     tree_growth::Function = depth_wise,
#     stopping_criterion::Function = max_depth_stopping_criterion(6);
#     #colsample_bynode = 1, colsample_bylevel = 1,
# 	kwargs...)


text:
module Progress

using Dates
import Parameters: @with_kw
import UUIDs: UUID

import Term: int, textlen, truncate, loop_last, get_file_format, update!
import ..Tprint: tprint, tprintln
import ..Style: apply_style
import ..Consoles:
    console_width,
    hide_cursor,
    show_cursor,
    move_to_line,
    cleartoend,
    change_scroll_region,
    console_height,
    up,
    down,
    erase_line,
    savecursor,
    restorecursor

import ..Renderables: AbstractRenderable
import ..Measures: Measure
import ..Segments: Segment
import ..Colors: RGBColor
import ..Layout: hLine

export ProgressBar,
    ProgressJob, addjob!, start!, stop!, update!, removejob!, with, @track, render

# ---------------------------------------------------------------------------- #
#                               PROGRESS BAR JOB                               #
# ---------------------------------------------------------------------------- #

# ------------------------------- constructors ------------------------------- #
"""
    ProgressJob

Single `job` whose progress we're tracking with a progress bar.

Each progress bar can have multiple jobs running at the same time, and more can
be added/removed at any time. Each `ProgressJob` keeps track of the state (progress)
of each of these jobs. `ProgressJob` is also rendered to create the visual display
of the progress bar.

Arguments:
- `id` specifies the job;s unique id
- `i`: keeps track of the progress (number of completed steps)
- `N`: total number of steps in the task at hand. Optional, set to `nothing` when not known
- `description`: a bit of text describing what the job is
- `columns`: set of `AbstractColumn` used to represent the job's information
- `columns_kwargs`: additional information passed to each column when created (e.g. to set style)
- `width`: width of the progress bar display
- `started`: whether the job was started or not.
- `stopped`: whether the job was stopped (finished job)
- `startime`: time at which job was started
- `stoptime`: time at which job was stopped
- `transient`: if truee the job's visual display disappears when `stopped=true`
"""
mutable struct ProgressJob
    id::Union{Int,UUID}

    i::Int   # keep track of progress
    N::Union{Nothing,Int}

    description::String
    columns::Any
    columns_kwargs::Dict
    width::Int

    started::Bool
    finished::Bool
    startime::Union{Nothing,DateTime}
    stoptime::Union{Nothing,DateTime}
    transient::Bool

    """
        ProgressJob(
            id::Union{Int, UUID},
            N::Union{Int, Nothing},
            description::String,
            columns::Vector{DataType},
            width::Int,
            columns_kwargs::Dict,
            transient::Bool,
        )

    Constructor for a `ProgressJob`.

    See also [`addjob!`](@ref), [`start!`](@ref), [`stop!`](@ref), [`update!`](@ref), [`redender`](@ref)
    """
    function ProgressJob(
        id::Union{Int,UUID},
        N::Union{Int,Nothing},
        description::String,
        columns::Vector{DataType},
        width::Int,
        columns_kwargs::Dict,
        transient::Bool,
    )
        return new(
            id,
            isnothing(N) ? 0 : 1,
            N,
            description,
            columns,
            columns_kwargs,
            width,
            false,
            false,
            nothing,
            nothing,
            transient,
        )
    end
end

Base.show(io::IO, ::MIME"text/plain", job::ProgressJob) =
    print(io, "Progress job $(job.id) \e[2m(started: $(job.started))\e[0m")

"""
    start!(job::ProgressJob)

Start a newly created `ProgressJob`.

When starting a job, take care of creating instance of `AbstractColumns`
to display the job's progress. If `N` is nothing, remove any `ProgressColumn`
"""
function start!(job::ProgressJob)
    job.started && return nothing

    # if the job doesn't have a defined `N`, we can't have a progress column display
    if isnothing(job.N)
        filter!(c -> c != ProgressColumn, job.columns)

        # but we should have a spinner column if there isnt one.
        !any(job.columns .== SpinnerColumn) && push!(job.columns, SpinnerColumn)
    end

    # create columns type instances, passing the appropriate keyword arguments
    csymbol(c) = Symbol(split(string(c), ".")[end])
    function makecol(c)
        return if haskey(job.columns_kwargs, csymbol(c))
            c(job; job.columns_kwargs[csymbol(c)]...)
        else
            c(job)
        end
    end
    job.columns = map(c -> makecol(c), job.columns)

    # if there's a progress column, set its width
    if !isnothing(job.N) && any(map(c -> c isa ProgressColumn, job.columns))
        # get the progress column width
        spaces = length(job.columns) - 1
        colwidths = sum(c -> c.measure.w, job.columns)
        bcol_width = job.width - colwidths - spaces

        # set width
        setwidth!.(job.columns, bcol_width)
    end

    # start job
    job.started = true
    job.startime = now()
    return nothing
end

"""
    update!(job::ProgressJob; i = nothing)

Update a job's progress `i` by setting its value or adding `+1`.
"""
function update!(job::ProgressJob; i = nothing)
    (!isnothing(job.N) && job.i ≥ job.N) && return stop!(job)

    job.i = isnothing(i) ? job.i + 1 : job.i + i
    return nothing
end

"""
    stop!(job::ProgressJob)

Stop a running job.
"""
function stop!(job::ProgressJob)
    job.stoptime = now()
    job.finished = true
    return nothing
end

# ---------------------------------------------------------------------------- #
#                                    COLUMNS                                   #
# ---------------------------------------------------------------------------- #
# load columns types definitions
include("_progress.jl")

# ---------------------------------------------------------------------------- #
#                                  PROGRESS BAR                                #
# ---------------------------------------------------------------------------- #
"""
    RenderStatus

Keep track of rendering information for a `ProgressBar`
"""
Base.@kwdef mutable struct RenderStatus
    rendered::Bool = false
    nlines::Int = 0
    maxnlines::Int = 0
    hline::String = ""
    scrollline::Int = 0
end

# ------------------------------- constructors ------------------------------- #
"""
    ProgressBar

Progressbar Type, stores information required
to render progress bar renderables for each `ProgressJob` assigned.

ProgressBar takes care of the work needed to setup/update progress bar(s)
visualizations. Each individual bar corresponds to a (running) `ProgressJob`
and `ProgressJob` itself is what actually creates the visuals.
Most of the work done by ProgressBar is to handling terminal stuff, move
cursor position, change scrolling regions, clear/update sections etc.

Arguments:
    - `jobs`: vector of `ProgressJob` to assign to the progress bar (more can added later)
    - `width`: width of the visualization, if `expand=false`
    - `columns`: which columns to show? Either a vector of column types or the name of a preset
    - `column_kwargs`: keyword arguments to pass to each column
    - `transient`: if true jobs disappear when done (and the whole pbar when all jobs are done)
    - `colors`: set of 3 `RGBColor` to change the color of the progress bar with progress
    - `Δt`: delay between refreshes of the visualization.
    - `buff`: `IOBuffer` used to render the progress bars
    - `running`: true if the progress bar is active
    - `paused`: false when the bar is running but briefly paused (e.g. to update `jobs`)
    - `task`: references a `Task` for updating the progress bar in parallel
    - `renderstatus`: a `RenderStatus` instance.
"""
mutable struct ProgressBar
    jobs::Vector{ProgressJob}
    width::Int

    columns::Vector{DataType}
    columns_kwargs::Dict

    transient::Bool
    colors::Vector{RGBColor}
    Δt::Float64

    buff::IOBuffer  # will be used to store temporarily re-directed stdout
    running::Bool
    paused::Bool
    task::Union{Task,Nothing}
    renderstatus::Any
end

"""
    ProgressBar(;
        width::Int=88,
        columns::Union{Vector{DataType}, Symbol} = :default,
        columns_kwargs::Dict = Dict(),
        expand::Bool=false,
        transient::Bool = false,
        colors::Vector{RGBColor} = [
            RGBColor("(1, .05, .05)"),
            RGBColor("(.05, .05, 1)"),
            RGBColor("(.05, 1, .05)"),
        ],
        refresh_rate::Int=60,  # FPS of rendering
    )

Create a ProgressBar instance.
"""
function ProgressBar(;
    width::Int = 88,
    columns::Union{Vector{DataType},Symbol} = :default,
    columns_kwargs::Dict = Dict(),
    expand::Bool = false,
    transient::Bool = false,
    colors::Vector{RGBColor} = [
        RGBColor("(1, .05, .05)"),
        RGBColor("(.05, .05, 1)"),
        RGBColor("(.05, 1, .05)"),
    ],
    refresh_rate::Int = 60,  # FPS of rendering
)
    columns = columns isa Symbol ? get_columns(columns) : columns

    # check that width is large enough
    width = expand ? console_width() - 5 : min(max(width, 20), console_width() - 5)

    return ProgressBar(
        Vector{ProgressJob}(),
        width,
        columns,
        columns_kwargs,
        transient,
        colors,
        1 / refresh_rate,
        IOBuffer(),
        false,
        false,
        nothing,
        RenderStatus(),
    )
end

function Base.show(io::IO, ::MIME"text/plain", pbar::ProgressBar)
    return print(io, "Progress bar \e[2m($(length(pbar.jobs)) jobs)\e[0m")
end

# ---------------------------------------------------------------------------- #
#                                    METHODS                                   #
# ---------------------------------------------------------------------------- #
# --------------------------------- edit pbar -------------------------------- #
"""
    addjob!(
            pbar::ProgressBar;
            description::String="Running...",
            N::Union{Int, Nothing}=nothing,
            start::Bool=true,
            transient::Bool=false,
            id=nothing
        )::ProgressJob

Add a new `ProgressJob` to a running `ProgressBar`

See also: [`removejob!`](@ref), [`getjob`](@ref)
"""
function addjob!(
    pbar::ProgressBar;
    description::String = "Running...",
    N::Union{Int,Nothing} = nothing,
    start::Bool = true,
    transient::Bool = false,
    id = nothing,
)::ProgressJob
    pbar.running && print("\n")

    # create Job
    pbar.paused = true
    id = isnothing(id) ? length(pbar.jobs) + 1 : id
    job = ProgressJob(
        id,
        N,
        description,
        pbar.columns,
        pbar.width,
        pbar.columns_kwargs,
        transient,
    )

    # start job
    start && start!(job)
    push!(pbar.jobs, job)
    pbar.paused = false
    return job
end

"""
    removejob!(pbar::ProgressBar, job::ProgressJob)

Remove a `ProgressJob` from a `ProgressBar`.

See also: [`addjob!`](@ref), [`getjob`](@ref)
"""
function removejob!(pbar::ProgressBar, job::ProgressJob)
    pbar.paused = true
    stop!(job)
    deleteat!(pbar.jobs, findfirst(j -> j.id == job.id, pbar.jobs))
    return pbar.paused = false
end

"""
    getjob(pbar::ProgressBar, id)

Get a `ProgressBar`'s `ProgressJob` by `id`.
"""
function getjob(pbar::ProgressBar, id)
    idx = findfirst(j -> j.id == id, pbar.jobs)
    isnothing(idx) && return nothing
    return pbar.jobs[idx]
end

"""
    start!(pbar::ProgressBar)

Start a `ProgressBar` and run a `Task` to update its visuals.

Starts a parallel `Task` for updating the `ProgressBar`
visualization while other code runs. The task is stopped
when the progress bar is set to have `running=false`

See also [`stop!`](@ref)
"""
function start!(pbar::ProgressBar)
    pbar.running = true
    print("\n"^(length(pbar.jobs)))
    return nothing
end

"""
    stop!(pbar::ProgressBar)

Stop a running `ProgressBar`.

Stops the `Task` updating the progress bar
visuals too, and if the progress bar was
`transient` it clears up the visuals.
"""
function stop!(pbar::ProgressBar)
    pbar.paused = true
    pbar.running = false

    # if transient, delete
    if pbar.transient
        # move cursor to stale scrollregion and clear
        move_to_line(stdout, console_height())
        for i in 1:(pbar.renderstatus.nlines + 2)
            erase_line(stdout)
            up(stdout)
        end
    else
        print("\n")
    end

    # restore scrollbar region
    change_scroll_region(stdout, console_height())
    show_cursor()
    pbar.transient || print("\n")
    return nothing
end

# --------------------------------- rendering -------------------------------- #
"""
    render(job::ProgressJob, pbar::ProgressBar)::String

Render a `ProgressJob`
"""
function render(job::ProgressJob, pbar::ProgressBar)::String
    color = jobcolor(pbar, job)
    return apply_style(join(update!.(job.columns, color), " "))
end

"""
    render(job::ProgressJob, pbar::ProgressBar)::String

Render a `ProgressJob`
"""
function render(job::ProgressJob)::String
    color = jobcolor(job)
    return apply_style(join(update!.(job.columns, color), " "))
end

"""
    render(job::ProgressJob, pbar::ProgressBar)::String

Render a `ProgressBar`.

When a progress bar is first rendered, this function uses
ANSI codes to change the scrolling region of the terminal
window to create a space at the bottom where the bar's visuals
can be displayed. This allows for thext printed to `stdout` to
still be visualized. On subsequent calls, this function
ensures that the height of the reserved space matches the
number of running jobs.

All fo this requires a bit of careful work in moving the
cursor around and doing ANSI magic.
"""
function render(pbar::ProgressBar)
    # check if running
    pbar.running || return nothing

    # remove completed, transient jobs
    for job in pbar.jobs
        if job.finished && job.transient
            removejob!(pbar, job)
        end
    end

    # get variables
    njobs, height = length(pbar.jobs) + 1, console_height()
    iob = pbar.buff

    # on the first render, create sticky region
    if !pbar.renderstatus.rendered
        print(iob, "\n"^(njobs))
        pbar.renderstatus.scrollline = height - njobs
        change_scroll_region(iob, pbar.renderstatus.scrollline)

        pbar.renderstatus.rendered = true
        pbar.renderstatus.hline =
            string(hLine(pbar.width, "progress"; style = "blue dim")) * "\n"
        pbar.renderstatus.nlines = njobs
        pbar.renderstatus.maxnlines = njobs

    elseif njobs > pbar.renderstatus.maxnlines
        # if we need more lines, scroll
        write(iob, "\n"^(njobs - pbar.renderstatus.maxnlines))

        # set scroll region
        pbar.renderstatus.maxnlines = njobs
        pbar.renderstatus.scrollline = height - pbar.renderstatus.maxnlines
        change_scroll_region(iob, pbar.renderstatus.scrollline)
    end

    # move cursor to scrollregion and clear
    move_to_line(iob, pbar.renderstatus.scrollline + 1)
    cleartoend(iob)

    # render the progressbars
    write(iob, pbar.renderstatus.hline)
    for (last, job) in loop_last(pbar.jobs)
        contents = render(job, pbar)
        coda = last ? "" : "\n"
        write(iob, contents * coda)
    end

    # restore position and write
    move_to_line(iob, pbar.renderstatus.scrollline)
    # write(stdout, take!(iob))

    return print(String(take!(iob)))
end

# ---------------------------------------------------------------------------- #
#                                     WITH                                     #
# ---------------------------------------------------------------------------- #
"""
    with(expr, pbar::ProgressBar)

Wrap an expression to run code in the context of a progress bar.

Ensures that the progress bar is correctly created and distroyed
even if code within the loop whose progress we're monitoring causes
and error. Since `render` changes the scrollregion in the terminal,
we need to make sure that we can restore things no matter what.

# Examples
```julia
pbar = ProgressBar()
with(pbar) do
    job = addjob!(pbar; description="Running")
    for i in 1:500
        update!(job)
        sleep(.0025)
    end
end
```
"""
function with(expr, pbar::ProgressBar)
    # @info expr expr.args

    val = nothing
    try
        start!(pbar)

        task = Threads.@spawn expr()
        while !istaskdone(task) && pbar.running
            pbar.paused || render(pbar)
            sleep(pbar.Δt)
        end
        stop!(pbar)
        var = fetch(task)
    catch err
        stop!(pbar)
        rethrow()
        quit()
    end
    stop!(pbar)
    return val
end

# ---------------------------------------------------------------------------- #
#                                     TRACK                                    #
# ---------------------------------------------------------------------------- #

"""
    @track(ex)

Macro to wrap a loop in a progress bar.

Like [`with`](@ref) it wraps the expression code
in a try/finally statement to ensure that the
progress bar is correctly stopped in any situation.

# Examples
```julia
@track for i in 1:10
    sleep(0.1)
end
```
"""
macro track(ex)
    iter = esc(ex.args[1].args[2])
    i = esc(ex.args[1].args[1])
    body = esc(ex.args[2])
    quote
        __pbar = ProgressBar()
        __pbarjob = addjob!(__pbar; N = length($iter))
        with(__pbar) do
            for $i in $iter
                update!(__pbarjob)
                $body
            end
        end
    end
end

# ------------------------------- general utils ------------------------------ #
"""
jobcolor(job::ProgressJob)

Get the RGB color of of a progress bar's bar based on progress.
"""
function jobcolor(pbar::ProgressBar, job::ProgressJob)
    isnothing(job.N) && return "white"

    α = 0.8 * job.i / job.N
    β = max(sin(π * job.i / job.N) * 0.7, 0.4)

    c1, c2, c3 = pbar.colors
    r = string(int((0.8 - α) * c1.r + β * c2.r + α * c3.r))
    g = string(int((0.8 - α) * c1.g + β * c2.g + α * c3.g))
    b = string(int((0.8 - α) * c1.b + β * c2.b + α * c3.b))
    return "(" * r * ", " * g * ", " * b * ")"
end

const PbarCol1 = RGBColor("(1, .05, .05)")
const PbarCol2 = RGBColor("(.05, .05, 1)")
const PbarCol3 = RGBColor("(.05, 1, .05)")

function jobcolor(job::ProgressJob)
    isnothing(job.N) && return "white"

    α = 0.8 * job.i / job.N
    β = max(sin(π * job.i / job.N) * 0.7, 0.4)

    c1, c2, c3 = PbarCol1, PbarCol2, PbarCol3
    r = string(int((0.8 - α) * c1.r + β * c2.r + α * c3.r))
    g = string(int((0.8 - α) * c1.g + β * c2.g + α * c3.g))
    b = string(int((0.8 - α) * c1.b + β * c2.b + α * c3.b))
    return "(" * r * ", " * g * ", " * b * ")"
end

end


text:
using StrongFieldApproximation
using Documenter

isdefined(Main, :NOPLOTS) && NOPLOTS || include("plots.jl")

makedocs(;
    modules=[StrongFieldApproximation],
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com>",
    repo="https://github.com/jagot/StrongFieldApproximation.jl/blob/{commit}{path}#{line}",
    sitename="StrongFieldApproximation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://www.tipota.org/StrongFieldApproximation.jl",
        assets=String[],
        mathengine = MathJax(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict(
                    :defd => "≝",
                    :abs => ["\\left|#1\\right|", 1],
                    :arcsinh => "\\operatorname{arcsinh}",
                    :ket => ["|#1\\rangle", 1],
                    :bra => ["\\langle#1|", 1],
                    :braket => ["\\langle#1|#2\\rangle", 2],
                    :ketbra => ["|#1\\rangle\\!\\langle#2|", 2],
                    :matrixel => ["\\langle#1|#2|#3\\rangle", 3],
                    :vec => ["\\mathbf{#1}", 1],
                    :mat => ["\\mathsf{#1}", 1],
                    :conj => ["#1^*", 1],
                    :im => "\\mathrm{i}",
                    :operator => ["\\mathfrak{#1}", 1],
                    :Hamiltonian => "\\operator{H}",
                    :hamiltonian => "\\operator{h}",
                    :Lagrangian => "\\operator{L}",
                    :fock => "\\operator{f}",
                    :lagrange => ["\\epsilon_{#1}", 1],
                    :vary => ["\\delta_{#1}", 1],
                    :onebody => ["(#1|#2)", 2],
                    :twobody => ["[#1|#2]", 2],
                    :twobodydx => ["[#1||#2]", 2],
                    :direct => ["{\\operator{J}_{#1}}", 1],
                    :exchange => ["{\\operator{K}_{#1}}", 1],
                    :diff => ["\\mathrm{d}#1\\,", 1],
                    :saddle => ["#1^{(\\textrm{st})}", 1],
                    :ce => "\\mathrm{e}"
                )
            )
        ))
    ),
    pages=[
        "Home" => "index.md",
        "Strong-field ionization" => "strong_field_ionization.md",
        "Reference" => "reference.md",
    ],
    doctest = false
)

deploydocs(;
    repo="github.com/jagot/StrongFieldApproximation.jl",
)


text:
"""
    standalone(b; full)

Internal function to return self-contained HTML with Javascript ready to be
displayed.
"""
function standalone(b::Board; full=false, preamble=true)
    s = ""
    if full
        s = """
            <!DOCTYPE html>

            """
    end
    s *= """

         """
    if full
        s *= """"""
    end
    return s
end

# NOTE: Juno/Atom does not allow interactive Javascript so obliged to go via
# a Blink window. Not ideal but ok.
#
# For IJulia, should manage to do it like
# https://github.com/queryverse/VegaLite.jl/blob/2208264fe0bfd38f563f26035dd00a0153bd0c61/src/rendering/render.jl#L74
# however this is black magic, I haven't  worked out how to  do it.

function Base.show(io::IO, b::Board)
    if isempty(b.objects)
        if isempty(b.functions)
            println(io, "Board $(b.name) (empty).")
        else
            println(io, "Board $(b.name) (no objects).")
        end
        return
    elseif isdefined(Main, :Atom) && Main.Atom.PlotPaneEnabled[]
        p = Blink.Page()
        Main.Atom.ploturl(Blink.localurl(p))
        wait(p)
        Blink.body!(p, standalone(b, full=true))
    else
        w = Blink.Window()
        Blink.body!(w, standalone(b, full=true))
    end
    return nothing
end


text:
# This file implements one solution method for the network decomposition
# using LagrangeDual from DualDecomposition.jl

function init_DD_algo(dn_model::NetDecModel)::DD.LagrangeDual
    algo = DD.LagrangeDual()

    partition = get_partition(dn_model)
    models = get_models(dn_model)
    for i in eachindex(partition)
        DD.add_block_model!(algo, i, models[i].model)
    end
    add_split_vars_to_algo!(dn_model, algo)
    return algo
end

function add_split_vars_to_algo!(dn_model::NetDecModel, algo::DD.AbstractMethod)
    coupling_vars = Vector{DD.CouplingVariableRef}()
    partition = get_partition(dn_model)
    split_vars = get_split_vars(dn_model)
    for i in eachindex(partition)
        vars_dict = split_vars[i]
        for (varname, dict) in vars_dict
            for (idx, vref) in dict
                # id = varname * "_" * string(idx)
                id = (varname, idx)
                push!(coupling_vars, DD.CouplingVariableRef(i, id, vref))
            end
        end
    end
    DD.set_coupling_variables!(algo, coupling_vars)
    return
end