text:
# postprocess_samples.jl
# 2020-01-17
# David Merrell
#
# Functions for postprocessing MCMC sample output.
# Specifically for computing convergence statistics
# (PSRF and N_eff).

module SamplePostprocess

using JSON
using ArgParse
using DataStructures

include("changepoint_vec.jl")
include("nested_data.jl")

export postprocess_sample_files, julia_main

"""
Read an MCMC samples output file
"""
function load_samples(samples_file_path::String)

    f = open(samples_file_path, "r")
    str = read(f, String)
    close(f)
    d = JSON.parse(str)

    d["parent_sets"] = [DefaultDict([], ps) for ps in d["parent_sets"]]
    V = length(d["parent_sets"])
    d["parent_sets"] = [[ps[string(p)] for p=1:V] for ps in d["parent_sets"]]

    return d
end

function seq_variogram(seq::ChangepointVec, t::Int64; is_binary::Bool=false)
    n = length(seq)
    if is_binary
        return mean(binop(!=, seq[1:n-t], seq[1+t:n]); is_binary=is_binary)
    else
        return mean(map(abs2, binop(-, seq[1:n-t], seq[1+t:n])); is_binary=is_binary)
    end
end

function seq_variogram(seq::Vector, t::Int64; is_binary::Bool=false)
    n = length(seq)
    if is_binary
        return sum(seq[1:n-t] .!= seq[1+t:n]) / (n - t)
    else
        return sum(abs2.(seq[1:n-t] .- seq[1+t:n])) / (n - t)
    end
end

function correlation_sum(seqs::Vector{ChangepointVec}, varplus::Float64;
                         is_binary::Bool=false, density_threshold::Float64=0.001)

    s = 0.0
    m = length(seqs)
    prev2_corr = 0.0
    prev1_corr = 0.0

    # Determine whether the sequences need to be converted to dense vectors.
    seq_densities = [length(seq.changepoints)/length(seq) for seq in seqs]
    seqs = [(seq_densities[i] < density_threshold ? seq : to_vec(seq) ) for (i, seq) in enumerate(seqs)]

    for t=1:length(seqs[1])-2

        variograms = [seq_variogram(seq, t; is_binary=is_binary) for seq in seqs]
        comb_variogram = sum(variograms) / m

        corr = 1.0 - (comb_variogram / (2.0*varplus))
        if t > 2
            if corr + prev1_corr >= 0.0
                s += prev2_corr
            else
                break
            end
        end
        prev2_corr = prev1_corr
        prev1_corr = corr

    end

    return s
end

n_eff_stat(corr_sum::Float64, m::Int64, n::Int64) = m*n/(1.0 + 2.0*corr_sum)

#################################################
# These functions are used for computing PSRF and
# N_eff (diagnostic scores for MCMC convergence)
#################################################

function seq_var(cpv::ChangepointVec, mean::Float64; is_binary::Bool=false)
    if is_binary
        N = length(cpv)
        return mean*(1 - mean)*N/(N - 1)
    else
        return sum(map( x->abs2(x-mean), cpv )) / (length(cpv) - 1)
    end
end

w_stat(variances::Vector{Float64}) = sum(variances) / length(variances)

b_stat(means::Vector{Float64}, n::Int64) = n*sum(abs2.(means .- sum(means)/length(means)))/length(means)

varplus(w_stat::Float64, b_stat::Float64, n::Int64) = (n-1)*w_stat/n + b_stat/n

psrf_stat(varplus, w_stat) = sqrt(varplus / w_stat)

"""
Compute the PSRF and N_eff for a collection of sequences.
In the context of MCMC, we assume the seqs have already been
split in half.

We assume the seqs all have equal length.
"""
function compute_psrf_neff(seqs::Vector{ChangepointVec}; is_binary::Bool=false)
    seq_means = [mean(seq; is_binary=is_binary) for seq in seqs]
    seq_variances = [seq_var(seq, seq_means[i]; is_binary=is_binary) for (i, seq) in enumerate(seqs)]
    m = length(seqs)
    n = length(seqs[1])
    B = b_stat(seq_means, n)
    W = w_stat(seq_variances)
    vp = varplus(W, B, n)
    psrf = psrf_stat(vp, W)

    corr_sum = correlation_sum(seqs, vp; is_binary=is_binary)
    n_eff = n_eff_stat(corr_sum, m, n)

    return psrf, n_eff
end

function is_binary(cpv)
    for pair in cpv.changepoints
        if !(pair[2] == 0 || pair[2] == 1)
            return false
        end
    end
    return true
end

"""
Compute convergence diagnostics for a collection of sequences,
at multiple stop indices
"""
function evaluate_convergence(whole_seqs::Vector{ChangepointVec}, stop_idxs;
                              burnin::Float64=0.5)

    # Determine whether the data is binary
    isbin = is_binary(whole_seqs[1])

    # We'll collect a list of (psrf, n_eff values)
    diagnostics = []

    for len in stop_idxs
        # Split the whole sequences into half sequences
        burnin_idx = Int(floor(burnin*len))
        split_idx = burnin_idx + div(len - burnin_idx, 2)
        half_seqs = Vector{ChangepointVec}()
        for ws in whole_seqs
            push!(half_seqs, ws[burnin_idx+1:split_idx])
            push!(half_seqs, ws[split_idx+1:len])
        end

        # compute the convergence diagnostics
        psrf, n_eff = compute_psrf_neff(half_seqs; is_binary=isbin)
        push!(diagnostics, (psrf, n_eff))
    end

    return diagnostics
end

function evaluate_summary(whole_seqs::Vector{ChangepointVec}, stop_idxs;
                          burnin::Float64=0.5)

    burnin_idxs = [Int(floor(burnin*length(seq))) for seq in whole_seqs]
    sums = [sum(cpv[(burnin_idxs[i]+1):length(cpv)]) for (i, cpv) in enumerate(whole_seqs)]
    ns = [length(cpv) - burnin_idxs[i] for (i, cpv) in enumerate(whole_seqs)]

    return sum(sums) / sum(ns)
end

function is_seq(node)
    if !(typeof(node) <: AbstractVector)
        return false
    end
    for (i, item) in enumerate(node)
        if !(length(item) == 2) || !(typeof(item[1]) <: Number)
            return false
        end
    end
    return true
end

is_cpv = x -> (typeof(x) <: ChangepointVec)

###############################################################
# CORE FUNCTIONS
###############################################################
function collect_stats(dict_vec, stop_points; burnin::Float64=0.5)

    results = Dict()
    results["stop_points"] = stop_points

    # Convert all of the sequences to ChangepointVecs
    dict_vec = [leaf_map(x -> ChangepointVec(x, d["n"]), d, is_seq) for d in dict_vec]

    results["summary"] = aggregate_trees(x -> evaluate_summary(x, stop_points; burnin=burnin),
                                         dict_vec, is_cpv)
    results["conv_stats"] = aggregate_trees(x -> evaluate_convergence(x, stop_points; burnin=burnin),
                                            dict_vec, is_cpv)

    results["parent_sets"] = results["summary"]["parent_sets"]
    results["edge_conf_key"] = "parent_sets"

    results["t_elapsed"] = [d["t_elapsed"] for d in dict_vec]
    results["n"] = [d["n"] for d in dict_vec]

    delete!(results["summary"], "parent_sets")
    return results
end

function get_max_n(dict_vec)
    n_vec = [d["n"] for d in dict_vec]
    return minimum(n_vec)
end

function postprocess_sample_files(sample_filenames, output_file::String, stop_points::Vector{Int},
                                  burnin::Float64)
    dict_vec = [load_samples(fname) for fname in sample_filenames]

    n_max = get_max_n(dict_vec)
    if length(stop_points) == 0
        stop_points = [n_max]
    else
        stop_points = [sp for sp in stop_points if sp <= n_max]
    end

    results = collect_stats(dict_vec, stop_points; burnin=burnin)

    f = open(output_file, "w")
    JSON.print(f, results)
    close(f)
    println("Saved postprocessed MCMC results to ", output_file)
end

function get_args(args::Vector{String})

    s = ArgParseSettings()
    @add_arg_table! s begin
        "--chain-samples"
            help = "mcmc sample files for one or more chains"
    	    required = true
            nargs = '+'
        "--stop-points"
            help = "points in the sequence where convergence should be analyzed."
            required = false
            nargs = '+'
            default = Vector{Int64}()
            arg_type = Int
        "--burnin"
            help = "fraction of sequence to discard as burnin."
            required = false
            arg_type = Float64
            default = 0.5
        "--output-file"
            help = "name of output JSON file containing convergence information"
    	    required = true
            arg_type = String
    end
    args = parse_args(args, s)

    arg_vec = [args["chain-samples"], args["output-file"], args["stop-points"],
               args["burnin"]]
    return arg_vec
end

# main function -- for purposes of static compilation
function julia_main()
    arg_vec = get_args(ARGS)
    postprocess_sample_files(arg_vec...)
    return 0
end

end

using .SamplePostprocess

julia_main()

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
function hidden_signal_loglikelihood_CIR(trajectory, times, δ, γ, σ)
    ll::Float64 = 0
    for k in 1:(length(trajectory)-1)
        ll += FeynmanKacParticleFilters.CIR_transition_logdensity(trajectory[k+1], trajectory[k], times[k+1]-times[k], δ, γ, σ)
    end
    return ll
end

function emitted_data_conditional_loglikelihood_CIR(data, trajectory, times, λ)
    ll::Float64 = 0
    for k in eachindex(times)
        t = times[k]
        ll += sum(logpdf.(Poisson(λ*trajectory[k]), data[t]))
    end
    return ll
end

function joint_loglikelihood_CIR(data, trajectory, times, δ, γ, σ, λ)
    return hidden_signal_loglikelihood_CIR(trajectory, times, δ, γ, σ) + emitted_data_conditional_loglikelihood_CIR(data, trajectory, times, λ)
end

function joint_sampler_CIR_pruning_precompute(data, λ, prior_logpdf, θ_init, niter, do_the_pruning::Function; final_chain_length = 1000, silence = false, jump_size = 0.5*Matrix(I,3,3), reparam = false)

    trajectory_chain = Array{Float64,2}(undef, niter, length(data |> keys))
    parameter_chain = Array{Float64,2}(undef, niter, 3)
    times = data |> keys |> collect |> sort

    #Not using the logfunction because not implemented yet
    function unnormalised_logposterior_reparam(a, b, σ_prime, trajectory)::Float64
        prior_contribution = prior_logpdf(a, b, σ_prime)
        if isinf(prior_contribution) #Avoids likelihood computation in this case
           return prior_contribution
        else
            return prior_contribution + reparam_joint_loglikelihood_CIR(data, trajectory, times, a, b, σ_prime, λ)::Float64
        end
    end
    function unnormalised_logposterior(δ, γ, σ, trajectory)::Float64
        prior_contribution = prior_logpdf(δ, γ, σ)
        if isinf(prior_contribution) #Avoids likelihood computation in this case
            return prior_contribution
        else
            return prior_contribution + joint_loglikelihood_CIR(data, trajectory, times, δ, γ, σ, λ)::Float64
        end
    end

    θ_it = θ_init

    Jtsym_rand(θ) = rand(MvNormal(zeros(3), jump_size)) .+ θ

    for it in 1:niter

        print_every = max(1, floor(niter/20))

        if !silence && mod(it, print_every) == 0
            @info "$it iterations out of $niter"
        end

        # Sample trajectory

        #Not using the logfunction because not implemented yet
        if reparam
            θ_it_δγ_param = inverse_reparam_CIR(θ_it...)
            Λ_of_t, wms_of_t, θ_of_t = filter_CIR_pruning(θ_it_δγ_param[1], θ_it_δγ_param[2], θ_it_δγ_param[3], λ, data, do_the_pruning, silence = true)
            Λ_of_t_pruned, wms_of_t_pruned =  prune_all_dicts(Λ_of_t, wms_of_t, do_the_pruning)
            # logwms_of_t = log.(wms_of_t)
            # logwms_pred_of_t = log.(wms_pred_of_t)
            # 1, 1, 1 are dummy arguments
            X_it = sample_1_trajectory_from_joint_smoothing_CIR_precompute(θ_it_δγ_param[1], θ_it_δγ_param[2], θ_it_δγ_param[3], Λ_of_t_pruned, wms_of_t_pruned, θ_of_t, 1, 1, 1, data)
        else
            Λ_of_t, wms_of_t, θ_of_t = filter_CIR_pruning(θ_it[1], θ_it[2], θ_it[3], λ, data, do_the_pruning, silence = true)
            Λ_of_t_pruned, wms_of_t_pruned =  prune_all_dicts(Λ_of_t, wms_of_t, do_the_pruning)
            # logwms_of_t = log.(wms_of_t)
            # logwms_pred_of_t = log.(wms_pred_of_t)
            # 1, 1, 1 are dummy arguments
            X_it = sample_1_trajectory_from_joint_smoothing_CIR_precompute(θ_it[1], θ_it[2], θ_it[3], Λ_of_t_pruned, wms_of_t_pruned, θ_of_t, 1, 1, 1, data)
        end

        @debug "iteration $it"
        @debug "trajectory" X_it

        # Metropolis-Hastings step

        if reparam
            θ_it = draw_next_sample(θ_it, Jtsym_rand, st -> unnormalised_logposterior_reparam(st[1], st[2], st[3], X_it))
        else
            θ_it = draw_next_sample(θ_it, Jtsym_rand, st -> unnormalised_logposterior(st[1], st[2], st[3], X_it))
        end
        # Record iteration

        trajectory_chain[it,:] = X_it
        parameter_chain[it,:] = collect(θ_it)

    end

    return parameter_chain, trajectory_chain
end

function joint_sampler_CIR_keep_fixed_number_precompute(data, λ, prior_logpdf, θ_init, niter, fixed_number::Int64; final_chain_length = 1000, silence = false, jump_size = 0.5*Matrix(I,3,3))
    # println("filter_WF_mem2")

    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    joint_sampler_CIR_pruning_precompute(data, λ, prior_logpdf, θ_init, niter, prune_keeping_fixed_number; final_chain_length = 1000, silence = silence, jump_size = jump_size)

end

ext: jl

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# KINEMATICS
@doc doc"""
    rigid_body_velocity(vel_P, ω, r_PQ)

Calculate the velocity of a point Q of a rigid solid given the velocity of a
point P (vel_P), the rotational velocity of the solid (ω) and the relative
position of Q wrt P.

If the reference frame 1 is attached to the solid and the velocity is
calculated with respect to reference frame 0:

``v_{10}^{Q} = v_{10}^{P} + \omega_{10} \times r^{PQ}``

being:
- ``v_{10}^{Q}`` the velocity of point Q, fixed to 1, wrt 0
- ``\omega_{10}`` the angular velocity of the solid 1 wrt 0
- ``r^{PQ}`` the position of Q wrt P (``r^{Q}-r^{P}``)

Every vector needs to be expressed in the same coordinate system.

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (Section 1.3, page 26)
"""
function rigid_body_velocity(vel_P, ω, r_PQ)
    vel_Q = vel_P + cross(ω, r_PQ)
    return vel_Q
end

@doc doc"""
    rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)

Calculate the acceleration of a point Q of a rigid solid given the acceleration
of a point P (acc_P), the rotational velocity of the solid (ω), the rotational
acceleration of the solid (ω_dot) and the relative position of Q wrt P.

``a_{10}^{Q} = a_{10}^{P} + \omega_{10} \times (\omega_{10} \times r^{PQ}) + \dot{\omega}_{10} \times r^{PQ}``

being:
- ``a_{10}^{Q}`` the acceleration of point Q, fixed to 1, wrt 0
- ``\omega_{10}`` the angular velocity of the solid 1 wrt 0
- ``\dot{\omega}_{10}`` the angular acceleration of the solid 1 wrt 0
- ``r^{PQ}`` the position of Q wrt P (``r^{Q}-r^{P}``)

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (Section 1.3, Formaula 1.3-14c, page 26)
"""
function rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)
    acc_Q = acc_P + cross(ω, cross(ω, r_PQ)) + cross(ω_dot, r_PQ)
    return acc_Q
end

# Mechanics
function steiner_inertia(point0, inertia0, mass, point1)
    r = point1 - point0
    inertia0 + mass * (dot(r, r) * I - r*r')
end


text:
# The number of Gaussian components in the PSF.
const psf_K = 2

"""
A single normal component of the point spread function.
All quantities are in pixel coordinates.

Args:
  alphaBar: The scalar weight of the component.
  xiBar: The 2x1 location vector
  tauBar: The 2x2 covariance

Attributes:
  alphaBar: The scalar weight of the component.
  xiBar: The 2x1 location vector
  tauBar: The 2x2 covariance (tau_bar in the ICML paper)
  tauBarInv: The 2x2 precision
  tauBarLd: The log determinant of the covariance
"""
immutable PsfComponent
    alphaBar::Float64  # TODO: use underscore
    xiBar::Vector{Float64}
    tauBar::Matrix{Float64}

    tauBarInv::Matrix{Float64}
    tauBarLd::Float64

    function PsfComponent(alphaBar::Float64, xiBar::Vector{Float64},
                          tauBar::Matrix{Float64})
        new(alphaBar, xiBar, tauBar, tauBar^-1, logdet(tauBar))
    end
end

"""
SDSS representation of a spatially variable PSF. The PSF is represented as
a weighted combination of eigenimages (stored in `rrows`), where the weights
vary smoothly across the image as a polynomial of the form

```
weight[k](x, y) = sum_{i,j} cmat[i, j, k] * (rcs * x)^i (rcs * y)^j
```

where `rcs` is a coordinate transformation and `x` and `y` are zero-indexed.
"""
immutable RawPSF
    rrows::Array{Float64,2}  # A matrix of flattened eigenimages.
    rnrow::Int  # The number of rows in an eigenimage.
    rncol::Int  # The number of columns in an eigenimage.
    cmat::Array{Float64,3}  # The coefficients of the weight polynomial

    function RawPSF(rrows::Array{Float64, 2}, rnrow::Integer, rncol::Integer,
                     cmat::Array{Float64, 3})
        # rrows contains eigen images. Each eigen image is along the first
        # dimension in a flattened form. Check that dimensions match up.
        @assert size(rrows, 1) == rnrow * rncol

        # The second dimension is the number of eigen images, which should
        # match the number of coefficient arrays.
        @assert size(rrows, 2) == size(cmat, 3)

        return new(rrows, Int(rnrow), Int(rncol), cmat)
    end
end

function get_psf_width(psf::Array{PsfComponent}; width_scale=1.0)
    # A heuristic measure of the PSF width based on an anology
    # with it being a mixture of normals.    Note that it is not an actual
    # mixture of normals, and in particular that sum(alphaBar) \ne 1.

    # The PSF is not necessarily centered at (0, 0), but we want a measure
    # of its maximal width around (0, 0), not around its center.
    # Approximate this by finding the covariance of a point randomly drawn
    # from a mixture of gaussians.
    alpha_norm = sum([ psf_comp.alphaBar for psf_comp in psf ])
    cov_est = zeros(Float64, 2, 2)
    for psf_comp in psf
        cov_est +=
            psf_comp.alphaBar * (psf_comp.xiBar * psf_comp.xiBar' + psf_comp.tauBar) /
            alpha_norm
    end

    # Return the twice the sd of the most spread direction, scaled by the total
    # mass in the PSF.
    width_scale * sqrt(eigvals(cov_est)[end]) * alpha_norm
end


text:
export
# Interpolations
    Interpolation,
    RefCube,
    RefTetrahedron,
    Lagrange,
    Serendipity,
    getnbasefunctions,

# Quadrature
    QuadratureRule,
    getweights,
    getpoints,

# FEValues
    CellValues,
    ScalarValues,
    VectorValues,
    CellScalarValues,
    CellVectorValues,
    FaceValues,
    FaceScalarValues,
    FaceVectorValues,
    reinit!,
    shape_value,
    shape_gradient,
    shape_symmetric_gradient,
    shape_divergence,
    shape_curl,
    function_value,
    function_gradient,
    function_symmetric_gradient,
    function_divergence,
    function_curl,
    spatial_coordinate,
    getnormal,
    getdetJdV,
    getnquadpoints,

# Grid
    Grid,
    Node,
    Cell,
    Line,
    Line2D,
    Line3D,
    QuadraticLine,
    Triangle,
    QuadraticTriangle,
    Quadrilateral,
    Quadrilateral3D,
    QuadraticQuadrilateral,
    Tetrahedron,
    QuadraticTetrahedron,
    Hexahedron,
    QuadraticHexahedron,
    CellIndex,
    FaceIndex,
    EdgeIndex,
    VertexIndex,
    getcells,
    getncells,
    getnodes,
    getnnodes,
    getcelltype,
    getcellset,
    getnodeset,
    getfaceset,
    getedgeset,
    getvertexset,
    getcoordinates,
    getcoordinates!,
    getcellsets,
    getnodesets,
    getfacesets,
    getedgesets,
    getvertexsets,
    onboundary,
    nfaces,
    addnodeset!,
    addfaceset!,
    addedgeset!,
    addvertexset!,
    addcellset!,
    transform!,
    generate_grid,
    compute_vertex_values,

# Dofs
    DofHandler,
    close!,
    ndofs,
    ndofs_per_cell,
    celldofs!,
    celldofs,
    create_sparsity_pattern,
    create_symmetric_sparsity_pattern,
    dof_range,
    renumber!,
    MixedDofHandler,
    FieldHandler,
    Field,

# Constraints
    ConstraintHandler,
    Dirichlet,
    update!,
    apply!,
    apply_rhs!,
    get_rhs_data,
    apply_zero!,
    add!,
    free_dofs,

# iterators
    CellIterator,
    UpdateFlags,
    cellid,

# assembly
    start_assemble,
    assemble!,
    end_assemble,

# VTK export
    vtk_grid,
    vtk_point_data,
    vtk_cell_data,
    vtk_nodeset,
    vtk_cellset,
    vtk_save,

# L2 Projection
    project,
    L2Projector


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

# ––––––––
# Emphasis
# ––––––––

mutable struct Italic
    text
end

@trigger '*' ->
function asterisk_italic(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "*")
    return result === nothing ? nothing : Italic(parseinline(result, md))
end

@trigger '_' ->
function underscore_italic(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "_")
    return result === nothing ? nothing : Italic(parseinline(result, md))
end

mutable struct Bold
    text
end

@trigger '*' ->
function asterisk_bold(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "**")
    return result === nothing ? nothing : Bold(parseinline(result, md))
end

@trigger '_' ->
function underscore_bold(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "__")
    return result === nothing ? nothing : Bold(parseinline(result, md))
end

# ––––
# Code
# ––––

@trigger '`' ->
function inline_code(stream::IO, md::MD)
    withstream(stream) do
        ticks = startswith(stream, r"^(`+)")
        result = readuntil(stream, ticks)
        if result === nothing
            nothing
        else
            result = strip(result)
            # An odd number of backticks wrapping the text will produce a `Code` node, while
            # an even number will result in a `LaTeX` node. This allows for arbitary
            # backtick combinations to be embedded inside the resulting node, i.e.
            #
            # `a`, ``a``, `` `a` ``, ``` ``a`` ```, ``` `a` ```, etc.
            #  ^     ^        ^            ^             ^
            #  C     L        L            C             C       with C=Code and L=LaTeX.
            #
            isodd(length(ticks)) ? Code(result) : LaTeX(result)
        end
    end
end

# ––––––––––––––
# Images & Links
# ––––––––––––––

mutable struct Image
    url::String
    alt::String
end

@trigger '!' ->
function image(stream::IO, md::MD)
    withstream(stream) do
        startswith(stream, "![") || return
        alt = readuntil(stream, ']', match = '[')
        alt ≡ nothing && return
        skipwhitespace(stream)
        startswith(stream, '(') || return
        url = readuntil(stream, ')', match = '(')
        url ≡ nothing && return
        return Image(url, alt)
    end
end

mutable struct Link
    text
    url::String
end

@trigger '[' ->
function link(stream::IO, md::MD)
    withstream(stream) do
        startswith(stream, '[') || return
        text = readuntil(stream, ']', match = '[')
        text ≡ nothing && return
        skipwhitespace(stream)
        startswith(stream, '(') || return
        url = readuntil(stream, ')', match = '(')
        url ≡ nothing && return
        return Link(parseinline(text, md), url)
    end
end

@trigger '[' ->
function footnote_link(stream::IO, md::MD)
    withstream(stream) do
        regex = r"^\[\^(\w+)\]"
        str = startswith(stream, regex)
        if isempty(str)
            return
        else
            ref = match(regex, str).captures[1]
            return Footnote(ref, nothing)
        end
    end
end

@trigger '<' ->
function autolink(stream::IO, md::MD)
    withstream(stream) do
        startswith(stream, '<') || return
        url = readuntil(stream, '>')
        url ≡ nothing && return
        _is_link(url) && return Link(url, url)
        _is_mailto(url) && return Link(url, url)
        return
    end
end

# This list is taken from the commonmark spec
# http://spec.commonmark.org/0.19/#absolute-uri
const _allowable_schemes = Set(split("coap doi javascript aaa aaas about acap cap cid crid data dav dict dns file ftp geo go gopher h323 http https iax icap im imap info ipp iris iris.beep iris.xpc iris.xpcs iris.lwz ldap mailto mid msrp msrps mtqp mupdate news nfs ni nih nntp opaquelocktoken pop pres rtsp service session shttp sieve sip sips sms snmp,soap.beep soap.beeps tag tel telnet tftp thismessage tn3270 tip tv urn vemmi ws wss xcon xcon-userid xmlrpc.beep xmlrpc.beeps xmpp z39.50r z39.50s
adiumxtra afp afs aim apt,attachment aw beshare bitcoin bolo callto chrome,chrome-extension com-eventbrite-attendee content cvs,dlna-playsingle dlna-playcontainer dtn dvb ed2k facetime feed finger fish gg git gizmoproject gtalk hcp icon ipn irc irc6 ircs itms jar jms keyparc lastfm ldaps magnet maps market,message mms ms-help msnim mumble mvn notes oid palm paparazzi platform proxy psyc query res resource rmi rsync rtmp secondlife sftp sgn skype smb soldat spotify ssh steam svn teamspeak
things udp unreal ut2004 ventrilo view-source webcal wtai wyciwyg xfire xri ymsgr"))

function _is_link(s::AbstractString)
    '<' in s && return false

    m = match(r"^(.*)://(\S+?)(:\S*)?$", s)
    m ≡ nothing && return false
    scheme = lowercase(m.captures[1])
    return scheme in _allowable_schemes
end

# non-normative regex from the HTML5 spec
const _email_regex = r"^[a-zA-Z0-9.!#$%&'*+/=?^_`{|}~-]+@[a-zA-Z0-9](?:[a-zA-Z0-9-]{0,61}[a-zA-Z0-9])?(?:\.[a-zA-Z0-9](?:[a-zA-Z0-9-]{0,61}[a-zA-Z0-9])?)*$"

function _is_mailto(s::AbstractString)
    length(s) < 6 && return false
    # slicing strings is a bit risky, but this equality check is safe
    lowercase(s[1:6]) == "mailto:" || return false
    return occursin(_email_regex, s[6:end])
end

# –––––––––––
# Punctuation
# –––––––––––

mutable struct LineBreak end

@trigger '\\' ->
function linebreak(stream::IO, md::MD)
    if startswith(stream, "\\\n")
        return LineBreak()
    end
end

@trigger '-' ->
function en_dash(stream::IO, md::MD)
    if startswith(stream, "--")
        return "–"
    end
end

const escape_chars = "\\`*_#+-.!{}[]()\$"

@trigger '\\' ->
function escapes(stream::IO, md::MD)
    withstream(stream) do
        if startswith(stream, "\\") && !eof(stream) && (c = read(stream, Char)) in escape_chars
            return string(c)
        end
    end
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
log_temperature(t::Real) = 1 / log(t)

constant_temperature(t::Real) = 1.0

function default_neighbor!(x::Array, x_proposal::Array)
    @assert size(x) == size(x_proposal)
    for i in 1:length(x)
        @inbounds x_proposal[i] = x[i] + randn()
    end
    return
end

immutable SimulatedAnnealing <: Optimizer
    neighbor!::Function
    temperature::Function
    keep_best::Bool # not used!?
end

SimulatedAnnealing(; neighbor!::Function = default_neighbor!,
                     temperature::Function = log_temperature,
                     keep_best::Bool = true) =
  SimulatedAnnealing(neighbor!, temperature, keep_best)

type SimulatedAnnealingState{T}
    @add_generic_fields()
    iteration::Int64
    x_current::Array{T}
    x_proposal
    f_x_current::T
    f_proposal::T
end

initial_state(method::SimulatedAnnealing, options, d, initial_x::Array) = initial_state(method, options, d.f, initial_x)

function initial_state{T}(method::SimulatedAnnealing, options, f::Function, initial_x::Array{T})
    # Count number of parameters
    n = length(initial_x)

    # Store f(x) in f_x
    f_x = f(initial_x)
    f_calls = 1

    # Store the best state ever visited
    best_x = copy(initial_x)
    best_f_x = f_x
    SimulatedAnnealingState("Simulated Annealing", n, copy(initial_x), f_x, f_calls, 0, 0, 1, copy(initial_x), best_x, best_f_x, f_x)
end

update_state!(d, state::SimulatedAnnealingState, method::SimulatedAnnealing) = update_state!(d.f, state, method)
function update_state!{T}(f::Function, state::SimulatedAnnealingState{T}, method::SimulatedAnnealing)

    # Determine the temperature for current iteration
    t = method.temperature(state.iteration)

    # Randomly generate a neighbor of our current state
    method.neighbor!(state.x_current, state.x_proposal)

    # Evaluate the cost function at the proposed state
    state.f_proposal = f(state.x_proposal)
    state.f_calls += 1

    if state.f_proposal <= state.f_x_current
        # If proposal is superior, we always move to it
        copy!(state.x_current, state.x_proposal)
        state.f_x_current = state.f_proposal

        # If the new state is the best state yet, keep a record of it
        if state.f_proposal < state.f_x
            state.f_x = state.f_proposal
            copy!(state.x, state.x_proposal)
        end
    else
        # If proposal is inferior, we move to it with probability p
        p = exp(-(state.f_proposal - state.f_x_current) / t)
        if rand() <= p
            copy!(state.x_current, state.x_proposal)
            state.f_x_current = state.f_proposal
        end
    end

    state.iteration += 1
    false
end


text:
struct ExponentialDec
    Di
    Qi
    Ti
end


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

@testset ExtendedTestSet "`Many`" begin
    enc = ImagePreprocessing()
    block = Many(Image{2}())
    FastAI.testencoding(enc, block)
end


text:
using TriMatrices
using TriMatrices: trinum, hasdiag, nelems, transpose_layout

LAYOUT_TYPES = [TriLower, TriUpper, TriSymmetric]

@testset "Type parameter" begin
	for L in LAYOUT_TYPES
		# Shouldn't throw any errors
		L{true}()
		L{false}()

		# Default to true
		@test typeof(L()) === L{true}

		# Non-bool parameter
		@test_throws ArgumentError L{Int}()
	end
end

@testset "Methods and type traits" begin
	n = 10

	for L in LAYOUT_TYPES
		for D in [true, false]
			layout = L{D}()

			@test hasdiag(layout) === D
			@test hasdiag(typeof(layout)) === D

			ne = trinum(D ? n : n - 1)
			@test nelems(layout, n) == ne
			@test nelems(typeof(layout), n) == ne
		end

		# Type parameter needed
		@test_throws MethodError hasdiag(L)
		@test_throws MethodError nelems(L, n)
	end
end

@testset "transpose_layout" begin
	@test transpose_layout(TriLower) === TriUpper
	@test transpose_layout(TriUpper) === TriLower
	@test transpose_layout(TriSymmetric) === TriSymmetric
	@test transpose_layout(TriLower{true}) === TriUpper{true}
	@test transpose_layout(TriUpper{true}) === TriLower{true}
	@test transpose_layout(TriSymmetric{true}) === TriSymmetric{true}
	@test transpose_layout(TriLower{false}) === TriUpper{false}
	@test transpose_layout(TriUpper{false}) === TriLower{false}
	@test transpose_layout(TriSymmetric{false}) === TriSymmetric{false}

	# Check transpose of instance matches transpose of its type
	for L in LAYOUT_TYPES
		for D in [true, false]
			layout = L{D}()
			@test typeof(transpose_layout(layout)) === transpose_layout(typeof(layout))
		end
	end
end


text:
__precompile__(false)

module CMPFit

using Printf
using Pkg

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("CMPFit not properly installed. Please run Pkg.build(\"CMPFit\")")
end

# Exported symbols
export cmpfit

######################################################################
# Private definitions
######################################################################

#---------------------------------------------------------------------
"Parameter info strutcture"
mutable struct Parinfo
    fixed::Cint                # 1 = fixed; 0 = free
    limited::NTuple{2, Cint}   # 1 = low/upper limit; 0 = no limit
    limits::NTuple{2, Cdouble} # lower/upper limit boundary value
    char::Ptr{Cchar}           # Name of parameter, or 0 for none
    step::Cdouble              # Step size for finite difference
    relstep::Cdouble           # Relative step size for finite difference
    side::Cint                 #= Sidedness of finite difference derivative
                                  0 - one-sided derivative computed automatically
                                  1 - one-sided derivative (f(x+h) - f(x)  )/h
                                 -1 - one-sided derivative (f(x)   - f(x-h))/h
                                  2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)
                                  3 - user-computed analytical derivatives
                               =#
    deriv_debug::Cint          #= Derivative debug mode: 1 = Yes; 0 = No;
                                  If yes, compute both analytical and numerical
                                  derivatives and print them to the console for
                                  comparison.

                                  NOTE: when debugging, do *not* set side = 3,
                                  but rather to the kind of numerical derivative
                                  you want to compare the user-analytical one to
                                  (0, 1, -1, or 2).
                               =#

    deriv_reltol::Cdouble      # Relative tolerance for derivative debug printout
    deriv_abstol::Cdouble      # Absolute tolerance for derivative debug printout

    "Create an empty `Parinfo` structure."
    Parinfo() = new(0, (0, 0), (0, 0), 0, 0, 0, 0, 0, 0, 0)
end

#---------------------------------------------------------------------
"Create a `Vector{Parinfo}` of empty `Parinfo` structures, whose length is given by `npar`."
function Parinfo(npar::Int)
    return [Parinfo() for i in 1:npar]
end

#---------------------------------------------------------------------
#Define sibling structure with toggled mutability
struct imm_Parinfo
    fixed::Cint
    limited::NTuple{2, Cint}
    limits::NTuple{2, Cdouble}
    char::Ptr{Cchar}
    step::Cdouble
    relstep::Cdouble
    side::Cint
    deriv_debug::Cint
    deriv_reltol::Cdouble
    deriv_abstol::Cdouble
end

function imm_Parinfo(p::Parinfo)
    imm_Parinfo(ntuple((i->begin
                        getfield(p, i)
                        end), nfields(p))...)
end

#---------------------------------------------------------------------
"CMPFit config structure"
mutable struct Config
    ftol::Cdouble        # Relative chi-square convergence criterium Default: 1e-10
    xtol::Cdouble        # Relative parameter convergence criterium  Default: 1e-10
    gtol::Cdouble        # Orthogonality convergence criterium       Default: 1e-10
    epsfcn::Cdouble      # Finite derivative step size               Default: eps()
    stepfactor::Cdouble  # Initial step bound                        Default: 100.0
    covtol::Cdouble      # Range tolerance for covariance calculation Default: 1e-14
    maxiter::Cint        # Maximum number of iterations.  If maxiter == MP_NO_ITER,
                         # then basic error checking is done, and parameter
                         # errors/covariances are estimated based on input
                         # parameter values, but no fitting iterations are done.
                         # Default: 200

    maxfev::Cint;        # Maximum number of function evaluations, or 0 for no limit
                         # Default: 0 (no limit)
    nprint::Cint;        # Default: 1
    douserscale::Cint    # Scale variables by user values?
                         # 1 = yes, user scale values in diag;
                         # 0 = no, variables scaled internally (Default)
    nofinitecheck::Cint  # Disable check for infinite quantities from user?
                         # 0 = do not perform check (Default)
                         # 1 = perform check
    iterproc::Cint       # Placeholder pointer - must set to 0

    Config() = new(1e-10, 1e-10, 1e-10, eps(), 100, 1e-14, 200, 0, 1, 0, 0, 0)
end

#---------------------------------------------------------------------
"CMPFit return structure (C side)."
mutable struct Result_C
    bestnorm::Cdouble           # Final chi^2
    orignorm::Cdouble           # Starting value of chi^2
    niter::Cint                 # Number of iterations
    nfev::Cint                  # Number of function evaluations
    status::Cint                # Fitting status code

    npar::Cint                  # Total number of parameters
    nfree::Cint                 # Number of free parameters
    npegged::Cint               # Number of pegged parameters
    nfunc::Cint                 # Number of residuals (= num. of data points)

    resid::Ptr{Cdouble}         # Final residuals nfunc-vector, or 0 if not desired
    perror::Ptr{Cdouble}        # Final parameter uncertainties (1-sigma) npar-vector, or 0 if not desired
    covar::Ptr{Cdouble}         # Final parameter covariance matrix npar x npar array, or 0 if not desired
    version::NTuple{20, Cchar}  # `mpfit` version string

    "Create an empty `Result_C` structure."
    Result_C() = new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NTuple{20, Cchar}("\0"^20))
end

#---------------------------------------------------------------------
"CMPFit return structure (Julia side)."
mutable struct Result
    bestnorm::Float64           # Final chi^2
    orignorm::Float64           # Starting value of chi^2
    niter::Int64                # Number of iterations
    nfev::Int64                 # Number of function evaluations
    status::Int64               # Fitting status code

    npar::Int64                 # Total number of parameters
    nfree::Int64                # Number of free parameters
    npegged::Int64              # Number of pegged parameters
    nfunc::Int64                # Number of residuals (= num. of data points)

    perror::Vector{Float64}     # Final parameter uncertainties (1-sigma)
    covar::Matrix{Float64}      # Final parameter covariance matrix
    version::String             # CMPFit version string
    param::Vector{Float64}      # Array of best fit parameters
    dof::Int                    # Degrees of freedom
    elapsed::Float64            # Elapsed time [s]
end

#---------------------------------------------------------------------
# Add Base.show method to show CMPFit results
function Base.show(stream::IO, res::Result)
    @printf "  CMPFit ver. = %s\n"  res.version
    @printf "      STATUS = %d\n"   res.status
    @printf "  CHI-SQUARE = %f    (%d DOF)\n"  res.bestnorm  res.dof
    @printf "        NPAR = %d\n"   res.npar
    @printf "       NFREE = %d\n"   res.nfree
    @printf "     NPEGGED = %d\n"   res.npegged
    @printf "       NITER = %d\n"   res.niter
    @printf "        NFEV = %d\n"   res.nfev
    @printf "Elapsed time = %f s\n" res.elapsed

    @printf "\n"
    for i in 1:res.npar
        @printf "  P[%d] = %f +/- %f\n"  i  res.param[i]  res.perror[i]
    end
end

#---------------------------------------------------------------------
mutable struct Wrap_A_Function
    funct::Function
end

#---------------------------------------------------------------------
# This function can not be nested into mpfit since this would lead to the error:
#   ERROR: closures are not yet c-callable
#
"Function called from C to calculate the residuals"
function julia_eval_resid(ndata::Cint, npar::Cint, _param::Ptr{Cdouble}, _resid::Ptr{Cdouble}, _derivs::Ptr{Ptr{Cdouble}}, _funct::Ptr{Wrap_A_Function})
    wrap  = unsafe_load(_funct)
    param = unsafe_wrap(Vector{Float64}, _param, npar)
    resid = unsafe_wrap(Vector{Float64}, _resid, ndata)

    pderiv = Vector{Int}()
    vderiv = Vector{Vector{Float64}}()
    if _derivs != C_NULL
        derivs = unsafe_wrap(Vector{Ptr{Cdouble}}, _derivs, npar)
        for ipar in 1:npar
            if derivs[ipar] != C_NULL
                push!(pderiv, ipar)
                push!(vderiv, unsafe_wrap(Vector{Float64}, derivs[ipar], ndata))
            end
        end
    end

    # Compute residuals
    try
        if length(pderiv) > 0
            jresid = wrap.funct(param, pderiv, vderiv)
        else
            jresid = wrap.funct(param)
        end
        resid .= reshape(jresid, length(jresid))
    catch err
        println("An error occurred during model evaluation: ")
        println(err)
        return Cint(-1)::Cint
    end

    return Cint(0)::Cint
end

#C-compatible address of the Julia `julia_eval_resid` function.
const c_eval_resid = @cfunction(julia_eval_resid, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Ptr{Cdouble}}, Ptr{Wrap_A_Function}))

######################################################################
# Public functions
######################################################################

#---------------------------------------------------------------------
"Main CMPFit function"
function cmpfit(funct::Function,
                _param::Vector{Float64};
                parinfo=nothing, config=nothing)

    # Compute elapsed time
    elapsedtime = Base.time_ns()

    # Use a local copy
    param = deepcopy(_param)

    # Check user function by evaluating it
    model = funct(param)

    res_C = Result_C()
    res_C.perror = Libc.malloc(length(param)   * sizeof(Cdouble) )
    res_C.covar  = Libc.malloc(length(param)^2 * sizeof(Cdouble) )

    if parinfo == nothing
        parinfo = Parinfo(length(param))
    end
    if length(parinfo) != length(param)
        error("Lengths of `param` and `parinfo` arrays are different")
    end
    imm_parinfo = map(imm_Parinfo, parinfo)

    if config == nothing
        config = Config()
    end

    wrap = Wrap_A_Function(funct)
    status = -999
    try
        status = ccall((:mpfit, libmpfit), Cint,
                       (Ptr{Nothing}, Cint         , Cint,          Ptr{Cdouble}, Ptr{imm_Parinfo}, Ptr{Config}, Ptr{Wrap_A_Function}, Ref{Result_C}),
                       c_eval_resid , length(model), length(param), param       , imm_parinfo     , Ref(config), Ref(wrap)           , res_C)
    catch err
        printstyled(color=:red, "An error occurred during `mpfit` call:\n")
        println(err)
        println("")
    end

    perror = Vector{Float64}();
    covar  = Array{Float64, 2}(undef, 0, 0)

    if status > 0
        covar = Vector{Float64}()
        for i in 1:length(param)  ; push!(perror, unsafe_load(res_C.perror, i)); end
        for i in 1:length(param)^2; push!(covar , unsafe_load(res_C.covar , i)); end
        covar  = reshape(covar, length(param), length(param))
    else
        @warn "mpfit returned status = " * string(status) * " < 0."
    end

    version = findall(x -> x != 0, collect(res_C.version))
    if length(version) == 0
        version = ""
    else
        version = join(Char.(res_C.version[version]))
    end

    result = Result(res_C.bestnorm,
                    res_C.orignorm,
                    res_C.niter,
                    res_C.nfev,
                    res_C.status,
                    res_C.npar,
                    res_C.nfree,
                    res_C.npegged,
                    res_C.nfunc,
                    perror,
                    covar,
                    version,
                    param,
                    res_C.nfunc - res_C.nfree,
                    0.)

    Libc.free(res_C.perror)
    Libc.free(res_C.covar )

    elapsedtime = Base.time_ns() - elapsedtime
    result.elapsed = convert(Float64, elapsedtime / 1e9)

    return result
end

#---------------------------------------------------------------------
function cmpfit(independentData::AbstractArray,
                observedData::AbstractArray,
                uncertainties::AbstractArray,
                funct::Function,
                guessParam::Vector{Float64};
                parinfo=nothing, config=nothing)
    function cmpfit_callback(param::Vector{Float64})
        model = funct(independentData, param)
        ret = (observedData - model) ./ uncertainties
        return ret
    end

    function cmpfit_callback(param::Vector{Float64}, pderiv::Vector{Int}, vderiv::Vector{Vector{Float64}})
        model = funct(independentData, param, pderiv, vderiv)
        ret = (observedData - model) ./ uncertainties
        for i in 1:length(vderiv)
            vderiv[i] ./= -uncertainties
        end
        return ret
    end

    cmpfit(cmpfit_callback, guessParam, parinfo=parinfo, config=config)
end

end # module


text:
using Bridge, Distributions, StaticArrays
using BenchmarkTools

suite = BenchmarkGroup()

suite["solver"] = BenchmarkGroup(["Float64"])
suite["solverSA"] = BenchmarkGroup(["SVector", "SMatrix"])

if isdefined(:NEW) && NEW
    ;
else
    BenchmarkTools.loadparams!(suite, BenchmarkTools.load(joinpath("test","perf","params.jld"), "suite"), :evals, :samples);
end

B = -0.5
mu = 0.2
sigma = 0.8
a = sigma*sigma'

P = Bridge.LinPro(B, mu, sigma)

u = 0.5
v = 1.0

t = 0.5
T = 2.0
n2 = 150
tt = range(t, stop=T, length=n2)
K = SamplePath(tt, zeros(length(tt)))
V = SamplePath(tt, zeros(length(tt)))
Mu = SamplePath(tt, zeros(length(tt)))
suite["solver"]["gpK!"] = @benchmarkable Bridge.gpHinv!(K, P)
suite["solver"]["gpV!"] = @benchmarkable Bridge.gpV!(V, v, P)
suite["solver"]["solve!(::R3, ...)"] = @benchmarkable Bridge.solve!(Bridge.R3(), Bridge._F, Mu, u, P)
suite["solver"]["solve!(::BS3, ...)"] = @benchmarkable Bridge.solve!(Bridge.BS3(), Bridge._F, Mu, u, P)

W = sample(tt, Wiener())
X = euler(u, W, P)
suite["solver"]["sample!(..., ::Wiener)"] = @benchmarkable sample!(W, Wiener())
suite["solver"]["euler!(..., ::Wiener"] = @benchmarkable euler!(X, u, W, Wiener())
suite["solver"]["euler!(..., ::LinPro)"] = @benchmarkable euler!(X, u, W, P);

if isdefined(:NEW) && NEW
    tune!(suite["solver"]);
    BenchmarkTools.save(joinpath("test","perf","params.jld"), "suite", BenchmarkTools.params(suite))
else
    BenchmarkTools.loadparams!(suite, BenchmarkTools.load(joinpath("test","perf","params.jld"), "suite"), :evals, :samples);
end

results = run(suite["solver"], verbose = true, seconds = 5) #  _ seconds per benchmark

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

SV = SVector{2,Float64}
SM = SArray{Tuple{2,2},Float64,2,4}
B = SM([-1 0.1; -0.2 -1])
mu = SV([0.2, 0.3])
sigma = SM(2*[-0.212887  0.0687025;
  0.193157  0.388997 ])
a = sigma*sigma'

P = Bridge.LinPro(B, mu, sigma)

u = SV([1.0, 0.0])
v = SV([.5, 0.0])

t = 0.5
T = 2.0
n2 = 150
tt = range(t, stop=T, length=n2)
K = SamplePath(tt, zeros(SM, length(tt)))
V = SamplePath(tt, zeros(SV, length(tt)))
Mu = SamplePath(tt, zeros(SV, length(tt)))
suite["solverSA"]["gpK!"] = @benchmarkable Bridge.gpHinv!(K, P)
suite["solverSA"]["gpV!"] = @benchmarkable Bridge.gpV!(V, v, P)
suite["solverSA"]["solve!(::R3, ...)"] = @benchmarkable Bridge.solve!(Bridge.R3(), Bridge._F, Mu, u, P)
suite["solverSA"]["solve!(::BS3, ...)"] = @benchmarkable Bridge.solve!(Bridge.BS3(), Bridge._F, Mu, u, P)

W = sample(tt, Wiener{SV}())
X = euler(u, W, P)
suite["solverSA"]["sample!(..., ::Wiener)"] = @benchmarkable sample!(W, Wiener{SV}())
suite["solverSA"]["euler!(..., ::Wiener"] = @benchmarkable euler!(X, u, W, Wiener{SV}())
suite["solverSA"]["euler!(..., ::LinPro)"] = @benchmarkable euler!(X, u, W, P);

if isdefined(:NEW) && NEW
    tune!(suite["solverSA"]);
end

resultsSA = run(suite["solverSA"], verbose = true, seconds = 5) #  _ seconds per benchmark

if isdefined(:NEW) && NEW
    BenchmarkTools.save(joinpath("test","perf","params.jld"), "suite", BenchmarkTools.params(suite))
end

BenchmarkTools.save(joinpath("test","perf","results$(Base.Dates.today()).jld"), "oldresults", results, "oldresultsSA", resultsSA)

oldresults, oldresultsSA = BenchmarkTools.load(joinpath("test","perf","results.jld"), "oldresults", "oldresultsSA")
display(median(results))
display(median(resultsSA))
display(judge(median(results), median(oldresults),time_tolerance = 0.01))
display(judge(median(resultsSA), median(oldresultsSA),time_tolerance = 0.01))

#  judge(memory(oldresults["solver"]), memory(results["solver"]))


text:
using JuMP
const MOIT = MOI.Test

MOIU.@model(NoFreeVariable,
            (), (MOI.EqualTo, MOI.LessThan, MOI.GreaterThan), (MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros), (),
            (), (MOI.ScalarAffineFunction,), (MOI.VectorOfVariables,), (MOI.VectorAffineFunction,))
# No free variables to make sure variable bridges are used to increase coverage
MOI.supports_add_constrained_variables(::NoFreeVariable, ::Type{MOI.Reals}) = false
MOI.supports_add_constrained_variables(::NoFreeVariable, ::Type{MOI.Nonnegatives}) = true
MOI.supports_add_constrained_variables(::NoFreeVariable, ::Type{MOI.Nonpositives}) = true
MOI.supports_add_constrained_variables(::NoFreeVariable, ::Type{MOI.Zeros}) = true

function bridged_mock(mock_optimize!::Function...;
                      model = MOIU.Model{Float64}())
    mock = MOI.Utilities.MockOptimizer(model)
    bridged = MOI.Bridges.full_bridge_optimizer(mock, Float64)
    MOI.Utilities.set_mock_optimize!(mock, mock_optimize!...)
    return bridged
end


text:
@doc """
Find all dipole in an activity map.

Determines the local maxima in a 3 dimensional array

#### Input

* s: Activity in 3d matrix
* window: Windowing to use in each dimension for min max filter
* x,y,z: Coordinates associated with s matrix

#### Output

* dips: An array of dipoles

""" ->
function find_dipoles{T <: Number}(s::Array{T, 3}; window::Array{Int}=[6,6,6], x::AbstractVector{T}=1:size(s,1),
                                   y::AbstractVector{T}=1:size(s,2), z::AbstractVector{T}=1:size(s,3))

    debug("Finding dipoles for 3d array")

    minval, maxval = extrema_filter(s, window)

    # Find the positions matching the maxima
    matching = s[2:size(maxval)[1]+1, 2:size(maxval)[2]+1, 2:size(maxval)[3]+1]
    matching = matching .== maxval

    # dipoles are defined as maxima locations and within 90% of the maximum
    peaks = maxval[matching]
    peaks = peaks[peaks .>= 0.1 * maximum(peaks)]

    # Store dipoles in an array
    dips = Dipole[]

    for l = 1:length(peaks)
        xidx, yidx, zidx = ind2sub(size(s), find(s .== peaks[l]))
        push!(dips, Dipole("Unknown", x[xidx[1]], y[yidx[1]], z[zidx[1]], 0, 0, 0, 0, 0, peaks[l]))
    end

    left_side = s[x .< 0, :, :]
    xidx, yidx, zidx = ind2sub(size(left_side), find(left_side .== maximum(left_side)))
    push!(dips, Dipole("Unknown", x[xidx[1]], y[yidx[1]], z[zidx[1]], 0, 0, 0, 0, 0, maximum(left_side)))

    right_side = s[x .> 0, :, :]
    xidx, yidx, zidx = ind2sub(size(right_side), find(right_side .== maximum(right_side)))
    x_tmp = x[x .> 0]
    push!(dips, Dipole("Unknown", x_tmp[xidx[1]], y[yidx[1]], z[zidx[1]], 0, 0, 0, 0, 0, maximum(right_side)))

    # Sort dipoles by size
    vec(dips[sortperm([dip.size for dip in dips], rev=true)])
end

#######################################
#
# Find the best dipoles from selection
#
#######################################

@doc """
Find best dipole relative to reference location.

Finds the largest dipole within a specified distance of a reference location

#### Input

* ref: Reference coordinate or dipole
* dips: Dipoles to find the best dipole from
* maxdist: Maximum distance a dipole can be from the reference

#### Output

* dip: The best dipole

""" ->
function best_dipole(ref::Union{Coordinate, Dipole}, dips::Array{Dipole}; maxdist::Number=0.30, min_dipole_size::Real=-Inf, kwargs...)

    Logging.info("Calculating best dipole for $(length(dips)) dipoles")

    dips = dips[find([d.size > min_dipole_size for d in dips])]

    if length(dips) > 0

      # Find all dipoles within distance
      dists = [euclidean(ref, dip) for dip=dips]
      valid_dist = dists .< maxdist

      if sum(valid_dist) >= 2
          # Valid dipoles exist find the largest one
          sizes = [dip.size for dip =dips]
          bestdip = maximum(sizes[valid_dist])
          dip = dips[find(sizes .== bestdip)]
          debug("$(sum(valid_dist)) dipoles within $(maxdist) m. ")

      elseif sum(valid_dist) == 1
          # Return the one valid dipole
          dip = dips[find(valid_dist)]
          debug("Only one dipole within $(maxdist) m. ")

      else
          # No dipoles within distance
          # Take the closest
          bestdip = minimum(dists)
          dip = dips[find(dists .== bestdip)]
          debug("No dipole within $(maxdist) m. ")

      end
      debug("Best = $(euclidean(ref, dip[1]))")

      return dip[1]

    else
      return NaN
    end
end


text:
using Plots
using SpikingNeuralNetworks
SNN.@load_units

E = SNN.HH(;N = 3200)
I = SNN.HH(;N = 800)
EE = SNN.SpikingSynapse(E, E, :ge; σ = 6nS, p = 0.02)
EI = SNN.SpikingSynapse(E, I, :ge; σ = 6nS, p = 0.02)
IE = SNN.SpikingSynapse(I, E, :gi; σ = 67nS, p = 0.02)
II = SNN.SpikingSynapse(I, I, :gi; σ = 67nS, p = 0.02)
P = [E, I]
C = [EE, EI, IE, II]

SNN.monitor(E, [(:v, [1, 10, 100])])
SNN.sim!(P, C; dt = 0.01ms, duration = 100ms)
SNN.vecplot(E, :v)


text:
# Autogenerated wrapper script for CFITSIO_jll for armv7l-linux-gnueabihf
export libcfitsio

using LibCURL_jll
JLLWrappers.@generate_wrapper_header("CFITSIO")
JLLWrappers.@declare_library_product(libcfitsio, "libcfitsio.so.9")
function __init__()
    JLLWrappers.@generate_init_header(LibCURL_jll)
    JLLWrappers.@init_library_product(
        libcfitsio,
        "lib/libcfitsio.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
fld-write.jl
write AVS .fld to file
2019-05-12, Jeff Fessler, University of Michigan
=#

export fld_write

"""
`fld_write(file, data, varargin)`

write data into AVS format `.fld` file

See ASPIRE user's guide (or AVS online documentation) for file format.

in
- `file` name of file typically ending in `.fld`
- `data` real data array

option
- `check::Bool`			report error if file exists; default `true`
- `dir::String`			directory name to prepend file name; default `""`
- `endian::`Symbol`		`:le` little endian (default), `:be` big endian
- `head::Array{String}`	comment information for file header
- `raw::Bool`			put raw data in `name.raw`, header in `name.fld`
						where `file` = `name.fld`; default `false`
"""

# * `how::Symbol`			file data type; default `:data` infer from data

function fld_write(
	file::String,
	data::AbstractArray{<:Real} ;
	check::Bool = true,
	dir::String = "",
	endian::Symbol = :le,
	head::Array{String} = empty([""]),
#	how::Symbol = :data,
	raw::Bool = false,
)

	data = fld_write_data_fix(data) # make data suitable for writing in .fld

	endian != :le && endian != :be && throw("endian '$endian' unknown")

	typedict = Dict([
		(Float32, endian === :le ? "float_le" : "xdr_float"),
		(Float64, endian === :le ? "double_le" : "xdr_double"),
		(UInt8, "byte"),
		(Int16, endian === :le ? "short_le" : "short_be"),
		(Int32, endian === :le ? "int_le" : "xdr_int"),
		])

	datatype = typedict[eltype(data)] # throws error if not there

	file = joinpath(dir, file)
#	@show file

	check && isfile(file) && throw("file '$file' exists")

	if raw # if writing data to separate ".raw" file, ensure it does not exist
		fileraw = file[1:(end-3)] * "raw"
		check && isfile(fileraw) && throw("file $fileraw exists")
	end

	# open output avs file for writing
	fid = open(file, "w")
	fraw = raw ? open(fileraw, "w") : fid

	# write header
	ndim = ndims(data)

	println(fid, "# AVS field file ($(basename(@__FILE__)))")
	for line in head
		println(fid, "# $line")
	end

	println(fid, "ndim=$ndim")
	for ii=1:ndim
		println(fid, "dim$ii=$(size(data,ii))")
	end
	println(fid, "nspace=$ndim")
	println(fid, "veclen=1")
	println(fid, "data=$datatype")
	println(fid, "field=uniform")

	if raw
		println(fid, "variable 1 file=$fileraw filetype=binary")
	else
		write(fid, "\f\f") # two form feeds: char(12)
	end

	# finally, write the binary data
	host_is_le = () -> ENDIAN_BOM == 0x04030201
	fun = (host_is_le() == (endian === :le)) ?
		identity : # host/file same endian
		(host_is_le() && (endian === :be)) ?
		hton :
		(!host_is_le() && (endian === :le)) ?
		write(fraw, htol.(data)) :
		throw("not done")
	write(fraw, fun.(data))

	close(fid)

	if raw
		close(fraw)
	end

end

"""
`data = fld_write_data_fix(data)`
convert data to format suitable for writing to .fld file
"""
function fld_write_data_fix(data::AbstractArray{BigFloat})
	@warn "BigFloat downgraded to Float64"
	return Float64.(data)
end
function fld_write_data_fix(data::AbstractArray{Float16})
	@warn "Float16 promoted to Float32"
	return Float32.(data)
end
function fld_write_data_fix(data::AbstractArray{T}) where {T <: Union{BigInt, Int64}}
	@warn "$T downgraded to Int32"
	return Int32.(data)
end
function fld_write_data_fix(data::AbstractArray{Bool})
	@warn "Bool promoted to UInt8"
	return UInt8.(data)
end
@inline fld_write_data_fix(data::AbstractArray{T}) where
	{T <: Union{Float32, Float64, UInt8, Int16, Int32}} = data
fld_write_data_fix(data::AbstractArray{T}) where {T <: Any} =
		throw(ArgumentError("unsupported type $T"))


text:
abstract type AbstractNNPolicy  <: Policy end

## NN Policy interface

"""
    getnetwork(policy)
    return the  value network of the policy
"""
function getnetwork end

"""
    resetstate!(policy)
reset the hidden states of a policy
"""
function resetstate! end

struct NNPolicy{P, Q, A} <: AbstractNNPolicy
    problem::P
    qnetwork::Q
    action_map::Vector{A}
    n_input_dims::Int64
end

NNPolicy(problem::MDPCommonRLEnv, qnetwork, action_map::Vector, n_input_dims::Int) = NNPolicy(convert(MDP, problem), qnetwork, action_map, n_input_dims)
NNPolicy(problem::POMDPCommonRLEnv, qnetwork, action_map::Vector, n_input_dims::Int) = NNPolicy(convert(POMDP, problem), qnetwork, action_map, n_input_dims)

function getnetwork(policy::NNPolicy)
    return policy.qnetwork
end

function resetstate!(policy::NNPolicy)
    Flux.reset!(policy.qnetwork)
end

actionmap(p::NNPolicy) = p.action_map

function _action(policy::NNPolicy, o)
    if ndims(o) == policy.n_input_dims
        obatch = reshape(o, (size(o)...,1))
        vals = policy.qnetwork(obatch)
        return policy.action_map[argmax(vals)]
    else
        throw("NNPolicyError: was expecting an array with $(policy.n_input_dims) dimensions, got $(ndims(o))")
    end
end

function _actionvalues(policy::NNPolicy{P,Q,A}, o::AbstractArray{T,N}) where {P,Q,A,T<:Real,N}
    if ndims(o) == policy.n_input_dims
        obatch = reshape(o, (size(o)...,1))
        return dropdims(policy.qnetwork(obatch), dims=2)
    else
        throw("NNPolicyError: was expecting an array with $(policy.n_input_dims) dimensions, got $(ndims(o))")
    end
end

function _value(policy::NNPolicy{P}, o::AbstractArray{T,N}) where {P,T<:Real,N}
    if ndims(o) == policy.n_input_dims
        obatch = reshape(o, (size(o)...,1))
        return maximum(policy.qnetwork(obatch))
    else
        throw("NNPolicyError: was expecting an array with $(policy.n_input_dims) dimensions, got $(ndims(o))")
    end
end

POMDPs.action(policy::NNPolicy, o) = _action(policy, o)
POMDPs.action(policy::NNPolicy{P}, s) where {P <: MDP} = _action(policy, POMDPs.convert_s(Vector{Float32}, s, policy.problem))
POMDPs.action(policy::NNPolicy{P}, a, o) where {P <: POMDP} = _action(policy, vcat(vec(POMDPs.convert_a(Vector{Float32}, a, policy.problem)), vec(POMDPs.convert_o(Vector{Float32}, o, policy.problem))))

POMDPPolicies.actionvalues(policy::NNPolicy, o) = _actionvalues(policy, o)
POMDPPolicies.actionvalues(policy::NNPolicy{P}, s) where {P<:MDP} = _actionvalues(policy, POMDPs.convert_s(Vector{Float32}, s, policy.problem))
POMDPPolicies.actionvalues(policy::NNPolicy{P}, a, o) where {P<:POMDP} = _actionvalues(policy, vcat(vec(POMDPs.convert_a(Vector{Float32}, a, policy.problem)), vec(POMDPs.convert_o(Vector{Float32}, o, policy.problem))))

POMDPs.value(policy::NNPolicy, o) = _value(policy, o)
POMDPs.value(policy::NNPolicy{P}, s) where {P <: MDP} = _value(policy, POMDPs.convert_s(Vector{Float32}, s, policy.problem))
POMDPs.value(policy::NNPolicy{P}, a, o) where {P <: POMDP} = _value(policy, vcat(vec(POMDPs.convert_a(Vector{Float32}, a, policy.problem)), vec(POMDPs.convert_o(Vector{Float32}, o, policy.problem))))


text:
using Merlin
using Merlin.Datasets.MNIST

const NEPOCHS = 50

function main()
    trainmode = true
    datapath = joinpath(@__DIR__, ".data")
    traindata = setup_data(MNIST.traindata(datapath)...)
    testdata = setup_data(MNIST.testdata(datapath)...)
    savefile = "mnist_epoch$(NEPOCHS).jld2"
    if trainmode
        # model = setup_model()
        train(traindata, testdata)
        #save(savefile, "model", model)
    else
        model = load(savefile, "model")
        test(model, testdata)
    end
end

function setup_data(x::Matrix{Float32}, y::Vector{Int})
    batchsize = 200
    xs = [x[:,i:i+batchsize-1] for i=1:batchsize:size(x,2)]
    xs = map(BACKEND, xs)
    y += 1 # Change label set: 0..9 -> 1..10
    y = Array{Int32}(y)
    ys = [y[i:i+batchsize-1] for i=1:batchsize:length(y)]
    ys = map(BACKEND, ys)
    collect(zip(xs,ys))
end

function NN()
    T = Float32
    hsize = 1000
    h = Linear(T,28*28,hsize)(Node(name="x"))
    h = relu(h)
    h = Linear(T,hsize,hsize)(h)
    h = relu(h)
    h = Linear(T,hsize,10)(h)
    BACKEND(Graph(h))
end

function train(traindata::Vector, testdata::Vector)
    nn = NN()
    opt = SGD(0.001)

    for epoch = 1:NEPOCHS
        println("epoch: $epoch")
        prog = Progress(length(traindata))
        loss = 0.0

        for (x,y) in shuffle!(traindata)
            z = nn(Var(x))
            z = softmax_crossentropy(Var(y), z)
            loss += sum(z.data)
            params = gradient!(z)
            foreach(opt, params)
            ProgressMeter.next!(prog)
        end
        loss /= length(traindata)
        println("Loss:\t$loss")

        golds = Int32[]
        preds = Int32[]
        for (x,y) in testdata
            append!(golds, Array(y))
            z = nn(Var(x))
            z = argmax(Array(z.data), 1)
            append!(preds, z)
        end
        @assert length(golds) == length(preds)
        acc = mean(i -> golds[i] == preds[i] ? 1.0 : 0.0, 1:length(golds))
        println("test accuracy: $acc")
        println()
    end
end

@time main()


text:
using TropicalYao.Reversible
using Test

@testset "LogLikeNumbers" begin
    include("LogLikeNumbers/LogLikeNumbers.jl")
end

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

@testset "instructions" begin
    include("instructs.jl")
end


text:
using Gtk

function controlwindow(win, lab)
    sleep(4)
    set_gtk_property!(lab, :label, "Hiding...")
    sleep(1)
    println("Hiding widow")
    set_gtk_property!(win, :visible, false)
    sleep(5)
    set_gtk_property!(lab, :label, "Showing...")
    println("Showing window")
    set_gtk_property!(win, :visible, true)
    sleep(5)
    set_gtk_property!(lab, :label, "Resizing...")
    println("Resizing window")
    resize!(win, 300, 300)
    sleep(4)
    set_gtk_property!(lab, :label, "Maximizing...")
    println("Maximizing window")
    sleep(1)
    maximize(win)
    set_gtk_property!(lab, :label, "Closing...")
    sleep(5)
    println("Closing window")
    destroy(win)
    sleep(2)
    exit(0)
end

function runwindow()
    win = GtkWindow("Window Control Test", 500, 30) |> (GtkFrame() |> (vbox = GtkBox(:v)))
    lab = GtkLabel("Window under external control")
    push!(vbox, lab)
    @async(controlwindow(win, lab))

    cond = Condition()
    endit(w) = notify(cond)
    signal_connect(endit, win, :destroy)
    showall(win)
    wait(cond)
end

runwindow()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Probably use FunctionWrappers.jl for interactive mode, to avoid
# excessive recompilation.

struct ProfileWrapper{P,F,T1,T2,R1,R2}#,C<:AbstractArray}
    f::F
    x::MutableFixedSizeVector{P,T1,R1,R1}
    y::MutableFixedSizeVector{P,T2,R2,R2}
    # z::C
    i::Base.RefValue{Int}
    v::Base.RefValue{T1}
end
function (pw::ProfileWrapper{P,F,T})(x::PaddedMatrices.AbstractMutableFixedSizeVector{Pm1,T}) where {P,Pm1,F,T}
    debug() && @show x
    @inbounds begin
        for i in 1:pw.i[]-1
            pw.x[i] = x[i]
        end
        # pw.x[pw.i[]] = pw.v[]
        for i in pw.i[]:Pm1
            pw.x[i+1] = x[i]
        end
    end
    debug() && @show pw.x
    - pw.f(pw.x)
end

function (pw::ProfileWrapper{P,F,T,T2})(y::PaddedMatrices.AbstractMutableFixedSizeVector{Pm1,T2}) where {P,Pm1,F,T,T2}
    # @show y
    @inbounds begin
        for i in 1:pw.i[]-1
            pw.y[i] = y[i]
        end
        # pw.x[pw.i[]] = pw.v[]
        for i in pw.i[]:Pm1
            pw.y[i+1] = y[i]
        end
    end
    # @show pw.y
    - pw.f(pw.y)
end

function set_profile_val_ind!(pw::ProfileWrapper, v, i::Int)
    pw.v[] = v
    pw.i[] = i
    @inbounds pw.x[i] = v
    @inbounds pw.y[i] = v
    nothing
end

struct Swap{P,F}
    f::F
    i::Base.RefValue{Int}
end
function (s::Swap{P})(x) where P
    @inbounds x[s.i[]], x[P] = x[P], x[s.i[]]
    - s.f(x)
end
Swap(f::F, ::Val{P}) where {F,P} = Swap{P,F}(f,Ref(P))

function swappable_twice_differentiable(f, ::Val{P}) where {P}
    TwiceDifferentiable(Swap(f, Val{P}()), Val{P}())
end

# struct Swappable{P,D <: DifferentiableObjects.AbstractTwiceDifferentiableObject{P}}
#     obj::D
#     i::Base.RefValue{Int}
# end


text:
module OptimNewton

export get_root;
export optimize;

"""
	Try to find a root of (fun). You need to provide the first derivative of fun (fder) and an initial x0 value.	Additional parameters are used as stop signals.	xtol, ytol: are the tolerances for x and f(x) respectively, if further iteration do not	improve the result in more that this tolerance values, the method returns. itn: is the max number of iterations that the method will do. If verbose is true, the method will print a status message in every iteration.
"""
function get_root(x0, fun, fder; xtol = 10e-10, ytol = 10e-10,
 itn = 10, verbose = false)

	r = x0;
	lastr = r;
	for i in 1:itn
		verbose && print("iteration $i: initial root $r,");
		r = r - fun(r)/fder(r);
		verbose && println(" oprox root $r");
		abs(lastr - r) <= xtol && break;
		lastr = r;
	end

	return r;

end # get_root

"""
	Try to find a optimus value of (fun). You need to provide the first derivative of fun (fder), de secund derivative of fun (sder) and an initial x0 value. Additional parameters are used as stop signals. xtol, ytol: are the tolerances for x and f(x) respectively, if further iteration do not	improve the result in more that this tolerance values, the method return. itn: is the max number of iterations that the method will do. If verbose is true, the method will print a status message in every iteration.
"""
function optimize(x0, fun, fder, sder; xtol = 10e-10, ytol = 10e-10, itn = 10, verbose = false)

	r = x0;
	lastr = r;
	for i in 1:itn
		verbose && print("iteration $i: before value $r,");
		r = r - fder(r)/sder(r);
		verbose && println(" after value $r");
		abs(lastr - r) <= xtol && break;
		abs(fun(r) - fun(lastr)) <= ytol && break;
		lastr = r;
	end

	return r;

end # optimize

end # module


text:
using DifferentialEquations, Plots
import StatsBase, JSON, Random

"""
This file contains the definition of the Adaptive Filter model as well as
the Lazar model as an ODE to be solved by the DifferentialEquations library.
"""

#=== meta execution ========#
exec = false            # if false, no solving or plotting from this file
algo = :naf             # options: :af, :naf, :lazar
task = "relu"           # options: recovery:     "sin", "mix2", "mix3",
                        # linear computation:    "1of3", "shift", "noise"
                        # nonlinear computation: "abs", "relu", "sign"
mode = "train"          # options: "train", "test" (disables learning rule)
show = "mse_evol"       # options: "out", "y", "y*w", "zyw", "imp", "impi",
                        # "w_evol", "mse_evol" (measured during training)
mark = ""               # any mark to put in the figure filename

#=== input/output setup ====#
lp = 3                  # lowest = shortest period (fastest frequency)
mp = 4                  # middle period (for mix3 and 1of3 task)
hp = 5                  # highest = longest period (slowest frequency)
phase = lp/4            # phase for task="shift" (additional to latency)
offset = 0              # y-offset of signal (important for nonlinearity)
latency = -.60          # accounts for latency of information due to spikes
noise_std = .5          # std for gaussian white noise

#=== encoder parameters ====#
rate = 10               # the resting firing rate (spikes/sec)
b = 3.0                 # bias (> |minimal input|, s.t. x(t) + b > 0)
θ = 1.0                 # threshold (when reached, we spike)
κ = (b+offset)/rate/θ   # capacitance (integrator weight)
Δ = 0.0                 # TODO refractory?
l = 0.0                 # TODO leak?

#=== decoder parameters ====#
n = 10                  # n filters in total (equal for both models)
η = 1e-4                # learning rate      (equal for both models)
α = 2                   # fastest decay rate
λ = 2                   # fastest decay rate
λs = range(1, length=n, stop=0.01)

ramp(z) = log(1 + ℯ^z)  # non-linearity (ramp)
sigmoid(z) = ℯ^z / (1 + ℯ^z) # derivative of ramp
sig_der(z) = ℯ^(-z) / (1 + ℯ^(-z))^2 # derivative of sigmoid

#=== simulation ============#
training_cycles = 10000  # how many total periods to see in training
test_cycles = 3          # and how many in the test phase

function init_task(task)
    """
    Returns an io tuple (of input x and target z) and the target period.
    The target period is just the lowest period (lp) for a sine-based task
    and the least common multiple of the periods in a mix task.
    """
    Random.Random.seed!(1) # rng seed, for model comparability

    io = () # a tuple of functions x (input) and z (target)

    if task == "sin"       # reconstruct sine
        let input
            input(t) = sin(2π/lp*t) + offset
            io = (x=input, z=t->input(t+latency))
        end
        tp = lp
    elseif task == "mix2"  # reconstruct a mix of 2 sines
        let input
            input(t) = sin(2π/hp*t) + sin(2π/lp*t) + offset
            io = (x=input,  z=t->input(t+latency))
        end
        tp = lcm([lp, hp]) # total period i
    elseif task == "mix3"  # reconstruct a mix of 3 sines
        let input
            input(t) = sin(2π/hp*t) + sin(2π/mp*t) + sin(2π/lp*t) + offset
            io = (x=input,  z=t->input(t+latency))
        end
        tp = lcm([lp, mp, hp]) # total period is lcm of periods
    elseif task == "1of3"  # demixing: extract 1 component from a mix3
        let input, target
            input(t) = sin(2π/hp*t) + sin(2π/mp*t) + sin(2π/lp*t) + offset
            target(t) = sin(2π/mp*t) + offset
            io = (x=input, z=t->target(t+latency))
        end
        tp = mp
    elseif task == "shift" # shift sine by phase
        let input
            input(t) = sin(2π/lp*t) + offset
            io = (x=input, z=t->input(t+latency+phase))
        end
        tp = lp
    elseif task == "noise" # denoise: reconstruct sine with white noise
        let input
            input(t) = sin(2π/lp*t) + offset
            noise = randn(1000) .* noise_std # noise with µ=0
            io = (x=t->input(t) + noise[Int(floor(t*100))%1000+1],
                  z=t->input(t+latency))
        end
        tp = lp
    elseif task == "abs"   # take absolute value
        let input
            input(t) = sin(2π/lp*t) + offset
            io = (x=input, z=t->abs(input(t+latency)-offset)+offset)
        end
        tp = lp
    elseif task == "relu"  # take maximum of 0 and signal
        let input
            input(t) = sin(2π/lp*t) + offset
            io = (x=input, z=t->Base.max(0., input(t+latency)-offset)+offset)
        end
        tp = lp
    elseif task == "sqr"   # square of input
        let input
            input(t) = sin(2π/lp*t) + offset
            io = (x=input, z=t->(input(t+latency)-offset)^2+offset)
        end
        tp = lp
    elseif task == "sign"  # binary transformation: sign of input
        let input
            input(t) = sin(2π/lp*t) + offset
            io = (x=input, z=t->sign(input(t+latency)-offset)+offset)
        end
        tp = lp
    end

    global tp = tp
    return io, tp
end

function init_times(cycles, tp)
    """
    `Cycles` is how many total periods (`tp`) of the signal are presented.
    According to the `show` variable, the saving times (`tsave`) are set.
    Returns the list of save points `tsave`.
    """
    tmax = tp * cycles

    if show in ["w_evol", "mse_evol"]
        tsave = 0.:tp:tmax    # save regularly (once / period)
    else
        tsave = [0.,tmax]     # save only first and last state
    end

    return tsave
end

function init_state(; range=.05)
    """
    Initializes taps and weights randomly from a normal disrtibution
    centered around zero with standard deviation of `range`.
    Note that there are n+1 weights, the last one is the bias weight.
    Returns the initial state being in the form [ys; ws; 0].
    """
    Random.Random.seed!(1)         # rng seed, for comparability
    w0 = (.5 .- rand(n+1))*2*range # random initial weights
    y0 = (.5 .- rand(n))*2*range   # initial taps
    u0 = [y0; w0; 0.0]             # initial state: ys, ws, X
    return u0
end

function init_training(task, cycles)
    """
    Wrapper function around `init_task`, `init_times` and `init_state`.
    """
    # some global settings
    global mode = "train"
    global λs = range(λ, length=n, stop=0.01) # for the lazar decoder

    io, tp = init_task(task)
    u0 = init_state()
    tsave = init_times(cycles, tp)
    return io, tp, u0, tsave
end

#=== ODE DEFINITIONS =======#

function af(du, u, io, t)
    """
    System of ODE to model communication between two AF neurons,
    one reduced on the encoding, the other reduced on the decoding part.

    input signal → enc → channel (spikes) → dec → reconstructed signal

    The state has the following form: (y_1, …, y_n, w_1, …, w_n, w_b, X).
    Note that the spiking part (X) belongs to the presynaptic
    and the filter bank (y_1, …, w_n, w_b) to the postsynaptic neuron.
    """

    # decoder (postsynaptic)
    banklen = n ÷ 2                  # two filter banks of the same size

    # filter bank 1
    du[1] = -α * u[1]                # no input here (generated by events)
    du[2:banklen] .= -α * u[2:banklen] .+ α * u[1:banklen-1]

    # decoder output ("potential")
    z = u[1:n]' * u[n+1:2n] + u[2n+1]

    # filter bank 2
    du[banklen+1] = -α * u[banklen+1] + α * z # input from previous output z
    du[banklen+2:n] .= -α * u[banklen+2:n] .+ α * u[banklen+1:n-1]

    # weight updates
    if mode == "train"
        du[n+1:2n] = 2 * η * (io.z(t) - z) .* u[1:n]
        du[2n+1]   = 2 * η * (io.z(t) - z) # bias weight
    end

    # encoder (presynaptic)
    du[end] = (io.x(t) + b) / κ    # rectangular integration of signal
end

function naf(du, u, io, t)
    """
    Nonlinear AF, same as AF but with ramp function.

    System of ODE to model communication between two AF neurons,
    one reduced on the encoding, the other reduced on the decoding part.

    input signal → enc → channel (spikes) → dec → reconstructed signal

    The state has the following form: (y_1, …, y_n, w_1, …, w_n, w_b, X).
    Note that the spiking part (X) belongs to the presynaptic
    and the filter bank (y_1, …, w_n, w_b) to the postsynaptic neuron.
    """

    # decoder (postsynaptic)
    banklen = n ÷ 2                  # two filter banks of the same size

    # filter bank 1
    du[1] = -α * u[1]                # no input here (generated by events)
    du[2:banklen] .= -α * u[2:banklen] .+ α * u[1:banklen-1]

    # decoder output ("potential")
    z = u[1:n]' * u[n+1:2n] + u[2n+1]

    # filter bank 2
    du[banklen+1] = -α * u[banklen+1] + α*z # input from previous output z
    du[banklen+2:n] .= -α * u[banklen+2:n] .+ α * u[banklen+1:n-1]

    # weight updates
    if mode == "train"
        w_b = u[end-1] # bias weight for shift of ramp
        du[n+1:2n] = (2 * η * (io.z(t) - ramp(z+w_b)) .* u[1:n]
                      * sigmoid(z+w_b))
        # update for w_b (multiplied with n for stronger learning)
        du[end-1]  = 2 * n * η * (io.z(t) - ramp(z+w_b)) * sigmoid(z+w_b)
    end

    # encoder (presynaptic)
    du[end] = (io.x(t) + b) / κ    # rectangular integration of signal
end

function lazar(du, u, io, t)
    """
    Implementation of [Lazar 2006: A simple model for spike processing].
    A zeroth-order approximation of the stimulus is computed
    by spike weighting and filtering the signal.

    A decoding filter bank is constructed from n parallel filters
    each of which gets the same input from the presynaptic neuron.

    Filters have different time constants λ_i stored in λs.
    The output is the weighted sum of the filters, ∑_i w_i · y_i.

    The state has the following form: (y_1, …, y_n, w_1, …, w_n, w_b, X).
    Note that the spiking part (X) belongs to the presynaptic
    and the filter bank (y_1, …, w_n, w_b) to the postsynaptic neuron.
    """

    # decoder (postsynaptic)
    du[1:n] .= -λs .* u[1:n]

    # output
    z = u[1:n]' * u[n+1:2n] + u[2n+1]

    # weight updates
    if mode == "train"
        du[n+1:2n] .= 2 * η * (io.z(t) - z) .* u[1:n]
        du[2n+1]    = 2 * η * (io.z(t) - z) # bias weight
    end

    # encoder (presynaptic)
    du[end] = (io.x(t) + b) / κ # rectangular integration of signal
end

#=== ACCESSOR UTILS ==========#

function z(sol, t)
    u = sol(t)
    y = u[1:n]
    w = u[n+1:2n]
    w_b = u[2n+1] # bias weight
    w' * y + w_b
end

function r(sol, t)
    u = sol(t)
    ramp(z(sol, t))
end

function out(sol, t)
    if algo == :naf
        r(sol, t)
    else
        z(sol, t)
    end
end

function y_i(sol, i)
    function y(t)
        sol(t)[i]
    end
end

function w_i(sol, i)
    function w(t)
        sol(t)[n+i]
    end
end

function mse(sol, io, period)
    """
    Calculate mean squared error (normalized by signal amplitude)
    of the solution output w.r.t. the target function
    over the given period (array of time points)
    and normalize it by the std of the target amplitude.
    """
    mse = sum([(out(sol,t) - io.z(t))^2 for t in period]) /length(period)

    # compute std of target (on 1000 points over tp) to normalize mse
    var = StatsBase.var([io.z(t) for t in period])
    error = mse / var
    round(error, digits=5) # round on 5 decimal digits
end

#=== SOLVER FUNCTIONS ======#

function run(io, u0, tsave, save_spikes=false)
    """
    Runs the model `algo` in `mode` (train or test), i.e.
    it solves the specified ODE over the time specified in `tsave`.
    Returns the solution of the solver and if appropriate
    (`save_spikes == true`), also a list of all spike times.
    """
    tspan = (tsave[1], tsave[end])       # run over the whole save period
    tspike = undef                       # time of last spike, no initial
    spikes = []                          # list to save all spike times
    function threshold(u, t, integrator) # callback condition
        return u[end] - θ
    end
    function spike!(integrator)          # callback handler
        if algo == :lazar
            if tspike != undef           # ignore first spike (weight issue)
                isi = integrator.t - tspike  # inter-spike interval
                c_k = κ * θ - b * isi        # weighting for spikes
                integrator.u[1:n] .+= c_k    # add spike upon each tap
            end
            tspike = integrator.t        # save time of spike
        else                             # for AF/nAF, just do:
            integrator.u[1]  += 1.0      # input to postsynaptic neuron
        end
        integrator.u[end] = 0.0          # reset X integrator
        if save_spikes
            push!(spikes, integrator.t)  # save spike time
        end
    end
    fire = ContinuousCallback(threshold, spike!, # not to save u at spikes:
                              save_positions=(false, false))
    prob = ODEProblem(eval(algo), u0, tspan, io, callback=fire)
    sol = solve(prob, saveat=tsave, maxiters=Inf)
    if save_spikes
        return (sol, spikes) # solution and all spikes
    end
    return sol
end

function test_error(sol, io, tp)
    """
    Tests the last state of a solution `sol`,
    saves for `100/tp * test_cycles` times and
    returns the new solution object and the error.
    """
    global mode = "test"
    u_end, t_end = sol[end], sol.t[end]
    tsave2 = t_end:.01:t_end+test_cycles*tp
    sol2   = run(io, u_end, tsave2)
    return sol2, mse(sol2, io, tsave2)
end

function impulse_response(sol; dur=10)
    """
    Compute impulse response of a trained model (`sol`).
    The model is fully specified by its weight vector.
    Returns a solution for the i.r. over the duration `dur`.
    """
    w_end = sol[end][n+1:2n+1]
    u = zeros(2n+2)       # make new fresh state
    u[n+1:2n+1] = w_end   # take learned model
    if algo == :lazar
        u[1:n] .= 1       # put initial spike into all taps
    else
        u[1] = 1          # initial spike y_1 = 1
    end
    global mode = "test"  # do not adapt weights in rerun
    io = (x=t->0, )       # use no input and no z, because mode == "test"
    # Note that no encoder is necessary, therefore no callback.
    prob = ODEProblem(eval(algo), u, (0.,dur), io)
    solve(prob, saveat=0:.01:dur, maxiters=Inf)
end

function running_mse(sol, io, save_periods=[])
    """
    Takes a solution object and tests the model for each stored state.
    `save_periods` is a list of periods for which the solution is saved.
    Returns a list of MSEs and a list of the saved solutions, if any.
    """
    mses   = []
    sols   = []
    period = 0
    global mode = "test"
    for (u,t) in zip(sol.u, sol.t)
        tsave2 = t:tp/100:t+test_cycles*tp
        sol2 = run(io, u, tsave2)
        push!(mses, mse(sol2, io, tsave2))
        if period in save_periods
            push!(sols, sol2)
        end
        period += 1
    end
    return mses, sols
end

#=== BASIC PLOTS =============#

function palgo()
    """ Prettifies the algo name for plots. """
    if algo == :af
        "AF"
    elseif algo == :lazar
        "Lazar"
    elseif algo == :naf
        "nAF"
    end
end

function evol_plot(show, sol, io)
    """ Returns a plot. """
    if !(show in ["w_evol", "y_evol", "mse_evol"]) return end
    if show == "w_evol"
        plot([w_i(sol, i) for i in 1:n+1], sol.t,
             labels=["w_$i" for i in 1:n+1],
             xlabel="time",
             title="$(palgo()), weights", titleloc=:left)
    elseif show == "y_evol"
        plot([w_i(sol, i) for i in 1:n+1], sol.t,
             labels=["y_$i" for i in 1:n+1],
             xlabel="time",
             title="$(palgo()), taps")
    elseif show == "mse_evol"
        mses, _ = running_mse(sol, io)
        plot(mses,
             leg=false, title="$algo (final MSE: $(mses[end]))",
             ylabel="mse", xlabel="periods")
    end
end

function end_plot(show, sol, io, tp; return_mse=false)
    """ Returns a plot and, if `return_mse==true` also the error. """
    if !(show in ["out", "zyw", "y*w", "y"]) return end

    sol2, error = test_error(sol, io, tp)

    plt = nothing
    if show == "out"
        plt = plot([t->out(sol2, t), io.z], sol2.t,
                   labels=["approx", "target"],
                   title="$algo (MSE: $error)")
    elseif show == "y"
        plt = plot([y_i(sol2, i) for i in 1:n], sol2.t,
                   labels=["y_$i" for i in 1:n], leg=:bottomright,
                   xticks=sol2.t[1]:tp:sol2.t[end], xlabel="time",
                   title="$(palgo()), taps", titleloc=:left)
    elseif show == "y*w"
        plt = plot([t->(y_i(sol2, i)(t) * w_i(sol2, i)(t)) for i in 1:n],
                   sol2.t,
                   labels=["y_$i" for i in 1:n], leg=:bottomright,
                   xticks=sol2.t[1]:tp:sol2.t[end], xlabel="time",
                   title="$(palgo()), taps")
    elseif show == "zyw"
        plt = plot([t->z(sol2, t), y_i(sol2, 1), w_i(sol2, 1),
                    io.z,          y_i(sol2, n), w_i(sol2, n)],
                   sol2.t,
                   label=["approx" "y_1" "w_1" "target" "y_n" "w_n"],
                   title=["$algo (MSE: $error)" "taps" "weights"],
                   layout=(3,1))
    end

    if return_mse
        return plot, error
    end
    return plt
end

function imp_plot(show, sol)
    """ Returns a plot. """
    if !(show in ["imp", "impi"]) return end
    sol_imp = impulse_response(sol)

    if show == "imp"
        plot(t->z(sol_imp, t),
             sol_imp.t,
             label="total imp", xlabel="t",
             title="$algo, impulse response")
    elseif show == "impi"
        plot([y_i(sol_imp, i) for i in 1:n],
             sol_imp.t,
             labels=["y$i" for i in 1:n], xlabel="t",
             title="$algo, individual impulse response")
    end
end

#=== RUNNING + SAVING ======#

function serialize(u_end, filename, error)
    dict = Dict("algo"    => algo,
                "task"    => task,
                "n"       => n,
                "eta"     => η,
                "alpha"   => α,
                "latency" => latency,
                "mse"     => error,
                "u"       => u_end
                )
    string = JSON.json(dict)
    open(filename, "w") do f
        write(f, string)
    end
end

if exec
    # init
    io, tp, u0, tsave = init_training(task, training_cycles)

    # train
    sol = run(io, u0, tsave)

    # plot
    evol_plot(show, sol, io)
    _, error = end_plot(show, sol, io, tp, return_mse=true)
    imp_plot(show, sol)
    println("MSE: ", error)

    # save figure
    fn = "$algo-$task-$show-$mark-alpha$α-latency$latency-eta$η-offset$offset-t$(tsave[end])"
    dir = task
    savefig("../figures/$dir/$fn.png")

    # serialize the model
    #if mode == "train"
    #    serialize(sol1, "../models/$fn", error)
    #end
end

max_stars_repo_name: toludwig/AdaptiveFilterNeuron

max_issues_repo_name: toludwig/AdaptiveFilterNeuron

max_forks_repo_name: toludwig/AdaptiveFilterNeuron


text:
# Example 1

using Base.Iterators
using Test
using Arbitrary

# Generate arbitrary values
xs = collect(take(arbitrary(BigInt), 100))
ys = collect(take(arbitrary(BigInt), 100))
zs = collect(take(arbitrary(BigInt), 100))

# Test commutativity
@test all(xs .+ ys .== ys .+ xs)
# Test associativity
@test all((xs .+ ys) .+ zs .== xs .+ (ys .+ zs))


text:
#######################################################################
function test_run_bpca(σ, seed, S, maxiter)
#######################################################################

    x, Σ = simulatedatapca(σ, seed)#
    #x, Σ = realresampleddataforpca()#

    posterior = bpca(x', Σ'; maxiter=maxiter, S=S)

    figure(1)

    cla()

    plot3D(x[:,1], x[:,2], x[:,3], "bo", alpha=0.3)

    for i in 1:30

        D = size(x, 2)

        param = rand(posterior)

        W, μ = param[1:D], param[D+1:2*D]

        project(x) = ((W'*W) \ (W'*(x - μ)))[1]

        fwd(z::Real) = vec(W*z + μ)

        z = sort([project(x[i,:]) for i=1:size(x,1)])

        zgrid = collect(LinRange(minimum(z), maximum(z), 100))

        X = reduce(hcat, fwd.(zgrid))'

        plot3D(X[:,1], X[:,2], X[:,3], "-")

    end

    return posterior

end


text:
"""
Creates ordering of length N from equivalent indices.
`eq_inds` can only contain unique values from 1:N.
When used for a spin system, N is the number of nuclei in the system.

Returns the indices for the N nuclei, and the number of degrees of freedom.

Example:
import Destruct
tmp = collect( createorderingfromeqinds(ME[i], A.N_spins_sys[i]) for i = 1:length(A.N_spins_sys) )
κs_β_ordering, κs_β_DOF = Destruct.destruct(tmp)
"""
function createorderingfromeqinds(eq_inds::Vector{Vector{Int}}, N::Int)

    j = 0 # degrees of freedom counter.
    out = zeros(Int, N)

    if isempty(eq_inds)
        return collect(1:N), N
    end

    # check if `eq_inds` only contains unique values from 1:N.
    if all( all( (eq_inds[k][l] in 1:N) for l = 1:length(eq_inds[k])) for k = 1:length(eq_inds) )
        for k = 1:length(eq_inds)
            j += 1

            out[eq_inds[k]] .= j
        end

        # fill the rest.
        for n = 1:N
            if out[n] == 0
                # out[n] is unassigned.
                j += 1
                out[n] = j
            end
        end

        return out, j
    end

    println("Invalid ME, using default.")
    return collect(1:N), N
end

"""
ordering must contain integers from 1:N, where N is the number of unique entries in ordering.

Example:
Δc_m_compound_i_l = randn(10)
κs_β_ordering_i = [ 2; 2; 2; 3; 4; 5;  1; 1; 1; 6]
y = condensenuclei(Δc_m_compound_i_l, κs_β_ordering_i)

"""
function condensenuclei(x::Vector{T}, ordering::Vector{Int}, N::Int)::Vector{T} where T
    @assert length(x) == length(ordering)
    @assert norm(collect(1:N) - sort(unique(ordering))) < 1e-14

    y = zeros(T, N)
    #Ns = zeros(Int, N)
    for i = 1:length(x)

        k = ordering[i]
        y[k] += x[i]
        #Ns[k] += 1
    end

    #return y ./ Ns
    return y
end

function condensenuclei(x::Vector{T}, ordering::Vector{Int})::Vector{T} where T
    @assert length(x) == length(ordering)

    N = length(unique(ordering))

    return condensenuclei(x, ordering, N)
end

"""
unused?
returns a reduced length (in the most nested level) version of Δc_m_compound.
ME is equivalence indices. Could be an entry of MEs, the return type of getmageqinfomixture().
Δc_m_compound is a field in the SHType composite type.
"""
function reduceΔc(Δc_m_compound::Vector{Vector{Vector{T}}},
    ME::Vector{Vector{Vector{Int}}},
    N_spins_sys::Vector{Int}) where T

    N_sys = length(Δc_m_compound)
    @assert length(ME) == N_sys == length(N_spins_sys)

    # prepare.
    κs_β_ordering = Vector{Vector{Int}}(undef, N_sys)
    κs_β_DOF = Vector{Int}(undef, N_sys)
    for i = 1:N_sys
        κs_β_ordering[i], κs_β_DOF[i] = createorderingfromeqinds(ME[i], N_spins_sys[i])
    end

    # condense.
    c = Vector{Vector{Vector{T}}}(undef, N_sys)

    for i = 1:N_sys # over spin systems.

        if !isempty(Δc_m_compound[i])

            c[i] = Vector{Vector{T}}(undef, length(Δc_m_compound[i]))

            for l = 1:length(Δc_m_compound[i]) # over existing groups in Δc_m_compound.

                for k = 1:length(ME[i])
                    c[i][l] = condensenuclei(Δc_m_compound[i][l], κs_β_ordering[i], κs_β_DOF[i])
                end
            end
        else
            c[i] = Δc_m_compound[i]
        end
    end

    return c
end

"""
returns a reduced length (in the most nested level) version of Δc_m_compound.
ME_m is equivalence indices for a spin system. Could be a spin group from an entry of MEs, the return type of getmageqinfomixture().
Δc_m is a field in a spin group of the SHType composite type.

Example: (work on later, see reduce/batch.jl)
c = NMRSpectraSimulator.reduceΔc(A.Δc_m_compound, ME, A.N_spins_sys)
c2 = NMRSpectraSimulator.reduceΔc(A.Δc_m_compound[1], ME[1], A.N_spins_sys[1])

"""
function reduceΔc(Δc_m::Vector{Vector{T}},
    ME_m::Vector{Vector{Int}},
    N_spins::Int) where T

    if isempty(ME_m)
        return Δc_m
    end

    # prepare.
    ordering, DOF = createorderingfromeqinds(ME_m, N_spins)

    # condense.
    c = Vector{Vector{T}}(undef, length(Δc_m))
    for l = 1:length(Δc_m) # over existing groups in Δc_m_compound.

        for k = 1:length(ME_m)
            c[l] = condensenuclei(Δc_m[l], ordering, DOF)
        end
    end

    return c
end


text:
ENV["GUROBI_HOME"] = "/home/simon/software/gurobi912/linux64/"
ENV["GRB_LICENSE_FILE"] = "/home/simon/software/gurobi912/gurobi.lic"

using Gurobi
using JuMP
using LinearAlgebra
using MeshCat
using Plots
using Test

include("structures.jl")
include("visuals.jl")

# probem data
T = 40 # horizon (T+1 states, T controls)
H = 380 # simulation horizon
dt = 0.04 # step size
x0 = [0.00, 0.0] # initial state
Q = 1.0 # cost on state
Qf = 50.0 # cost on state
R = 1.0 # cost on control
β = 1e3 # extreme values in MIQP

# dynamics model
model = WallPendulum12(2, 1, 1.0, 1.0, 10.0, 1e4, 0.1) # m l g k d

# dynamics model
A1, B1, c1 = dynamics_model(model, dt, mode = :none)
A2, B2, c2 = dynamics_model(model, dt, mode = :left)
A3, B3, c3 = dynamics_model(model, dt, mode = :right)

A = [A1, A2, A3]
B = [B1, B2, B3]
c = [c1, c2, c3]

# domains
C1 = domain(model, mode = :none)
C2 = domain(model, mode = :left)
C3 = domain(model, mode = :right)
C = [C1, C2, C3]

# problem
prob = WallProblem16(model, T, x0, Q, Qf, R, β, A, B, c, C)

# disturbances
i1 = 20
i2 = i1 + 100
i3 = i2 + 40
i4 = i3 + 100
i5 = i4 + 60
ind = [i1, i2, i3, i4, i5]
w = [[-15.5], [+15.5], [+15.5], [-11.5], [-10.5]]

# simulation
x, u, s = simulate!(prob, x0, H, w = w, iw = ind)
sum(s) / (H * dt) # 4.52
maximum(s) / dt #17.6
mean(s) # 0.18
sqrt(mean((s .- mean(s)).^2)) #

# visualization
DTx = range(0.0, step = dt, length = size(x)[1])
DTu = range(0.0, step = dt, length = size(u)[1])
plot(DTu, hcat(u...)', label = "u")
plot(DTx[1:end], hcat(x...)'[1:end,:], label = ["θ" "θd"])

# vis = Visualizer()
# open(vis)
build_robot!(vis, model)
set_robot!(vis, model, [0.10, 0.1])
anim = visualize_robot!(vis, model, x, h = dt)

pθ_right = generate_pusher_traj(ind, w, x, side=:right)
pθ_left  = generate_pusher_traj(ind, w, x, side=:left)
visualize_disturbance!(vis, model, pθ_right, anim=anim, sample=1, offset=0.02, name=:PusherRight)
visualize_disturbance!(vis, model, pθ_left,  anim=anim, sample=1, offset=0.02, name=:PusherLeft)

# filename = "wall_pendulum_with_pusher"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/video/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/video/$filename.mp4",
#     "/home/simon/Documents/video/$filename.gif", overwrite=true)

# system simplified compared to PushBot
    # only 1 control and 2 states
    # linearized about equilibrium point
    # only 3 contact modes considered (no stick-slip modes, no friction considered)
# performance compared to PushBot with LCI-MPC
    # on average 5 time slower than real-time on the same computer
    # especially struggling with contact switches (steps are 15x slower than real-time MPC)


text:
export IFeature

abstract type IFeature end

function layers(this::IFeature, available::Vector{String})::Vector{String}
    return []
end

function instanceExtensions(this::IFeature, available::Vector{String})::Vector{String}
    return []
end

function deviceExtensions(this::IFeature)::Vector{String}
    return []
end

function onInstanceCreated(this::IFeature, vkInstance::vk.VkInstance)
end

function beforeInstanceDestruction(this::IFeature)
end

function onLogicalDeviceCreated(this::IFeature, device)
end

function onPhysicalDeviceSelected(this::IFeature, phy::vk.VkPhysicalDevice)
end

function beforeDeviceDestruction(this::IFeature)
end

function supportsDevice(this::IFeature, dev::vk.VkPhysicalDevice)::Bool
    return true
end

function queueRequests(this::IFeature, families::Vector{vk.VkQueueFamilyProperties})
    return []
end

function addPhysicalDeviceFeatures(this::IFeature, outDeviceFeatures::VkExt.VkPhysicalDeviceFeatures)
end


text:
module TestGBIFRecords

  using GBIF
  using Test

  qpars = Dict("scientificName" => "Mus musculus", "year" => 1999, "hasCoordinate" => true)
  set = occurrences(qpars)
  @test typeof(set) == GBIFRecords
  @test length(set) == 20

  set2 = occurrences()
  @test typeof(set2) == GBIFRecords
  @test length(set2) == 20

end


text:
using Sockets
using Printf

const file_to_lint = @sprintf "%s\n" Base.ARGS[1]
const TIMEOUT = 300 # Number of seconds to wait before closing

"""
    read_timeout(tcp)

Read line over TCP with a timeout. Copied from
`https://github.com/JuliaLang/julia/issues/36217#issuecomment-823344326`.
"""
function read_timeout(tcp::TCPSocket)
    # create timeout timer
    status_before = tcp.status
    tm = Timer(TIMEOUT) do
        lock(tcp.cond)
        tcp.status = Base.StatusEOF
        notify(tcp.cond)
        unlock(tcp.cond)
    end
    # read synchronously
    out = readline(tcp)
    # check timer
    if isopen(tm) # readline worked
        close(tm)
        return out
    else # got EOF because of timer
        # Read timed out!
        tcp.status = status_before
        return nothing
    end
end

# Set to server
try
    # Send a file path to server for linting
    conn = connect(1111)
    write(conn, file_to_lint)
    # Receive error messages back
    while true
        response = read_timeout(conn)
        if isnothing(response) || response == "<<end_lint>>"
            break
        end
        @printf "%s\n" response
    end
catch err
    # No server
end


text:
# interface to C++ library "figtree" from
#   Vlad I. Morariu, Balaji Vasan Srinivasan, Vikas C. Raykar,
#   Ramani Duraiswami, and Larry S. Davis. Automatic online tuning for
#   fast Gaussian summation. Advances in Neural Information Processing
#   Systems (NIPS), 2008.
# and website that hosts code is
#   http://www.umiacs.umd.edu/~morariu/figtree/

using LinearAlgebra, Plots
gr(size=(900,300))

#include("../../src/part/fgt.jl")  # loads in functions to be tests

function fgtTest()
    d = 2                           # dimension of data
    M = 1_000                      # number of targets
    N = 2
    #N = 1_000                      # number of sources
    h = 0.8                         # bandwidth
    ε = 1e-2                        # max tolerance
    x = [-0.5,0.5]
    #x = rand(N*d)                   # source vectors, reshape'd to 1D array
    y = rand(M*d)                   # target vectors, reshape'd to 1D array
    W = 1                           # number of sources being evaluated
    q=[0.5,1.0]
    #q = rand(N*W)                   # weights for sources
    v1 = rand(M)                     # output values array
    v2 = rand(M)                    # output 2
    #v3 = rand(M)                    # output 3

    temp = norm(v1-v2)

    fgt!(v1,d,M,N,h,ε,x,y,q,W)
    dgt!(v1,d,M,N,h,ε,x,y,q,W)
    #mydgt!(v3,d,M,N,h,ε,x,y,q,W)

    println("fast (C++) transform: ")
    @time fgt!(v1,d,M,N,h,ε,x,y,q,W)
    println("slow (C++) transform: ")
    @time dgt!(v2,d,M,N,h,ε,x,y,q,W)
    #println("slow (Julia) transform: ")
    #@time mydgt!(v3,d,M,N,h,ε,x,y,q,W)
    println()

    println("(check that two algorithms give same result)")
    println("initial norm: ",temp)
    println("later norm:   ",norm(v1-v2))
end

function fgt1Dtest()
    d = 1                           # dimension of data
    M = 10_000                      # number of targets
    #N = 2
    N = 10_000                      # number of sources
    h = 0.2                         # bandwidth
    ε = 1e-2                        # max tolerance
    #x = [-0.5,0.5]
    x = rand(N*d)                   # source vectors, reshape'd to 1D array
    y = collect(range(-2,stop=2,length=M))                   # target vectors, reshape'd to 1D array
    W = 1                           # number of sources being evaluated
    #q=[0.5,1.0]
    q = rand(N*W)                   # weights for sources
    v1 = rand(M)                     # output values array
    v2 = rand(M)                    # output 2
    v3 = rand(M)                    # output 3

    @time fgt!(v1,d,M,N,h,ε,x,y,q,W)
    @time dgt!(v2,d,M,N,h,ε,x,y,q,W)
    @time mydgt!(v3,d,M,N,h,ε,x,y,q,W)

    p1 = plot(y,v1,label="fgt (c++)")
    p2 = plot(y,v2,label="dgt (c++)")
    p3 = plot(y,v3,label="mydgt (julia)")

    l = @layout [a a a]
    p = plot(p1,p2,p3,layout = l)
    gui(p)
end