text:
# surrogate_ground_truth_ais-multivariatetpseudohuber
#
# This script runs a long chain based on the
# Metropolis-adjusted Riemannian Langevin algorithm (MARLA). This is used as
# a gold standard as this procedure yields an asymptotically unbiased chain
# and thus is used to compare the efficiency of the biased chains produced
# in the experiment compare-hyperparameters-multivariatetpseudohuber.

using DataFrames

using SteinDistributions: SteinMultivariateStudentTRegressionPseudoHuberPrior
using SteinSamplers: runmarla

include("experiment_utils.jl")
include("ais/ais_utils.jl")

# Set random number generator seed
@parseintcli seed "s" "seed" 7
# set the delta param for the Huber prior
@parsestringcli clinu "n" "nu" "10.0"
nu = float(clinu)
# set the nu param degrees of freedom for the d-distribution
@parsestringcli clidelta "d" "delta" "1.0"
delta = float(clidelta)

# distname is multivariatetpseudohuber
distname = "multivariatetpseudohuber"
###            Pull in data here             ###
### ais dataset is available in sn R package ###
data_set = "ais"
yfeature = "LBM"
dat = readtable("src/experiments/ais/ais.tsv")
X = convert(Array{Float64,2}, dat[[:RCC, :WCC, :Fe]])
# standardize columns of X
X = X .- mean(X,1)
X = [ones(size(X,1)) X]
(N, d) = size(X)
column_norms = zeros(1, d)
for ii in 1:d
    column_norms[1,ii] = norm(X[:,ii])
end
X = X ./ column_norms
y = convert(Array{Float64,1}, dat[:LBM])
y = (y - mean(y)) / std(y)
# number of samples to make
n = 5000000
# the thinning rate
thinby = 20
# create the target distribution
sigma = Diagonal(ones(size(X,1)))
target = SteinMultivariateStudentTRegressionPseudoHuberPrior(X, y, nu, sigma, zeros(d), delta)
# starting beta0
beta0 = getapproxmode(target)
# the epsilons to run the experiemnt
epsilons = [1e-2, 1e-3, 1e-4, 1e-5]
# define the covariance terms for the Riemannian-Langevin diffusion
function volatility_covariance(x::Array{Float64,1})
    p = length(x)
    sqrt(1 + sum(x.^2) / delta^2) * eye(p)
end
function grad_volatility_covariance(x::Array{Float64,1})
    delta^(-2) * x ./ sqrt(1 + sum(x.^2) / delta^2)
end

for epsilon = epsilons
    # Draw MARLA samples
    (betas_marla, numgrad_marla, acceptance_ratio) =
        @setseed runmarla(target,
                          beta0,
                          volatility_covariance,
                          grad_volatility_covariance;
                          epsilonfunc=t -> epsilon,
                          verbose=true,
                          numiter=n*thinby)
    # Thin the samples
    betas_marla_thinned = betas_marla[thinby:thinby:(thinby*n),:]

    # Package and save results
    summary_data = Dict{Any, Any}(
        "n" => n,
        "seed" => seed,
        "distname" => distname,
        "epsilon" => epsilon,
        "nu" => nu,
        "delta" => delta,
        "betas" => betas_marla_thinned,
        "numgrad" => numgrad_marla,
        "acceptance_ratio" => acceptance_ratio,
        "sampler" => "ais-marla"
    )
    save_json(
        summary_data;
        dir="compare-hyperparameters-multivariatetpseudohuber-MARLA",
        distname=distname,
        dataset=data_set,
        n=n,
        thinby=thinby,
        epsilon=epsilon,
        nu=nu,
        delta=delta,
        d=d,
        seed=seed
    )
end


text:
############
# Priors
############
#
# need logprior(theta) = p(.)*p(.)
#
# assume indpendence for now

# priors = Dict{Symbol, Union{Function, Vector{Function}}()

"""
$(SIGNATURES)

Combines the independent component-wise specified priors in
`component_wise_priors` with any free-form priors in
`full_theta_priors`.

`component_wise_priors` is a `Dict{Symbol,Any}` which maps the name of
the parameter to a probability function `p`.  This `p` takes the value
of the parameter and returns the log-prior.

`full_theta_priors` is a tuple of functions.  Each function takes a
`theta` vector and returns the log-prior for it.

If only theta is passed in then it uses `def_component_wise_priors()`
and `def_full_theta_priors()` as defaults.

Return:
- the function `logprior(theta)`
- a function `logprior_debug(theta)` which will print the values of the individual
  components of the prior and return them
"""
function make_logprior(theta0::Theta0,
                       component_wise_priors::Dict{Symbol}=def_component_wise_priors(),
                       full_theta_priors::Dict{Symbol}=def_full_theta_priors(theta0))
    cwps = Dict{Symbol,Any}()
    for (i,name) in enumerate(keys(theta0.names))
        if haskey(component_wise_priors, name)
            cwp = component_wise_priors[name]
            inds = theta0.names[name]
            f = _make_one_logprior(inds, cwp)
            cwps[name] = f
        end
    end
    len_theta0 = length(theta0.th0)

    # return function which adds all prior-functions
    logprior = function (theta::Vector)
        # make sure the length of theta is correct
        @assert length(theta)==len_theta0 "$(length(theta)) != $len_theta0"
        out = 0.0
        for fn in values(cwps)
            out += (fn(theta))::Float64
        end
        for fn in values(full_theta_priors)
            out += (fn(theta))::Float64
        end
        out
    end

    logprior_debug = function (theta::Vector, verbose=true)
        # make sure the length of theta is correct
        @assert length(theta)==len_theta0 "$(length(theta)) != $len_theta0"

        out = 0.0
        out2 = Float64[]
        verbose && println("** Component-wise priors:")
        for (k,fn) in cwps
            out += (fn(theta))::Float64
            push!(out2, fn(theta))
            verbose && println("$k: $(fn(theta))")
        end
        verbose && println("** Full priors:")
        for (k,fn) in full_theta_priors
            out += (fn(theta))::Float64
            push!(out2, fn(theta))
            verbose && println("$k: $(fn(theta))")
        end
        out, out2
    end

    return logprior, logprior_debug
end

function _make_one_logprior(theta_inds, prior::Function)
    function (theta::Vector)
        out = 0.0
        for i in theta_inds
            out += prior(theta[i])
        end
        out
    end
end

## Useful prior functions

"""
    logprior_scalar(val, expectation, sigma, range_abs, range_rel, range_rel_cutoff)

For a given value `val` and its `expectation`, return a
log-probability according to the sigma, range_abs, range_rel,
range_rel_cutoff (the fields of SomeData).  It assumes a truncated
Normal distribution.
"""
function logprior_scalar(val, expectation, sigma, range_abs, range_rel, range_rel_cutoff)
    # check that val is in abs-range and rel-range:
    range_abs[1]<=val<=range_abs[2] || return -Inf
    if abs(expectation)>range_rel_cutoff
        range_rel[1]<=val/expectation<=range_rel[2] || return -Inf
    end

    # Gaussian
    return -(val-expectation)^2/(2*sigma^2) - log(sigma)
end

"""
    logprior_scalar_subnormal(val, expectation, sigma, range_abs, range_rel, range_rel_cutoff)

For a given value `val` and its `expectation`, return a
log-probability according to the sigma, range_abs, range_rel,
range_rel_cutoff (the fields of SomeData).  If val is within the range,
then assume a Normal distribution, outside a super-Normal with power 10.

I.e. once outside the range the probability falls sharply.
"""
function logprior_scalar_subnormal(val, expectation, sigma, range_abs, range_rel=(-Inf,Inf), range_rel_cutoff=Inf)
    # Normal dist
    out = -(val-expectation)^2/(2*sigma^2) - log(sigma)

    # if val is outside abs-range or rel-range use sub-Normal, aka sub-Gaussian
    exponent = 10
    if val<range_abs[1]
        out += -(val-range_abs[1])^exponent
    elseif range_abs[2]<val
        out += -(val-range_abs[2])^exponent
    end

    @assert  range_rel_cutoff==Inf # not implemented

    return out
end

"""
    logprior_1dfields(theta, theta0::Theta0, name)

Priors for the parameters encoding 1d fields (aside-reminder: those
parameters are offsets).  Note that we do the priors on the parameters
and not the field in order to avoid correlations.

The values for the prior distribution come from the `SomeData`
structure of the gl.bdot, gl.dhdt, gl.fsl, gl.temp fields.

Somewhat inconsistently:
- the SomeData-ranges are checked for the value at all elevation bands
- the theta-values themselves are then assumed to have a Gaussian PDF
  with SomeData.sigma and mu=0 (as they are offsets).

This is a full_theta_priors.
"""
function logprior_1dfield(theta, theta0::Theta0, name)
    gb = theta0.gb
    gl = gb.gl
    field1d = make_fields1d(theta, theta0, name)
    field1d_expectation = getfield1d(gb, name)
    if name==:btilde # combo of bdot and dhdt
        # here we need to get
        sigma = sqrt(gl.bdot.sigma^2 +
                     gl.dhdt.sigma^2)
        @assert all(gl.bdot.range_rel==gl.dhdt.range_rel) "range_rel need to be equal for bdot and dhdt"
        @assert all(gl.bdot.range_rel_cutoff==gl.dhdt.range_rel_cutoff) "range_rel_cutoff need to be equal for bdot and dhdt"
        @unpack range_rel, range_rel_cutoff = gl.bdot

        range_abs = (gl.bdot.range_abs[1]+gl.dhdt.range_abs[1],
                     gl.bdot.range_abs[2]+gl.dhdt.range_abs[2])
    else
        @unpack sigma, range_abs, range_rel, range_rel_cutoff = getfield(gl, name)
    end

    # check that all are in abs-range and rel-range:
    @inbounds for (i,fd) in enumerate(field1d)
        range_abs[1]<=fd<=range_abs[2] || return -Inf
        if abs(field1d_expectation[i])>range_rel_cutoff
            range_rel[1]<=fd/field1d_expectation[i]<=range_rel[2] || return -Inf
        end
    end

    # Gaussian distribution: run only over actual parameters as their
    # value on the elevation bands will be correlated.
    inds = theta0.names[name]
    out = 0.0
    @inbounds @fastmath for i in inds
        thi = theta[i]
        out -= thi*thi
    end
    return out/(2*sigma^2)
end

"""
    calc_expected_terminus_flux(terminus_flux::Tuple, gb::Bands)
    calc_expected_terminus_flux(fkind::FluxKind._FluxKind, vals, gb::Bands)

Return expected flux at terminus in m^3/s.  NOTE, this value will still be scaled in
terminus_flux_logprior
"""
calc_expected_terminus_flux(terminus_flux::Tuple, gb::Bands) =
    calc_expected_terminus_flux(terminus_flux[1], terminus_flux[2].vals, gb::Bands)
function calc_expected_terminus_flux(fkind::FluxKind._FluxKind, vals, gb::Bands)
    tflux_expectation = if fkind==FluxKind.abs
        # absolute flux
        return vals
    elseif fkind==FluxKind.rel_bdot_teminus
        flux_n(gb, gb.bdot)[end] # flux at terminus due to expectation bdot
    elseif fkind==FluxKind.rel_bdot_max
        out = maximum(flux_n(gb, gb.bdot)) # maximum flux due to expectation bdot (== total accumulation rate)
        if out<=0
            # if all bdot is negative, i.e. no accumulation area, use the maximum apparent mass balance:
            maximum(flux_n(gb, gb.bdot-gb.dhdt))
        else
            out
        end
    elseif fkind==FluxKind.rel_btilde_teminus
        (flux_n(gb, gb.bdot-gb.dhdt)[end]) # expected flux at terminus
    else
        error("Unknown flux-kind: $fkind")
    end
    # this needs to be bigger than zero because sigma is relative... hack it!
    #@assert tflux_expectation>0 "Expectation value of terminus flux zero or negative: $tflux_expectation"
    if !(tflux_expectation>0)
        tflux_expectation = -1.0
    end
    return tflux_expectation
end

function _calc_flux(theta, theta0::Theta0)
    gb = theta0.gb
    # calculate 1d fields
    btilde = make_fields1d(theta, theta0, :btilde) # value at current MCMC step
    # flux
    return flux_n(gb, btilde)
end

"""
    terminus_flux_logprior(theta::Vector, theta0::Theta0)

A prior for the whole of btilde-1D to enforce that the total flux at
the terminus is according to specs.  See `Glacier` type-def.

Notes:
- this is a full_theta_prior
"""
function terminus_flux_logprior(theta::Vector, theta0::Theta0)
    gb = theta0.gb
    gl = gb.gl
    fkind, terminus_flux = gl.terminus_flux
    fkind==FluxKind.none && return 0.0
    @unpack vals, sigma, range_abs, range_rel, range_rel_cutoff = terminus_flux
    tflux = _calc_flux(theta, theta0)[end]
    # expected flux at terminus
    tflux_expectation = calc_expected_terminus_flux(fkind, vals, gb)

    # hack to catch glaciers with no accumulation area
    # TODO make this proper
    if fkind!=FluxKind.abs && tflux_expectation==-1
        tflux_expectation = 0.0
        # some heuristic:
        heuristic_max_flux = area(gl)/1e6/100
        sigma = heuristic_max_flux / 10
        fkind = FluxKind.abs
    end

    if fkind==FluxKind.abs
        return logprior_scalar_subnormal(tflux, tflux_expectation, sigma, range_abs, range_rel, range_rel_cutoff)
    elseif fkind==FluxKind.rel_bdot_teminus || fkind==FluxKind.rel_btilde_teminus || fkind==FluxKind.rel_bdot_max
        # all SomeData-stuff is relative here...
        @assert range_rel == (-Inf,Inf)
        @assert range_rel_cutoff == Inf
        # convert sigma etc to relative values
        return logprior_scalar_subnormal(tflux, tflux_expectation*vals, tflux_expectation.*sigma,
                                         tflux_expectation.*range_abs, range_rel, range_rel_cutoff)
    else
        error("FluxKind $fkind not implemented")
    end
end

"""
    continuous_glacier_logprior(theta::Vector, theta0::Theta0)

Prior to make sure that the glacier is not in two or more pieces,
i.e. that h=0 for some stretch and the h>0 again downstream.
Only takes 0 or -Inf as values.

Note: this will impact the allowed values of btilde only.
"""
function continuous_glacier_logprior(theta::Vector, theta0::Theta0)
    flux = _calc_flux(theta, theta0)
    @assert length(flux)>0
    ifirst = findfirst(flux.>0)
    ilast = findlast(flux.>0)
    if ifirst==0 || ilast==0
        # there is some odd bug I try to catch...
        println("\n========================================\n=== continuous_glacier_logprior bug ===")
        @show flux
        @show theta
        @show theta0
        println("\n========================================\n")
        return -Inf
    end
    return any(flux[ifirst:ilast].<=0) ? -Inf : 0.0
end

"""
    continuous_glacier_from_top_logprior(theta::Vector, theta0::Theta0)

Prior to make sure that the glacier starts at the top and is not in
two or more pieces, i.e. that h=0 for some stretch and the h>0 again
downstream.  Only takes 0 or -Inf as values.

Note: this will impact the allowed values of btilde only.
"""
function continuous_glacier_from_top_logprior(theta::Vector, theta0::Theta0)
    flux = _calc_flux(theta, theta0)
    @assert length(flux)>0
    ilast = findlast(flux.>0)
    ilast==0 && return -Inf
    @assert flux[1]==0
    return any(flux[2:ilast].<=0) ? -Inf : 0.0
end

"""
    def_component_wise_priors()

Default component-wise priors.  Probably mostly for MPara.
"""
function def_component_wise_priors()
    Dict(:dist_exp => x-> 0.1<x<=1 ? 0.0 : -Inf,
         :alpha_min => x -> 0<x ? 0.0 : -Inf,
         :iv_h_exp => x -> 0.1<x<5 ? 0.0 : -Inf,
         :iv_dist_exp1 => x -> 0<=x<=1 ? 0.0 : -Inf, # Huss-type
         :iv_dist_exp2 => x -> 0<x<20 ? -(x-4)^2/(2*0.5^2) : -Inf, # Nye-type; theory says n+1, effects for >20 are minimal, thus cap it.
         :sigma_h_model => x -> 0<x<200 ? 0.0 : -Inf, # -(x-20)^2/(2*10^2)
         :sigma_iv_model => x -> 0<x<200 ? 0.0 : -Inf, # -(x-10)^2/(2*10^2)
         )
end

"""
    def_full_theta_priors(theta0)

Default full-theta priors.  Includes all priors on the 1D fields and
 the prior on the flux at the terminus `terminus_flux_logprior`.
"""
function def_full_theta_priors(theta0)
    # extra priors, namely the terminus_flux and that the glacier is continuous:
    out = Dict{Symbol,Function}(
        :teminus_flux => theta -> terminus_flux_logprior(theta, theta0), # uses values of gl.terminus_flux
        :continuous_glacier => theta -> continuous_glacier_from_top_logprior(theta, theta0))
    # priors on the 1d-fields originating from the fields gl.bot, gl.dhdt, gl.temp, gl.fsl
    for n in [:btilde, :temp, :fsl]
        if n in keys(theta0.names)
            out[n] = theta -> logprior_1dfield(theta, theta0, n)
        end
    end
    return out
end

ext: jl


text:
__precompile__(true)

module Discretizers

import Base: extrema
import StatsBase: iqr, moment
import Statistics: std
import SpecialFunctions: lfactorial

export
    # generic types
    AbstractDiscretizer,
    DiscretizationAlgorithm,

    # discretizer types
    CategoricalDiscretizer, # maps 1D continuous values to/from discrete values
    LinearDiscretizer,      # maps discrete values to/from discrete values
    HybridDiscretizer,      # encodes categorical values first and then applies linear discretization

    # discretization algorithms
    DiscretizeUniformWidth,
    DiscretizeUniformCount,
    DiscretizeMODL_Optimal,
    DiscretizeMODL_Greedy,
    DiscretizeMODL_PostGreedy,
    DiscretizeBayesianBlocks,

    # sampling methods
    AbstractSampleMethod,
    SampleUniform,                    # sample uniform from bin bounds
    SampleBinCenter,                  # return bin center
    SAMPLE_UNIFORM,
    SAMPLE_BIN_CENTER,

    # methods
    encode,             # map from natural to discretized state
    decode,             # map or sample from discretized to natural state
    nlabels,            # number of discrete labels
    encoded_type,       # obtain the encoded type, typically discrete
    decoded_type,       # obtain the decoded type, often continuous
    supports_encoding,  # whether discretizer supports encoding given variable
    supports_decoding,  # whether discretizer supports decoding given variable
    binedges,           # compute bin edges using a discretization algorithm
    binwidth,           # the width of the ith bin
    binwidths,          # an array of bin widths

    get_discretization_counts, # get a list of counts; can be used to generate a histogram
    get_histogram_plot_arrays, # returns arr_x and arr_y for plotting a histogram given binedges and counts

    datalineardiscretizer, # build a linear discretizer which maps Inf to a bin

    # linear discretizer
    totalwidth,        # total width of set, max-mina
    bincenters,        # Vector{Float64} list of bin centers
    bincenter          # Float64 center of a particular bin

### source files

include("common.jl")

include("sample_methods.jl")

include("categorical_discretizer.jl")
include("linear_discretizer.jl")
include("hybrid_discretizer.jl")

include("disc_uniformwidth.jl")
include("disc_uniformcount.jl")
include("disc_MODL.jl")
include("disc_bayesianblocks.jl")

end # module


text:
using REPL
using REPL.LineEdit

const orange = "\e[38;5;166m"
const prompt = "MATLAB> "

function repl_eval(script::String, stdout::IO, stderr::IO)
    eval_string(script)
end

function simple_showerror(io::IO, er)
    Base.with_output_color(:red, io) do io
        print(io, "ERROR: ")
        showerror(io, er)
        println(io)
    end
end

function createMatlabRepl(repl, main)
    mat_mode = LineEdit.Prompt(prompt;
    prompt_prefix=orange,
    prompt_suffix=main.prompt_suffix,
    sticky=true);

    hp = main.hist
    hp.mode_mapping[:matlab] = mat_mode
    mat_mode.hist = hp
    mat_mode.on_enter = (s) -> true
    mat_mode.on_done = (s, buf, ok) -> begin
        if !ok
            return REPL.transition(s, :abort)
        end
        script = String(take!(buf))
        if !isempty(strip(script))
            REPL.reset(repl)
            try
                repl_eval(script, repl.t.out_stream, repl.t.err_stream)
            catch y
                # should never reach
                simple_showerror(repl.t.err_stream, y)
            end
        end
        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || REPL.transition(s, main)
    end

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, mat_mode)

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    mk = REPL.mode_keymap(main)

    b = Dict{Any,Any}[
        skeymap, mk, prefix_keymap, LineEdit.history_keymap,
        LineEdit.default_keymap, LineEdit.escape_defaults
    ]
    mat_mode.keymap_dict = LineEdit.keymap(b)

    mat_mode

end

function initializeMatlabRepl(repl = Base.active_repl, initkey = ')')
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    main_mode = mirepl.interface.modes[1]
    mat_mode = createMatlabRepl(mirepl, main_mode)
    push!(mirepl.interface.modes,mat_mode)

    mat_prompt_keymap = Dict{Any,Any}(
        initkey => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, mat_mode) do
                    LineEdit.state(s, mat_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, initkey)
            end
        end
    )

    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, mat_prompt_keymap);
    nothing
end

function matlabReplInitialized(repl)
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    any(:prompt in fieldnames(typeof(m)) && m.prompt == prompt for m in mirepl.interface.modes)
end


text:
using SV, Econometrics, LinearAlgebra, Statistics, DelimitedFiles

function main()
    # these are the true params that generated the data
    σe = exp(-0.736/2.0)
    ρ = 0.9
    σu = 0.363
    θtrue = [σe, ρ, σu] # true param values, on param space
    n = 1000 # sample size
    burnin = 100
    S = 100 # number of simulations
    # load the example data
    y = readdlm("svdata.txt")
    m =  sqrt(n)*aux_stat(y[:]) # generate the statistic
    # or generate some new data, if you prefer
    #m = sqrt(n)*aux_stat(SVmodel([σe, ρ, σu], n, burnin)[1])
    # set up MCMC
    shocks_u = randn(n+burnin,S) # fixed shocks for simulations
    shocks_e = randn(n+burnin,S) # fixed shocks for simulations
    tuning = [0.01, 0.01, 0.01] # fix this somehow
    lb = [0.0, 0.0, 0.0]
    ub = [2.0, 0.99, 1.0]
    θinit = (ub+lb)./2.0
    lnL = θ -> logL(θ, m, n, shocks_u, shocks_e)
    Prior = θ -> prior(θ, lb, ub) # uniform, doesn't matter
    # define things for MCMC
    verbosity = true
    burnin = 100
    ChainLength = 1000
    # initial proposal moves one at a time
    Proposal = θ -> proposal1(θ, tuning, lb, ub)
    chain = mcmc(θinit, ChainLength, burnin, Prior, lnL, Proposal, verbosity)
    # keep every 10th
    i = 1:size(chain,1)
    keep = mod.(i,10.0).==0
    chain = chain[keep,:]
    # now use a MVN random walk proposal
    Σ = cov(chain[:,1:3])
    tuning = 1.0
    for j = 1:2
        P = (cholesky(Σ)).U
        Proposal = θ -> proposal2(θ,tuning*P, lb, ub)
        θinit = vec(mean(chain[:,1:3],dims=1))
        if j == 2
            ChainLength = 10000
        end
        chain = mcmc(θinit, ChainLength, 0, Prior, lnL, Proposal, verbosity)
        accept = mean(chain[:,end])
        #println("tuning: ", tuning, "  acceptance rate: ", accept)
        if accept > 0.35
            tuning *= 1.5
        elseif accept < 0.25
            tuning *= 0.5
        end
        # keep every 10th
        i = 1:size(chain,1)
        keep = mod.(i,10.0).==0
        θinit = vec(mean(chain[:,1:3],dims=1))
        Σ = 0.5*Σ + 0.5*cov(chain[:,1:3])
    end
    # keep every 10th to reduce autocorrelation
    i = 1:size(chain,1)
    j = mod.(i,10.0).==0
    chain = chain[j,:]
    # plain MCMC fit
    posmean = vec(mean(chain[:,1:3],dims=1))
    inci = zeros(3)
    for i = 1:3
        lower = quantile(chain[:,i],0.05)
        upper = quantile(chain[:,i],0.95)
        inci[i] = θtrue[i] >= lower && θtrue[i] <= upper
    end
    p1 = npdensity(chain[:,1]) # example of posterior plot
    p2 = npdensity(chain[:,2]) # example of posterior plot
    p3 = npdensity(chain[:,3]) # example of posterior plot
    display(plot(p1,p2,p3))
    prettyprint([posmean inci])
    return chain
end
main()

max_stars_repo_name: Vkamalvanshi/Econometrics

max_issues_repo_name: Vkamalvanshi/Econometrics

max_forks_repo_name: Vkamalvanshi/Econometrics


text:
# To run:
# include("DI_design.jl")

# DESCRIPTION
# Trims and linearizes the aircraft model at incremental flight speeds, and
# generates the gains of the dynamic inversion controller. Always run after
# modifying the aircraft parameters in the F16_constants.jl file.
#
# ------------------------------------------------------------------------------

# include packages
#using DifferentialEquations
using LinearAlgebra
using SparseArrays
using ControlSystems
using MAT
using Interpolations
using ApproxFun
using Printf
#ENV["MPLBACKEND"]="tkagg"
using Plots
using LaTeXStrings
using JLD
#using PyPlot
#using GR

# include functions
include("atan2.jl")
include("linearize.jl")
include("trimmer.jl")
#include("finp_F16_loop.jl")
include("THRUST.jl")
include("DAMP.jl")
include("DNDR.jl")
include("DNDA.jl")
include("DLDR.jl")
include("DLDA.jl")
include("CN.jl")
include("CL.jl")
include("CM.jl")
include("CZ.jl")
include("CY.jl")
include("CX.jl")
include("RTAU.jl")
include("TGEAR.jl")
include("PDOT.jl")
include("ADC.jl")
include("fix.jl")
include("F16.jl")
# load F-16 constants
include("F16_constants.jl")

# vector of trim speeds [kts]
speed_vec=range(200,stop=800,length=7)
# initialize DI matrices
CATAB = zeros(3,3,length(speed_vec))
CBinvTAB = zeros(3,3,length(speed_vec))
x0_mat = zeros(NSTATES,length(speed_vec))
u0_mat = zeros(NCTRLS,length(speed_vec))
# trim at different speeds
for iv=1:length(speed_vec)
    # set velocity [ft/s]
    speed=speed_vec[iv]*1.688
    # vertical flight path angle [deg]
    gamma=0.0
    # horizontal flight path angle [deg]
    chi=0.0
    # set initial altitude [ft]
    alt=2000.0
    # set initial North-South and East-West position [ft]
    p_north=0.0
    p_east=0.0
    # turn rate [deg/s]
    turnrate=0.0
    # pull-up rate [deg/s]
    pullup=0.0
    # roll rate [deg/s]
    rollrate=0.0
    # initial guess for state and control vectors
    phi_est=turnrate*D2R*speed/32.17
    pndot=speed*cos(gamma*D2R)*cos(chi*D2R)
    pedot=speed*cos(gamma*D2R)*sin(chi*D2R)
    hdot=speed*sin(gamma*D2R)
    x0=[speed; 0.0; 0.0; phi_est; 0; chi*D2R; 0.0; 0.0; 0.0; p_north; p_east; alt;
        0.0]
    u0=[0.2; 0; 0; 0]
    targ_des=[0; 0; 0; rollrate*D2R; pullup*D2R; turnrate*D2R; 0; 0; 0; pndot;
              pedot; hdot; 0; 0]
    # trim aircraft
    x0, u0, iter = trimmer(F16,x0,u0,targ_des)
    # store trim state and controls
    x0_mat[:,iv] = x0
    u0_mat[:,iv] = u0
    # linearized aircraft dynamics
    A, Bmat, C, D = linearize(F16,x0,u0)
    # dynamic inverse matrices
    C_DI=[1 0 0
          0 1 0
          0 0 1]
    CATAB[:,:,iv]=C_DI*A[7:9,7:9]
    CBinvTAB[:,:,iv]=inv(C_DI*Bmat[7:9,[2,3,4]])
end

# command filters time constants
tau_p=1/10
tau_q=1/10
tau_r=1/4
# washout filter for controls
tau_u=5.0;
# disturbance rejection frequency and damping and integrator poles
wnroll_d=10.0
dmproll_d=1.0
wnpitch_d=10.0
dmppitch_d=1.0
wnyaw_d=4.;
dmpyaw_d=1.;
# feedback gains
kp_roll=2*wnroll_d*dmproll_d
ki_roll=wnroll_d^2
kp_pitch=2*wnpitch_d*dmppitch_d
ki_pitch=wnpitch_d^2
kp_yaw=2*dmpyaw_d*wnyaw_d
ki_yaw=wnyaw_d^2;
# pilot gains
k_a=(290*pi/180)/21.5
k_e=(9-1)*G/25 ./(500*1.688)
k_r=1/30
# save DI matrices and gains into a file
save("DI.jld","CATAB",CATAB,"CBinvTAB",CBinvTAB,"kp_roll",kp_roll,"ki_roll",
    ki_roll,"kp_pitch",kp_pitch,"ki_pitch",ki_pitch,"kp_yaw",kp_yaw,"ki_yaw",
    ki_yaw,"k_a",k_a,"k_e",k_e,"k_r",k_r,"tau_p",tau_p,
    "tau_q",tau_q,"tau_r",tau_r,"tau_u",tau_u)

# plot trim pitch attitude and controls
gr(size=(800,600))
p10=plot(speed_vec,x0_mat[5,:]*180/pi,label="",line=(3, :royalblue, :solid), marker = (:circ, :royalblue, 8), markerstrokecolor = :royalblue)
xaxis!((200, 800), 200:50:800)
yaxis!("θ [deg]",yguidefontsize=14)
p11=plot(speed_vec,u0_mat[1,:]*100,label="",line=(3, :royalblue, :solid), marker = (:circ, :royalblue, 8), markerstrokecolor = :royalblue)
xaxis!((200, 800), 200:50:800)
yaxis!(L"\delta_t [\%]",yguidefontsize=14)
p12=plot(speed_vec,u0_mat[3,:],label="",line=(3, :royalblue, :solid), marker = (:circ, :royalblue, 8), markerstrokecolor = :royalblue)
xaxis!((200, 800), 200:50:800)
yaxis!(L"\delta_e [\mathrm{deg}]",yguidefontsize=14)
p13=plot(p10, p12, p11, layout=(3,1), leftmargin=3Plots.mm)
display(p13)
savefig(".\\Plots\\trim_F16.svg")


text:
using TensorBoardLogger
using TraceGraph

function foo(a)
    b = a+10
    c = b+10
    return (a,b,c)
end

tg = tracegraph(foo, 100)
lgr = TBLogger("tblog", tb_overwrite)
log_graph(lgr, "demo", tg.g, step = 1, nodelabel = tg.nodelabel, nodeop = tg.nodeop, nodevalue = tg.nodevalue)


text:
using DrWatson
@quickactivate "ExercisesRepo"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

# %% permutation entropy
using Random, Combinatorics, Statistics
Random.seed!(356)
o = 3
N = 8
x = rand(N)

fig = figure(figsize = (figx, figx/3.5))
ax1 = subplot2grid((1, 3), (0, 0))
ax1.plot(1:N, x, marker = "o", color = "C0", mfc = "C1", ms = 15, zorder = 99)
ax1.set_title("timeseries", size = 28)
ax1.spines["bottom"].set_position("center")
ax1.spines["right"].set_color("none")
ax1.spines["top"].set_color("none")
ax1.set_xticks([])
ax1.set_yticks([])
ax1.set_xlabel("\$t\$")
ax1.set_ylabel("\$x\$", rotation = 0)
ax1.xaxis.set_label_coords(1.0, 0.5)
ax1.yaxis.set_label_coords(-0.05, 0.9)
ax1.set_ylim(-0.1, 1.1)
# ax1.axis("off")

# add pattern indications
for (s, n) in zip(("2", "3", "3"), (2, 4, 6))
    xs = [n, n, n+o-1, n+o-1]
    ym = minimum(x[xs[1]:xs[end]]) - 0.1
    yd = ym - 0.02
    ys = [ym, yd, yd, ym] #.+ 0.05
    ax1.plot(xs, ys, color = "k", lw = 1)
    ax1.text(mean(xs), minimum(ys)-0.02, "#=$s", ha="center", va = "top")
end

# ax1.annotate("4", (2, -0.2), xytext = (4, -0.2), ha="center")

ax2 = subplot2grid((1, 3), (0, 1), colspan = 2)
p = permutations([0, 0.5, 1], o) |> collect
counts = [0, 1, 3, 1, 0, 1]
for (i, a) in enumerate(p)
    ax2.plot((1:o) .+ i*o .+ 1, a, marker = "o",
    color = "C0", mfc = "C2", ms = 10)
    ax2.text(o÷2 + i*o + 1, 1.2, "$i")
    ax2.text(o÷2 + i*o + 1, -0.5, "$(counts[i])")
end
ax2.text(o+1, 1.2, "#", ha = "right")
ax2.text(o+1, 0.5, "pattern", ha = "right")
ax2.text(o+1, -0.5, "count", ha = "right")
ax2.set_ylim(-1, 1.5)
ax2.set_xlim(0, ax2.get_xlim()[2])
ax2.set_title("relative amplitude permutations, \$d=$o\$", size = 28)
ax2.axis("off")

fig.tight_layout()
fig.subplots_adjust(top = 0.9, bottom = 0.02, right = 0.95, left = 0.05, wspace=0.1, hspace = 0.1)

# fsave(joinpath(figdir, "permentropy"), fig)

# %% Koch snowflake
linepoints = [[0.0; 0.0], [1.0; 0.0]]
flakepoints = [[0.0; 0.0], [0.5; sqrt(3)/2], [1; 0.0], [0.0; 0.0]]
function koch(points, maxk, α = sqrt(3)/2)
  Q = [0 -1; 1 0]
  for k = 1:maxk
    n = length(points)
    new_points = Vector{Float64}[]
    for i = 1:n-1
      p1, p2 = points[i], points[i+1]
      v = (p2 - p1) / 3
      q1 = p1 + v
      q2 = p1 + 1.5v + α * Q * v
      q3 = q1 + v
      append!(new_points, [p1, q1, q2, q3])
    end
    push!(new_points, points[end])
    points = new_points
  end
  return points
end

# Plot construction
fig, axes = subplots(2,3, figsize = (2figx/3,figx/2 - 2))
toprow = axes[1:2:5]
botrow = axes[2:2:6]

for (o, ax) in enumerate(toprow)
    ax.clear()
    ax.set_aspect("equal")
    k1 = koch(flakepoints, o-1)
    k2 = koch(flakepoints, o)

    ax.plot([a[1] for a in k1], 0.5 .+ [a[2] for a in k1], COLORS[o])
    ax.plot([a[1] for a in k2], 0.5 .+ [a[2] for a in k2], COLORS[o+1], label="Step $(o)", alpha = 0.75)
    # ax.legend(loc = "upper left", framealpha = 1.0, ncol = 2, handlelength = 1)
    ax.text(-0.1, 0.88, "step $(o)", transform = ax.transAxes, size=18)
    # ax.set_ylim(ys...)
    # ax.set_xlim(xs...)
    ax.set_xticklabels([])
    ax.set_yticklabels([])
    ax.axis("off")
end

# Plot versus circle
points1 = [[0.5; 0.865], [2/3, 0.577361]]
points2 = [[0.537233; 0.801228], [0.555887, 0.768912]]
allpoints = [flakepoints, points1, points2]

for i in 1:3
    botrow[i].clear()
    largekoch = koch(allpoints[i], 6)
    kx, ky = [a[1] for a in largekoch], [a[2] for a in largekoch]
    botrow[i].plot(kx, ky, lw = 0.5)
    botrow[i].set_aspect("equal")
    botrow[i].add_artist(plt.Circle((0.5,0.28900), 0.5, lw=2.0, color="C1", alpha = 0.75, fill=false))
    botrow[i].axis("off")
end
for i in 1:2
    origin = botrow[i]
    zoomin = botrow[i+1]
    zbox = allpoints[i+1]
    axis_zoomin!(zoomin, origin, zbox, zbox, "C$(i+1)", 0.4)
end

fig.tight_layout()
fig.subplots_adjust(wspace = 0.01, hspace = 0.01, bottom = 0.05, left = 0.05)
fsave(joinpath(figdir, "koch"), fig)

# %% Henon fractal zoom
fig, axes = subplots(1,3, figsize = (figx,figx/4))
botrow = axes

he = Systems.henon()
N = 10000
tr = trajectory(he, N, Ttr = 100)
integ = integrator(he)
data1 = Dataset([integ.u])
data2 = Dataset([integ.u])
n, m = 1, 1

zbox1 = ((0.93, 0.05), (1.17, 0.13))
zbox2 = ((1.1275, 0.078), (1.1425, 0.085))
while n < N
    step!(integ)
    if zbox1[1][1] < integ.u[1] < zbox1[2][1] && zbox1[1][2] < integ.u[2] < zbox1[2][2] && m < N
        push!(data1, integ.u)
        global m+=1
    end
    if zbox2[1][1] < integ.u[1] < zbox2[2][1] && zbox2[1][2] < integ.u[2] < zbox2[2][2]
        push!(data2, integ.u)
        global n+=1
    end
end

for ax in botrow[1:2]
    ax.clear()
    ax.plot(tr[:, 1], tr[:, 2], ls = "None", marker = ".", color = COLORS[1], ms = 1)
end
botrow[1].set_xticks([])
botrow[1].set_yticks([])
botrow[2].clear()
botrow[2].plot(data1[:, 1], data1[:, 2], ls = "None", marker = ".", color = COLORS[1], ms = 1)
botrow[2].axis("off")
botrow[3].clear()
botrow[3].plot(data2[:, 1], data2[:, 2], ls = "None", marker = ".", color = COLORS[1], ms = 1)
botrow[3].axis("off")

for (i, zbox) in enumerate((zbox1, zbox2))
    axis_zoomin!(botrow[i+1], botrow[i], zbox, zbox, COLORS[i+1])
    botrow[i+1].set_xlim(zbox[1][1], zbox[2][1])
    botrow[i+1].set_ylim(zbox[1][2], zbox[2][2])
end

fig.tight_layout()
fig.subplots_adjust(top = 0.98, bottom = 0.05, right = 0.95, left = 0.05, wspace=0.1, hspace = 0.1)
# fig.savefig(joinpath(figdir, "henon_zoom.png"))

# %% dimension: Entropy and correlation
he = Systems.henon()
N, α = 10000, 2
X = trajectory(he, N, Ttr = 10)
ε = 10 .^ (-5.0:0.1:1)
H = genentropy(α, ε, X)
C = correlationsum(X, ε)

fig, axs = subplots(2, 1; figsize = (figx/2, figx/2), sharex = true)

axs[1].plot(log.(ε), -H, COLORS[1])

x = -5:0
D = 1.22

axs[2].plot(log.(ε), log.(C), COLORS[1])

for (ax, x) in zip(axs, (-5:0, -10:0))
    Dl = D .* x
    Dl = Dl .- 2
    ax.plot(x, Dl, c = COLORS[3])
    ax.axvspan(x[1], x[end], color = COLORS[2], alpha = 0.25)
    ax.grid("on")
    ax.text(x[2] + 1.5, Dl[2], "\$D_2\$ = $D", color = COLORS[3], size = 30)

end

axs[2].set_xticks(-10:2:2)
axs[1].tick_params(labelbottom=false)
axs[2].set_xlabel("\$\\log(\\varepsilon)\$")
axs[1].set_ylabel("\$-H_$α\$", labelpad = -20)
axs[2].set_ylabel("\$\\log(C)\$", labelpad = -20)
axs[1].set_yticks(-10:2:0)
axs[1].set_yticklabels(["-10", "", "", "", "", "0"])
axs[2].set_yticks(0:-2:-12)
axs[2].set_yticklabels(["0", "", "", "",  "","", "-12"])
fig.tight_layout()
fig.subplots_adjust(left=0.14, bottom = 0.15, hspace = 0.1)
# fsave(joinpath(figdir, "dimension"), fig)

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
# create Vec
@testset "Vec{$ST}" begin
  vtype = PETSc.C.VECMPI
  vec = PETSc.Vec(ST, vtype)
  resize!(vec, 4)
  @test_throws ArgumentError resize!(vec)
  len_ret = length(vec)

  @test length(vec) == 4
  @test size(vec) == (4,)
  @test lengthlocal(vec) == 4
  @test sizelocal(vec) == (4,)
  @test PETSc.gettype(vec) == PETSc.C.VECMPI

  vt = complex(2.,2)  # use vt to hold temporary values
  vec[1] = RC(vt)
  val_ret = vec[1]
  @test vec[1] == RC(vt)

  vec2 = similar(vec,ST)
  PETSc.AssemblyBegin(vec2)
  PETSc.AssemblyEnd(vec2)

  @test isassembled(vec2)
  val2_ret = vec2[1]

  @test val2_ret != val_ret

  if gettype(vec2) == PETSc.C.VECSEQ
    lv2 = localpart(vec2)
    @test lv2 == vec2
  end

  vec_tmp = Vec([1., 2, 3])
  @test PETSc.isfinalized(vec_tmp) == false
  PETSc.PetscDestroy(vec_tmp)
  @test PETSc.isfinalized(vec_tmp) == true

  vec3 = similar(vec, ST, 5)
  @test length(vec3) == 5

  vec4 = copy(vec)
  @test vec4 ≈ vec

  idx = [1,3, 4]
  vt = RC(complex(2.,2))
  vec4[idx] = vt
  vals_ret = vec4[idx]
  @test vals_ret == fill(vt,length(idx))

  vt = RC(complex(3.,3))
  fill!(vec4, vt)

  @test vec4 ≈ fill(vt,length(vec4))

  vt = RC(complex( 4.,4))
  vec4[1:2] = vt

  @test vec4[1:2] == [vt, vt]

  vals = [RC(complex(1,1.)), RC(complex(3.,3)), RC(complex(4., 3))]
  vec4[idx] = vals

  @test vec4[idx] == vals

  vec5 = Vec(Float64, 4)
  varr = LocalVector(vec5)
  @test length(vec5) == 4
  @test length(varr) == length(vec5)
  @test stride(varr, 1) == 1
  vec5j = [1., 2, 3, 4]
  for i=1:length(vec5)  varr[i] = vec5j[i] end

  @test varr[1] == vec5j[1]
  @test varr == vec5j

  varr2 = similar(varr)
  T2 = eltype(varr)
  @test typeof(varr2) == Array{eltype(T2), 1}
  ptr = Base.unsafe_convert(Ptr{T2}, varr)
  @test ptr == varr.ref[]

  restore(varr)

  @test vec5 == vec5j

  varr = LocalVector_readonly(vec5)
  for i=1:length(vec5) @test varr[i] ==  vec5[i] end
  restore(varr)

  # test mlocal constructor
  vec5 = Vec(ST, mlocal=3)
  @test length(vec5) == 3

  @testset "testing logical indexing" begin
      logicals = Array(Bool, length(vec4))
      for i=eachindex(logicals)
        logicals[i] = false
      end
      logicals[2] = true

      vt = RC(complex(5,5.))
      vec4[logicals] = vt

      @test vec4[2] ≈ vt
      @test vec4[1] != vt

      vt = RC(complex(rand(), rand()))
      vals = [vt]
      vec4[logicals] = vals
      @test vec4[2] ≈ vals[1]
      @test vec4[1] != vals[1]

      # reset vec4
      vec4_j = zeros(ST, length(vec4))
      for i=1:length(vec4)
        vec4[i] = RC(complex(Float64(-i), Float64(-i)))
        vec4_j[i] = RC(complex(Float64(-i), Float64(-i)))
      end
  end
  @testset "testing math functions" begin

    @testset "testin chop" begin
       jvec = RC([complex(1.0, 1.0), complex(2.0, 2.0), complex(3.0, 3.0)])
       pvec = Vec(jvec)
       chop!(pvec, RT(1.5))
       jvec[1] = 0.0
       @test pvec ≈ jvec
    end

    vec4_j = zeros(ST, length(vec4))
    for i=1:length(vec4)
      vec4[i] = RC(complex(Float64(-i), Float64(-i)))
      vec4_j[i] = RC(complex(Float64(-i), Float64(-i)))
    end
    @testset "testing abs" begin
      vec4_j = abs(vec4_j)
      absv4  = abs(vec4)
      abs!(vec4)
      if VERSION >= v"0.5.0-dev+0"
          @test real(vec4) ≈ vec4_j
          @test real(absv4) ≈ vec4_j
          @test imag(vec4) ≈ zeros(vec4_j)
          @test imag(absv4) ≈ zeros(vec4_j)
      else
          @test vec4 == vec4_j
          @test absv4 == vec4_j
      end
    end
    @testset "testing exp" begin
      vec4_j = exp(vec4_j)
      exp!(vec4)
      @test vec4 ≈ vec4_j
    end
    @testset "testing log" begin
      vec4_j = log(vec4_j)
      log!(vec4)
      @test vec4 ≈ vec4_j
    end
    onevec = PETSc.Vec(ST, vtype)
    resize!(onevec, 4)
    PETSc.AssemblyBegin(onevec)
    PETSc.AssemblyEnd(onevec)
    for i=1:length(onevec)
        onevec[i] = one(ST)
    end

    @testset "testing norm" begin
      @test_throws ArgumentError norm(onevec,3)
      @test norm(onevec,Inf) == 1
      normvec = copy(onevec)
      PETSc.normalize!(normvec)
      @test norm(normvec,2) == one(ST)
    end
    if ST <: Real
      @testset "testing max and min" begin
        maxvec = copy(onevec)
        maxvec[1] = ST(2)
        @test maximum(maxvec) == 2
        @test findmax(maxvec) == (2.0,1)
        minvec = copy(onevec)
        minvec[1] = ST(0)
        @test minimum(minvec) == 0
        @test findmin(minvec) == (0.0,1)
      end
    end

    @testset "testing pointwise max, min, /" begin
      div1vec = 2*copy(onevec)
      div2vec = 4*copy(onevec)
      @test max(div1vec,div2vec) == div2vec
      @test min(div1vec,div2vec) == div1vec
      @test div1vec .* div2vec == 8*onevec
      @test div2vec ./ div1vec == div1vec
    end
    @testset "testing scale! and negation" begin
      scalevec = scale!(copy(onevec),2)
      @test scalevec == fill(2,length(onevec))
      minusvec = -onevec
      @test minusvec == -onevec
    end

    @testset "testing sum, +, -, *, and /" begin
      @test sum(onevec) == length(onevec)
      multvec = copy(onevec)
      multvec = multvec * 2 * 3 * 4
      @test multvec == 24*onevec
      multvec = copy(onevec)
      multvec = 2 .* multvec
      @test multvec == 2*onevec
      divvec = copy(onevec)
      divvec = divvec * 2 * 3
      divvec = divvec ./ 2
      @test divvec == 3*onevec
      divvec = 3 .\ divvec
      @test divvec == onevec

      divvec = 2*copy(onevec)
      divvec = 2 ./ divvec
      @test divvec == onevec
      addvec = copy(onevec)
      addvec = addvec + 2
      addvec = addvec - 2
      @test addvec == onevec
      addvec = copy(onevec)
      addvec = 2 - addvec
      addvec = 2 + addvec
      @test addvec == 3*onevec
    end
  end

  @testset "testing dot product" begin
    val = dot(vec4, vec)
    val_j = dot(vec4, vec)
    @test val == val_j
  end
  # make copies of vecs 1 2 4

  @testset "testing level 1 Blas" begin

    vecj = zeros(ST, length(vec))
    vec2j = zeros(ST, length(vec))
    vec4j = zeros(ST, length(vec))

    for i=1:length(vec)
      vecj[i] = vec[i]
      vec2j[i] = vec2[i]
      vec4j[i] = vec4[i]
    end

    @testset "testing axpy" begin
      vt = RC(complex(2.,2))
      axpy!(vt, vec, vec2)
      vec2j = vt*vecj + vec2j
      @test vec2j == vec2

      @testset "testing 4 argument axpy" begin
        axpy!(vt, vec, vec2, vec4)
        vec4j = vt*vecj + vec2j
        @test vec2j == vec2
      end

      @testset "testing aypx" begin
        aypx!(vec, vt, vec2)
        vec2j = vt*vec2j + vec
        @test vec2j == vec2
      end

      vt2 = RC(complex(3.,3))
      vt3 = RC(complex(4.,4))
      @testset "testing axpby" begin
        axpby!(vt, vec, vt2, vec2)
        vec2j = vt*vecj + vt2*vec2j
        @test vec2j == vec2

        axpbypcz!(vt, vec, vt2, vec2, vt3, vec4)
        vec4j = vt*vecj + vt2*vec2j + vt3*vec4j
        @test vec4j == vec4
      end

      vecs = Array(typeof(vec), 2)
      vecs[1] = vec
      vecs[2] = vec2
      alphas = [vt2, vt3]
      axpy!(vec4, alphas, vecs)
      vec4j = vec4j + vt2*vecj + vt3*vec2j
      @test vec4j == vec4
    end
    @testset "testing .*, ./, .^" begin
      vec5 = Vec(ST, 3, vtype=PETSc.C.VECMPI)
      vec6 = similar(vec5)
      vec5j = zeros(ST, 3)
      vec6j = zeros(ST, 3)

      for i=1:3
        i_float = Float64(i)

        vec5[i] = RC(complex(i_float, i_float))
        vec6[i] = RC(complex(i_float+3, i_float+3))
        vec5j[i] = RC(complex(i_float, i_float))
        vec6j[i] = RC(complex(i_float +3, i_float+3))
      end

      vec7 = vec5.*vec6
      vec7j = vec5j.*vec6j
      @test vec7 ≈ vec7j

      vec8 = vec5./vec6
      vec8j = vec5j./vec6j
      @test vec8 ≈ vec8j

      vec9 = vec5.^3
      vec9j = vec5j.^3
      @test vec9 ≈ vec9j

      vec10 = vec5 + vec6
      vec10j = vec5j + vec6j
      @test vec10 ≈ vec10j

      vec11 = vec5 - vec6
      vec11j = vec5j - vec6j
      @test vec11 ≈ vec11j
    end
    @testset "test unconjugated dot product" begin
      x = Vec(ST, 2)
      y = Vec(ST, 2)
      copy!(y, [1, 1])
      if ST <: Complex
          copy!(x, [1, im])
          @test (x'*y)[1] == 1-im
          @test (x.'*y)[1] == 1+im
      else
          copy!(x, [2, 3])
          @test (x'*y)[1] == 5
          @test (x.'*y)[1] == 5
      end
    end
  end
  let x = rand(ST, 7)
    @test Vec(x) == x
  end

  @testset "map" begin
    x = rand(3)
    y = Vec(x)
    map!(sin, x)
    map!(sin, y)
    @test x ≈ y
    x2 = map(sin, x)
    y2 = map(sin, y)
    @test x2 ≈ y2

    function myfunc(a, b)
      return a + b
    end

    x3 = copy(x2)
    y3 = copy(y2)
    map!(myfunc, x3, x2, x)
    map!(myfunc, y3, y2, y)
    @test x3 ≈ y3
  end

  @testset "advanced indexing" begin
    x = zeros(ST, 5)
    y = Vec(ST, 5)
    idxs = Int32[0, 1]
    vals = ST[1, 2]
    set_values!(x, idxs, vals)
    set_values!(y, idxs, vals)
    assemble(x)
    assemble(y)
    for i=1:length(idxs)
      @test x[idxs[i]+1] ≈ vals[i]
      @test y[idxs[i]+1] ≈ vals[i]
    end

    vals = ST[2,3]
    ltog = local_to_global_mapping(y)
    set_local_to_global_mapping(y, ltog)
    set_values_local!(x, idxs, vals)
    set_values_local!(y, idxs, vals)
    assemble(x)
    assemble(y)
    for i=1:length(idxs)
      @test x[idxs[i]+1] ≈ vals[i]
      @test y[idxs[i]+1] ≈ vals[i]
    end

    y2 = Vec(ST, 4, bs=2)
    @test get_blocksize(y2) == 2
    idxs = Int32[0]
    vals = ST[1, 2]
    set_values_blocked!(y2, idxs, vals)
    @test y2[idxs[1]+1] ≈ vals[1]
    @test y2[idxs[1]+2] ≈ vals[2]

    rng = localpart(y2)
    ltog = local_to_global_mapping(y2)
    set_local_to_global_mapping(y2, ltog)
    idx = Int32[1]
    vals = ST[2,3]
    set_values_blocked_local!(y2, idxs, vals)
    @test y2[idxs[1]+1] ≈ vals[1]
    @test y2[idxs[1]+2] ≈ vals[2]

  end

end


text:
#__precompile__()
module S3Dicts

using JSON
using AWSCore
#using AWSSDK.S3
using AWSS3
#using Retry
import HTTP
import BigArrays.BackendBase: AbstractBigArrayBackend, get_info, get_scale_name

global const NEUROGLANCER_CONFIG_FILENAME = "info"
global const CONTENT_TYPE = "binary/octet-stream"
global const GZIP_MAGIC_NUMBER = UInt8[0x1f, 0x8b, 0x08]

if haskey(ENV, "AWS_ACCESS_KEY_ID")
    global AWS_CREDENTIAL = AWSCore.aws_config()
elseif isfile("/secrets/aws-secret.json")
    d = JSON.parsefile("/secrets/aws-secret.json")
    global AWS_CREDENTIAL = AWSCore.aws_config(creds=AWSCredentials(d["AWS_ACCESS_KEY_ID"], d["AWS_SECRET_ACCESS_KEY"]))
else
    @warn("did not find AWS credential! set it in environment variables.")
end

export S3Dict

struct S3Dict <: AbstractBigArrayBackend
    bkt         ::String
    keyPrefix   ::String
end

"""
    S3Dict( dir::String )
construct S3Dict from a directory path of s3
"""
function S3Dict( path::String )
    @assert startswith(path, "s3://")
    path = replace(path, "s3://" => "")
    bkt, keyPrefix = split(path, "/", limit = 2)
    keyPrefix = rstrip(keyPrefix, '/')
    S3Dict(bkt, keyPrefix)
end

function get_info(self::S3Dict)
    #data = S3.get_object(AWS_CREDENTIAL; Bucket=self.bkt,
    #              Key=joinpath(dirname(self.keyPrefix), "info"))
    data = s3_get(AWS_CREDENTIAL, self.bkt, joinpath(dirname(self.keyPrefix), "info"))
    return String(data)
end

function get_scale_name(self::S3Dict)  basename( self.keyPrefix ) end

function Base.show( self::S3Dict )  show( joinpath(self.bkt, self.keyPrefix) ) end

function Base.setindex!(h::S3Dict, v::Array, key::AbstractString)
    #@assert startswith(h.dir, "s3://")
    data = reinterpret(UInt8, v[:]) |> Array
    local contentEncoding::String
    if all(data[1:3].== GZIP_MAGIC_NUMBER)
        contentEncoding = "gzip"
    else
        contentEncoding = ""
    end
    #arguments = Dict("Bucket"   => h.bkt,
    #             "Key"      => joinpath(h.keyPrefix, key),
   #              "Body"     => data,
   #              "Content-Type"  => CONTENT_TYPE,
   #              "Content-Encoding" => contentEncoding)
    #@show arguments
    #resp = S3.put_object(AWS_CREDENTIAL, arguments)
    resp = s3_put(AWS_CREDENTIAL, h.bkt, joinpath(h.keyPrefix, key), data;
                  metadata = Dict("Content-Type" => CONTENT_TYPE,
                                  "Content-Encoding" => contentEncoding))
    nothing
end

function Base.getindex(h::S3Dict, key::AbstractString)
    try
        #data = S3.get_object(AWS_CREDENTIAL;
        #                     Bucket=h.bkt,
        #                     Key=joinpath(h.keyPrefix, key))
        data = s3_get(AWS_CREDENTIAL, h.bkt, joinpath(h.keyPrefix, key))
        return data
    catch err
        @show err
        if isa(err, AWSCore.AWSException) && err.code == "NoSuchKey"
            throw(KeyError("NoSuchKey in AWS S3: $key"))
        elseif isa(err, HTTP.ClosedError)
            display(err.e)
            rethrow()
        else
            rethrow()
        end
    end
    nothing
end

function Base.delete!( h::S3Dict, key::AbstractString)
    #S3.delete_object(AWS_CREDENTIAL; Bucket=h.bkt, Key=joinpath(h.keyPrefix, key))
    s3_delete(AWS_CREDENTIAL, h.bkt, joinpath(h.keyPrefix, key))
end

function Base.keys( h::S3Dict )
    #S3.list_objects_v2(AWS_CREDENTIAL; Bucket=h.bkt, prefix=h.keyPrefix)
    s3_list_objects(AWS_CREDENTIAL, h.bkt, h.keyPrefix)
end

function Base.values(h::S3Dict)
    error("normally values are too large to get them all to RAM")
end

function Base.haskey(h::S3Dict, key::String)
    #resp = S3.list_objects_v2(AWS_CREDENTIAL; Bucket=h.bkt, prefix=joinpath(h.keyPrefix, key))
    #return Meta.parse( resp["KeyCount"] ) > 0
    s3_exists(AWS_CREDENTIAL, h.bkt, joinpath(h.keyPrefix, key))
end

end # end of module S3Dicts


text:
using SFRDSP
using Test
using DSP

@testset "SFRDSP.jl" begin
    @testset "stft" begin
        @testset "basic round-trip (L=$L)" for L in 10:17
            x = randn(L)
            X = stft2(x, 8, 8) # window centers at 1, 9, 17
            @test size(X) == (5,3)
            x2 = istft2(X, 8, 8)
            @test length(x2) == 17
            @test x2 ≈ [x; zeros(17-L)]
        end

        # TODO: overlap-add still doesn't work quite right with non-1/2 hop
        # sizes (need to handle the endpoints better)
        @testset "overlap-add (windows=($awin, $swin))" for (awin, swin, scale) in
            [(rect, rect, 2),
             (cosine, cosine, 1),
             (hanning, rect, 1),
             (rect, hanning, 1)]
            x = randn(32)
            X = stft2(x, 8, 4, window=awin)
            x2 = istft2(X, 8, 4, window=swin)
            @test x2 ≈ scale * [x; zeros(length(x2)-length(x))]
        end
    end

    @testset "xcorr2" begin
        # setting the length of `u` and `v` to 6 forces the desired nfft to be 11,
        # which gets bumped to 12 because of `nextfastfft`, so these test that all
        # the index twiddling still holds.
        @testset "basic cross-correlation" begin
            @test xcorr2([1; zeros(5)],
                         [1; zeros(5)];
                         center=false,
                         unbiased=false) ≈ [zeros(5); 1; zeros(5)]
            @test xcorr2([zeros(5); 1],
                         [1; zeros(5)];
                         center=false,
                         unbiased=false) ≈ [zeros(10); 1]
            @test xcorr2([1; zeros(5)],
                         [zeros(5); 1];
                         center=false,
                         unbiased=false) ≈ [1; zeros(10)]
        end
        @testset "xcorr with lagbounds" begin
            u = [1; zeros(5)]
            v = [1; zeros(5)]
            @test xcorr2(u,v; center=false, lagbounds=(-1,1)) ≈ [0, 1, 0]
        end
        @testset "xcorr with short u" begin
            u = [1, 0]
            v = [1, 0, 0]
            @test xcorr2(u,v; center=false, unbiased=false) ≈ [0, 0, 1, 0]
        end
        @testset "xcorr with short v" begin
            u = [1, 0, 0]
            v = [1, 0]
            @test xcorr2(u,v; center=false, unbiased=false) ≈ [0, 1, 0, 0]
        end
    end

    @testset "padding" begin
        @test pad([1, 2, 3], :reflect_at, 2, 2) == [3,2,1,2,3,2,1]
    end

    @testset "resample" begin
        x = randn(1000)
        @test length(resample2(x, 2)) == 2000
        @test length(resample2(x, 0.5)) == 500

    end

    @testset "xc_align" begin
        x = [zeros(3); 1; zeros(7)]
        y = [zeros(5); 1; zeros(5)]
        vx, vy, os = xc_align(x, y, (-3,3))
        @test vx ≈ vy
        @test os == 2

        vx, vy, os = xc_align(y, x, (-3,3))
        @test vx ≈ vy
        @test os == -2

        x = [zeros(3); 1; zeros(7)]
        y = [zeros(5); 1; zeros(7)]
        vx, vy, os = xc_align(x, y, (-3,3))
        @test vx ≈ vy
        @test os == 2

        vx, vy, os = xc_align(y, x, (-3,3))
        @test vx ≈ vy
        @test os == -2
    end

end

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
include("Trees.jl")

typealias Params Dict{Any,Any}

abstract RouteNode

# StringNodes store all static route parts
# For example, route "/foo/bar/" is comprised of 2 StringNodes, "foo" and "bar"
#
immutable StringNode <: RouteNode
    val::String
end
# Equality & Matching for building / traversing the routing table
isequal(s1::StringNode, s2::String)     = s1.val == s2
isequal(s1::StringNode, s2::StringNode) = s1.val == s2.val
ismatch(s1::StringNode, s2::String)     = isequal(s1, s2)

# Returns the name part of a RegEx or DataType dynamic route
# Will return `foo` in `<foo::String>`
#
getname(p::String) = Base.ismatch(r"^<([^\:]*)::", p) ? match(r"^<([^\:]*)::", p).captures[1] : nothing

# DynamicNodes
#
# Route search uses `regex` to validate matches,
# then adds key `name` to `req.state[:params]` with unique value from
# resource after it is ( optionally ) processed through `convert`.
#
# DynamicNodes are used for three kinds of nodes:
#
#   - NamedParamNodes:     "/foo/<bar>"
#            => matches    "/foo/*", adds `req.params[:bar]::String`
#
#   - RegexNodes:          "/foo/<bar::%[0-9][0-9]-[A-Z]{5}>"
#            => matches    "/foo/00-ABCDE", adds `req.params[:bar]::String`
#            => wont match "/foo/00-ABCDEF", "foo/0-ABCDE"
#
#   - DataTypeNodes:       "/foo/<bar::Int>"
#            => matches    "/foo/10", adds `req.params[:bar]::Int`
#            => wont match "/foo/10.0", "foo/abc"
#
immutable DynamicNode <: RouteNode
    name::String
    regex::Regex
    convert::Union(Nothing, Function)
end
DynamicNode(n::String, r::Regex)    = DynamicNode(n,r,nothing)
DynamicNode(p::String, c::Function) = DynamicNode(getname(p), Regex("^$(match(r":%([^>]*)>", p).captures[1])\$"), c)

# Equality & Matching for building / traversing the routing table
isequal(s1::DynamicNode, s2::String)      = Base.ismatch( s1.regex, s2 )
isequal(s1::DynamicNode, s2::DynamicNode) = s1.regex.pattern == s2.regex.pattern && s1.name == s2.name
ismatch(s1::DynamicNode, s2::String)      = isequal( s1, s2 )

# NamedParam & Regex node constructors, these are NOT types.
NamedParamNode(p::String) = DynamicNode(p[2:length(p)-1], r".*")
RegexNode(p::String)      = DynamicNode(getname(p), Regex("^$(match(r":%([^>]*)>", p).captures[1])\$"))

# All valid DataTypeNode types -- provides `regex` validator + converter.
# Should be exposed to users for extension with custom DataTypes?
#
const DataTypeNodeBuilders = (String => (Regex, Function))[
    "Int"   => (r"^[0-9]*$", int),
    "Float" => (r"^[0-9]*\.[0-9]*$", float),
    "String"=> (r"^.*$", string)
]
# DataType node constructor, also NOT a type.
DataTypeNode(p::String, t::String) = DynamicNode(getname(p), DataTypeNodeBuilders[t]...)

# Build the params dataset – dispatch needed for both route types.
extend_params(params::Params, v::RouteNode, p::String) = params
function extend_params(params::Params, v::DynamicNode, p::String)
    params[v.name] = ((v.convert == nothing) ? p : v.convert(p))
    params
end

typealias Route (RouteNode,Union(Function,Nothing)) # ('/about', function()...)
isequal(r::Route, v) = isequal(r[1], v)
ismatch(r::Route, v) = ismatch(r[1], v)
isequal(node::RouteNode, route::Route) = isequal(node, route[1])

typealias RoutingTable Tree
RoutingTable() = RoutingTable((StringNode("/"), nothing))

ismatch(node::String, resource_chunk::String) = node == resource_chunk

# Regex matchers for determining DynamicRoute types =>
# Functions for building matching type from `part`.
#
const dynamic_route_dispatch = (Regex => Function)[
    r"^<[^:%]*::%"          => RegexNode,
    r"^<[^:>]*::[^%>]*>$"   => part -> DataTypeNode(part, match(r"^<[^:>]*::([^%>]*)>$", part).captures[1]),
    r"^<[^:>]*>$"           => NamedParamNode
]

# Run for each part of a route, builds RouteNodes.
function parse_part(part::String)
    if length(part) > 0 && Base.ismatch(r"^<[^>]*>$", part)
        for v in dynamic_route_dispatch
            if Base.ismatch(v[1], part)
                return v[2](part)
            end
        end
        throw("$part is an invalid route part.")
    end
    StringNode(part != "" ? part : "/")
end

# `path_to_handler` returns an array of `(RouteNode,Union(Nothing,Function))`
# pairs. Each element will hold a `nothing` except for the final element which
# will contain the handler function. e.g.:
#
#   path_to_handler("/hello/world", ()->"")
#
# returns:
#
#   (StringNode("hello"),nothing)
#   (StringNode("world"),# function)
#
function path_to_handler(route::String, handler::Function)
    path = Route[(parse_part(part),nothing) for part in split(strip(route, '/'), "/")]
    path[end] = (path[end][1], handler)
    path
end

# `register!` inserts a handler into a `RoutingTable`. If it is for the root
# resource, "/", then it overwrites the route node rather than inserting a new
# one.
#
function register!(table::RoutingTable, resource::String, handler::Function)
    path = path_to_handler(resource, handler)
    # NOTE: a bit hack-ey, but fixes the root routing problem
    if resource == "/"
        table.value = path[1]
    else
        insert!(table, path)
    end
end

# It is easiest to understand the behavior of `searchroute` by example.
# When passed an array of url/resource components, e.g. for "/hello/world" the
# array `["hello", "world"], `searchroute` returns a function, `searchpred` that
# takes a single argument.
#
# When passed matching components `searchpred` returns `false` until the final
# element of `paths` is matched, e.g.:
#
#   > searchpred = searchroute(["hello", "world"])
#   > searchpred("hello")
#   false
#   > searchpred("world")
#   true
#
# However if a non-matching component is passed to `search` it returns `PRUNE`:
#
#   > searchpred = searchroute(["hello", "world"])
#   > searchpred("goodbye")
#   PRUNE
#
# This is used to indicate that it is not neccessary to continue searching a
# given branch of the `RoutingTable`.
#
function searchroute(parts::Array, params::Params)
    function searchpred(val)
        if ismatch(val, parts[1])
            params = extend_params(params, val[1], parts[1])
            if length(parts) == 1
                true
            else
                parts = parts[2:end]; false
            end
        else
            PRUNE
        end
    end
end

# `match_route_handler` looks up a handler in `table` when given a route to a
# resource array form (e.g. "/hello/world" would be `["hello", "world"]`). If
# no match is found then it returns `nothing`.
#
function match_route_handler(table::RoutingTable, parts::Array)
    params = Params()
    result = search(table, searchroute(parts, params))
    ((result != nothing ? result[2] : nothing), params)
end


text:
"""
    GitHubPages(; assets::Vector{<:AbstractString}=String[]) -> GitHubPages

Add `GitHubPages` to a template's plugins to add [`Documenter`](@ref) support via GitHub
Pages, including automatic uploading of documentation from [`TravisCI`](@ref). Also
adds appropriate badges to the README, and updates the `.gitignore` accordingly.

# Keyword Arguments
* `assets::Vector{<:AbstractString}=String[]`: Array of paths to Documenter asset files.

!!! note
    If deploying documentation with Travis CI, don't forget to complete the required
    configuration (see
    [here](https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#SSH-Deploy-Keys-1)).
"""
struct GitHubPages <: Documenter
    gitignore::Vector{String}
    assets::Vector{String}

    function GitHubPages(; assets::Vector{<:AbstractString}=String[])
        for file in assets
            if !isfile(file)
                throw(ArgumentError("Asset file $(abspath(file)) does not exist"))
            end
        end
        # Windows Git recognizes these paths as well.
        new(["/docs/build/", "/docs/site/"], abspath.(assets))
    end
end

function badges(::GitHubPages, user::AbstractString, pkg_name::AbstractString)
    return [
        format(Badge(
            "Stable",
            "https://img.shields.io/badge/docs-stable-blue.svg",
            "https://$user.github.io/$pkg_name.jl/stable"
        )),
        format(Badge(
            "Dev",
            "https://img.shields.io/badge/docs-dev-blue.svg",
            "https://$user.github.io/$pkg_name.jl/dev"
        )),
    ]
end

function gen_plugin(p::GitHubPages, t::Template, pkg_name::AbstractString)
    invoke(gen_plugin, Tuple{Documenter, Template, AbstractString}, p, t, pkg_name)

    if haskey(t.plugins, TravisCI)
        docs_src = joinpath(t.dir, pkg_name, "docs", "src")
        open(joinpath(dirname(docs_src), "make.jl"), "a") do file
            write(
                file,
                """

                deploydocs(;
                    repo="$(t.host)/$(t.user)/$pkg_name.jl",
                )
                """
            )
        end
    end
    return ["docs/"]
end

function interactive(::Type{GitHubPages})
    print("GitHubPages: Enter any Documenter asset files (separated by spaces) []: ")
    return GitHubPages(; assets=string.(split(readline())))
end


text:
if get(ARGS, 1, "Array") == "CuArray"
    using CUDA
    ArrayType = CUDA.CuArray
    CUDA.allowscalar(false)
else
    ArrayType = Array
end

using Test
using RootSolvers
using StaticArrays

include("test_helper.jl")

@testset "Convergence reached" begin
    for problem in problem_list
        FT = typeof(problem.x̃)
        for sol_type in [CompactSolution(), VerboseSolution()]
            for tol in get_tolerances(FT)
                for method in get_methods(problem.x_init, problem.x_lower, problem.x_upper, problem.f′)
                    if problem.x_init isa AbstractArray
                        sol = RootSolvers.find_zero.(Ref(problem.f), method, sol_type, tol)
                        converged = map(x -> x.converged, sol)
                        X_roots = map(x -> x.root, sol)
                        @test isbits(method)
                        @test all(converged)
                        @test eltype(X_roots) == eltype(problem.x_init)
                        @test all(X_roots .≈ problem.x̃)
                        test_verbose!(sol_type, sol, problem, tol, all(converged))
                    else
                        sol = find_zero(problem.f, method, sol_type, tol)
                        @test isbits(method)
                        @test sol.converged
                        @test sol.root isa FT
                        @test sol.root ≈ problem.x̃
                        test_verbose!(sol_type, sol, problem, tol, sol.converged)
                    end
                end
            end
        end
    end
end

@testset "Convergence not reached" begin
    for problem in problem_list
        FT = typeof(problem.x̃)
        for sol_type in [CompactSolution(), VerboseSolution()]
            for tol in get_tolerances(FT)
                for method in get_methods(problem.x_init, problem.x_lower, problem.x_upper, problem.f′)
                    if problem.x_init isa AbstractArray
                        sol = RootSolvers.find_zero.(Ref(problem.f), method, sol_type, tol, 1)
                        converged = map(x -> x.converged, sol)
                        X_roots = map(x -> x.root, sol)
                        @test isbits(method)
                        @test !any(converged)
                        @test eltype(X_roots) == eltype(problem.x_init)
                        test_verbose!(sol_type, sol, problem, tol, any(converged))
                    else
                        sol = find_zero(problem.f, method, sol_type, tol, 1)
                        @test isbits(method)
                        @test !sol.converged
                        @test sol.root isa FT
                        test_verbose!(sol_type, sol, problem, tol, sol.converged)
                    end
                end
            end
        end
    end
end

include("runtests_kernel.jl")


text:
mutable struct EvalError{T}
  err::T
  trace::StackTrace
end

EvalError(err) = EvalError(err, StackTrace())

# Stacktrace fails on empty traces
EvalError(err, bt::Vector{Ptr{Nothing}}) = EvalError(err, isempty(bt) ? [] : stacktrace(bt))

macro errs(ex)
  :(try
      $(esc(ex))
    catch e
      EvalError(isa(e, LoadError) ? e.error : e, stacktrace(catch_backtrace()))
    end)
end

errtrace(e::EvalError) = errtrace(e.err, e.trace)

errmsg(e::EvalError) = errmsg(e.err)

function Base.show(io::IO, err::EvalError)
  print(io, errmsg(err))
  println(io)
  for frame in errtrace(err)
    show(io, frame)
    println(io)
  end
end

function cliptrace(trace::StackTrace)
  ind = findlast(frame -> (frame.func == :include_string &&
                           frame.file == Symbol(joinpath(".", "loading.jl"))) ||
                          (frame.func == :render′ &&
                           endswith(string(frame.file), joinpath("display", "errors.jl"))), trace)
  trace[1:(ind == nothing ? end : ind - 1)]
end

highlights(trace::StackTrace) =
  @>> trace filter(x->x.line > 0) map(x->d(:file=>fullpath(string(x.file)),:line=>x.line)) unique

highlights(e::EvalError) = highlights(e.trace)

function locationshading(file)
  f, p = expandpath(file)

  # error in base
  occursin(r"^base(/|\\).*", f) && return ".dark"
  # error in package
  # should probably be a bit smarter with figuring out if thats the actually package code
  occursin(joinpath(homedir(), ".julia"), p) && return ".medium"
  # error in "user code"
  return ".bright"
end

function renderbt(trace::StackTrace)
  span(".error-trace",
       [div(".trace-entry $(locationshading(string(frame.file)))",
            [fade("in "),
             frame.linfo == nothing || frame.linfo isa Core.CodeInfo ?
               string(frame.func) :
               replace(sprint(Base.show_tuple_as_call, frame.linfo.def.name, frame.linfo.specTypes),
                       r"\(.*\)$" => ""),
             fade(" at "),
             render(Inline(), Copyable(baselink(string(frame.file), frame.line))),
             fade(frame.inlined ? " <inlined>" : "")])
        for frame in reverse(trace)])
end

rendererr(err) = strong(".error-description", err)

function render(::Editor, e::EvalError)
  header = errmsg(e.err)
  header = split(header, '\n')
  trace = cliptrace(errtrace(e))
  view = if isempty(trace)
    if length(header) == 1
      rendererr(header[1])
    else
      Tree(rendererr(header[1]), [rendererr(join(header[2:end], '\n'))])
    end
  else
    Tree(rendererr(header[1]), [rendererr(join(header[2:end], '\n')); renderbt(trace)])
  end

  d(:type => :error,
    :view => render(Inline(),
                    Copyable(view, string(e))),
    :highlights => highlights(trace))
end

render(::Console, e::EvalError) =
  @msg result(render(Editor(), e))

mutable struct DisplayError
  obj
  err
end

errmsg(d::DisplayError) =
  "Error displaying $(typeof(d.obj)): $(errmsg(d.err))"

function render′(display, obj)
  try
    render(display, obj)
  catch e
    render(display, EvalError(DisplayError(obj, e), stacktrace(catch_backtrace())))
  end
end

render′(x) =
  render′(Media.getdisplay(Media.primarytype(x)), x)


text:
# DAY 1 PROBLEM 2 ADVENT OF CODE
# Get the number of increases in depth by window.

depths = """156
176
175
176
183
157
150
153
154
170
162
167
170
188
190
194
196
198
202
203
187
189
194
213
216
217
224
217
216
224
237
251
259
260
292
312
313
319
315
316
324
330
331
346
361
388
373
363
338
342
337
331
333
328
329
335
336
334
332
333
337
338
345
344
342
345
344
320
343
346
347
339
363
349
362
359
373
362
363
356
360
367
368
369
370
362
366
367
382
376
393
404
403
405
404
437
447
448
447
448
456
457
459
461
459
462
450
452
453
461
469
464
470
474
478
497
500
524
527
522
525
526
524
527
547
548
549
532
518
555
558
556
575
586
585
596
592
593
583
584
598
604
626
629
635
636
641
644
646
639
643
617
616
617
619
629
630
625
624
628
636
638
639
658
665
670
667
641
645
653
680
689
665
640
633
635
641
632
635
637
638
641
653
661
671
679
678
675
677
681
682
703
704
705
710
712
719
722
709
710
712
715
721
719
720
750
739
717
714
715
716
723
722
720
724
723
717
708
725
723
709
715
722
711
694
695
702
700
724
721
720
749
748
749
754
756
763
773
753
755
754
758
761
762
760
762
765
768
788
812
813
814
828
833
848
851
853
858
863
865
866
864
863
865
866
879
888
889
912
917
910
912
919
924
913
927
928
927
931
939
943
946
948
952
951
952
965
968
970
964
946
947
952
954
962
965
967
981
982
992
993
995
997
987
990
1005
1018
1000
999
1000
1001
983
978
1004
999
1002
1024
1022
1030
1037
1061
1071
1081
1083
1077
1076
1086
1101
1116
1137
1153
1158
1160
1167
1168
1184
1204
1214
1250
1271
1273
1272
1274
1297
1303
1304
1331
1349
1350
1355
1356
1357
1358
1371
1391
1396
1395
1399
1367
1368
1330
1327
1335
1336
1340
1331
1332
1329
1337
1324
1325
1326
1327
1313
1326
1328
1331
1307
1309
1310
1320
1324
1325
1324
1353
1354
1356
1357
1359
1360
1361
1357
1360
1362
1366
1373
1380
1399
1403
1414
1435
1439
1424
1433
1439
1446
1443
1444
1438
1435
1434
1426
1441
1451
1452
1463
1472
1468
1461
1489
1485
1487
1492
1489
1492
1489
1497
1468
1469
1473
1477
1476
1477
1502
1499
1501
1506
1502
1505
1512
1543
1541
1556
1557
1555
1557
1556
1530
1548
1545
1552
1559
1569
1570
1555
1563
1562
1565
1574
1573
1574
1573
1563
1572
1576
1571
1573
1565
1570
1571
1572
1579
1590
1591
1595
1597
1600
1596
1597
1596
1605
1607
1606
1607
1610
1612
1614
1620
1625
1639
1636
1639
1634
1633
1624
1644
1645
1642
1651
1654
1658
1663
1665
1666
1668
1671
1676
1663
1656
1659
1656
1655
1650
1656
1657
1667
1673
1674
1675
1676
1680
1673
1674
1677
1683
1690
1692
1697
1693
1694
1701
1714
1689
1691
1692
1684
1686
1685
1683
1685
1686
1692
1706
1694
1692
1702
1703
1723
1728
1726
1728
1729
1730
1738
1737
1703
1701
1714
1729
1730
1728
1729
1730
1729
1707
1714
1755
1760
1786
1787
1807
1814
1815
1813
1808
1810
1794
1797
1801
1800
1807
1802
1803
1804
1819
1821
1819
1828
1816
1817
1852
1839
1840
1837
1838
1877
1876
1868
1879
1886
1885
1886
1859
1864
1872
1905
1907
1912
1915
1920
1936
1935
1943
1945
1968
1963
1954
1969
1968
1976
1977
1978
1963
1961
1966
1968
1983
1984
1983
1982
1954
1957
1959
1970
1980
1986
1996
1997
1964
1957
1958
1957
1969
1982
1988
1985
1986
1982
1981
1980
2007
1999
2002
2033
2039
2040
2028
2048
2049
2047
2046
2051
2049
2036
2028
2041
2038
2041
2044
2052
2058
2060
2063
2068
2071
2070
2076
2075
2055
2054
2049
2035
2038
2048
2071
2090
2097
2102
2086
2085
2083
2085
2090
2107
2106
2138
2144
2151
2142
2153
2154
2155
2156
2159
2155
2178
2181
2182
2189
2190
2194
2197
2198
2194
2198
2201
2202
2201
2199
2171
2169
2164
2178
2181
2186
2187
2188
2191
2190
2193
2211
2226
2241
2227
2253
2259
2281
2282
2283
2282
2264
2263
2278
2295
2299
2290
2291
2294
2307
2306
2315
2314
2316
2315
2340
2339
2350
2362
2365
2322
2339
2340
2353
2368
2363
2361
2370
2371
2376
2368
2369
2371
2372
2394
2408
2415
2417
2399
2405
2408
2410
2417
2418
2425
2426
2424
2423
2453
2456
2467
2475
2478
2480
2487
2503
2506
2502
2508
2494
2497
2502
2505
2511
2512
2509
2516
2511
2516
2523
2493
2500
2502
2506
2507
2483
2486
2484
2454
2456
2467
2470
2473
2485
2484
2491
2492
2496
2497
2518
2483
2484
2511
2516
2485
2486
2487
2482
2483
2484
2481
2508
2512
2503
2511
2542
2536
2538
2547
2540
2531
2532
2535
2526
2525
2552
2562
2563
2564
2563
2562
2561
2564
2565
2567
2568
2563
2564
2576
2577
2582
2573
2569
2564
2570
2568
2571
2582
2564
2566
2564
2579
2588
2587
2590
2592
2598
2611
2634
2653
2652
2648
2649
2650
2636
2641
2660
2661
2659
2661
2653
2654
2640
2638
2636
2635
2633
2635
2636
2656
2657
2658
2667
2679
2707
2697
2708
2702
2703
2714
2715
2722
2726
2730
2720
2724
2725
2733
2745
2748
2750
2746
2736
2737
2739
2740
2720
2729
2734
2726
2739
2734
2736
2735
2736
2735
2737
2727
2694
2692
2702
2703
2710
2713
2708
2710
2721
2731
2722
2723
2724
2742
2750
2745
2735
2745
2748
2778
2770
2772
2773
2772
2773
2768
2779
2778
2779
2802
2801
2809
2810
2814
2833
2835
2844
2859
2860
2861
2865
2870
2872
2870
2876
2875
2877
2883
2872
2871
2872
2891
2889
2890
2868
2866
2864
2867
2862
2863
2864
2867
2870
2874
2875
2895
2887
2886
2887
2886
2889
2890
2891
2890
2867
2888
2895
2896
2891
2894
2896
2886
2888
2892
2878
2871
2876
2877
2875
2874
2875
2887
2875
2901
2920
2923
2924
2925
2907
2908
2931
2932
2933
2937
2951
2952
2951
2947
2967
2966
2962
2972
2953
2941
2962
2988
2995
3014
3008
3009
3024
3034
3029
3038
3058
3059
3060
3037
3040
3039
3030
3046
3058
3059
3067
3069
3082
3086
3088
3089
3086
3097
3105
3116
3114
3148
3152
3154
3164
3196
3198
3200
3212
3217
3224
3247
3249
3250
3251
3256
3243
3247
3249
3250
3249
3254
3255
3256
3273
3274
3275
3273
3274
3276
3263
3271
3272
3290
3289
3290
3292
3279
3284
3273
3277
3284
3289
3292
3301
3300
3338
3340
3375
3384
3382
3375
3377
3373
3376
3375
3382
3385
3404
3400
3403
3415
3408
3420
3421
3422
3419
3423
3432
3434
3435
3436
3438
3431
3426
3434
3460
3462
3463
3438
3445
3446
3451
3448
3468
3460
3453
3459
3464
3465
3459
3465
3462
3464
3474
3489
3501
3502
3503
3506
3507
3515
3505
3508
3509
3508
3510
3494
3495
3494
3495
3496
3497
3500
3513
3512
3513
3514
3531
3533
3531
3523
3524
3528
3550
3549
3538
3532
3529
3535
3526
3529
3528
3563
3533
3534
3541
3543
3544
3543
3565
3576
3581
3622
3624
3630
3622
3628
3640
3649
3644
3645
3655
3654
3648
3649
3644
3642
3643
3668
3672
3688
3718
3724
3734
3737
3758
3757
3760
3761
3777
3775
3770
3767
3768
3767
3774
3771
3778
3776
3782
3790
3791
3803
3819
3824
3822
3828
3819
3805
3809
3808
3809
3815
3817
3819
3817
3821
3842
3846
3844
3842
3857
3858
3876
3869
3878
3848
3850
3840
3862
3861
3859
3869
3872
3865
3864
3859
3861
3854
3853
3854
3850
3851
3860
3866
3850
3846
3847
3864
3852
3859
3857
3860
3859
3860
3862
3877
3880
3881
3914
3911
3913
3923
3924
3934
3936
3927
3934
3939
3940
3942
3936
3937
3936
3956
3936
3937
3942
3943
3951
3954
3966
3978
3980
3992
3997
4001
3991
3990
3994
4000
4007
4010
4014
3990
3991
3992
3995
4005
4006
4009
4022
4023
4026
4003
4017
4018
4020
4021
4023
4026
4028
4027
4016
4015
4011
4001
4027
4031
4033
4045
4059
4065
4063
4064
4065
4070
4075
4076
4069
4052
4054
4024
4017
4018
4017
4018
4015
4011
4014
4015
4022
4009
4008
4009
3988
4010
4025
3990
3998
4002
4008
4006
4008
4015
4013
4010
4019
4035
4034
4036
4044
4073
4070
4067
4068
4086
4087
4088
4085
4084
4083
4086
4083
4082
4083
4080
4077
4049
4050
4052
4053
4057
4065
4068
4065
4083
4084
4101
4110
4114
4113
4111
4112
4116
4094
4101
4105
4104
4107
4106
4120
4138
4147
4148
4152
4157
4160
4161
4163
4155
4153
4156
4178
4179
4184
4203
4206
4208
4212
4196
4197
4169
4173
4159
4160
4161
4162
4163
4160
4162
4164
4165
4166
4168
4166
4167
4157
4166
4157
4159
4175
4173
4162
4153
4137
4142
4144
4145
4146
4152
4158
4173
4171
4167
4168
4173
4172
4173
4178
4181
4161
4164
4160
4153
4155
4154
4149
4160
4162
4178
4172
4174
4183
4184
4192
4198
4201
4204
4210
4211
4216
4222
4232
4233
4248
4250
4254
4255
4270
4275
4276
4282
4289
4287
4288
4290
4291
4319
4334
4352
4360
4367
4377
4384
4407
4408
4407
4419
4404
4407
4423
4425
4424
4428
4429
4451
4446
4444
4445
4449
4472
4474
4475
4486
4496
4470
4473
4478
4477
4478
4480
4491
4492
4487
4486
4494
4499
4521
4520
4532
4536
4510
4516
4526
4534
4515
4509
4510
4511
4472
4473
4483
4477
4478
4497
4522
4524
4550
4553
4558
4559
4543
4534
4533
4521
4525
4526
4514
4520
4515
4503
4506
4511
4529
4530
4536
4538
4551
4553
4558
4561
4545
4542
4546
4534
4533
4535
4539
4536
4532
4541
4545
4538
4535
4516
4514
4513
4538
4549
4556
4562
4563
4564
4568
4570
4575
4576
4591
4577
4589
4590
4571
4576
4574
4575
4576
4574
4577
4574
4575
4588
4580
4598
4595
4610
4603
4604
4605
4593
4594
4580
4583
4604
4613
4614
4631
4629
4647
4628
4615
4634
4635
4636
4640
4646
4652
4651
4653
4649
4643
4645
4642
4643
4639
4644
4655
4658
4659
4660
4668
4669
4670
4680
4692
4682
4686
4688
4691
4690
4697
4698
4700
4703
4709
4713
4718
4720
4728
4742
4744
4743
4751
4753
4741
4747
4750
4751
4733
4732
4734
4723
4724
4738
4739
4767
4773
4775
4778
4782
4786
4790
4754
4757
4769
4763
4768
4770
4773
4798
4818
4819
4820
4819
4820
4822
4830
4828
4830
4815
4821
4828
4826
4829
4828
4824
4823
4827
4809
4816
4818
4819
4821
4822
4833
4837
4834
4835
4838
4848
4855
4866
4864
4859
4879
4882
4880
4878
4860
4847
4850
4835
4836
4840
4827
4849
4848
4849
4855
4847
4842
4846
4849
4807
4809
4815
4842
4847
4863
4846
4845
4844
4846
4849
4861
4841
4861
4862
4881
4880
4890
4892
4895
4916
4917
4923
4943
4945
4944
4946
4957
4955
4943
4939
4955
4957
4961
4962
4970
4972
4973
4979
4964
4963
4960
4961
4955
4946
4951
4949
4962
4964
4968
4976
4965
4967
4961
4965
4966
4953
4923
4928
4929
4932
4934
4936
4951
4949
4950
4976
4975
4996
4999
5007
5011
5010
5004
5017
5022
5016
5015
5022
5016
5025
5040
5033
5037
5024
5025
5028
5030
5026
5025
5028
5029
5035
5034
5058
5061
5060
5065
5049
5056
5055
5050
5068
5069
5074
5090
5083
5068
5080
5079
5080
5081
5084
5099
5100
5115
5114
5119
5120
5128
5120
5131
5130
5120
5113
5126
5130
5156
5158
5159
5163
5162
5125
5126
5135
5162
5166
5173
5169
5177
5183
5199
5200
5188
5191
5192
5184
5193
5203
5206
5214
5219
5215
5230
5231
5208
5207
5208
5181
5161
5162
5164
5189
5190
5170"""

"""
    get_depths(depths<:AbstractString) ::Vector{<:Integer}

Get the list of depths.
"""
function get_depths(depths::AbstractString) ::Vector{<:Integer}
    return [parse(Int16, x) for x in split(depths, "\n")]
end

"""
    window(sequence::Vector{<:Any}, win_size::Integer, step::Integer=1) ::Vector{Vector{<:Any}}

Iterate over the sequence by a window of the given size and step.
"""
function window(sequence::Vector, win_size::Integer, step::Integer=1) ::Vector{Vector}
    # Adjust for 1 indexing.
    win_size -= 1
    # get last valid index we can get win_size items after.
    last::Integer = length(sequence) - win_size
    # get vector of item vectors.
    return [sequence[i:i+win_size]
            for i=1:step:last]
end

"""
    windowed_depth_increases(depths::Vector{Vector{<:Integer}}) ::Integer

Get the number of depth increases.
"""
function windowed_depth_increases(depths::Vector{Vector{<:Number}}) ::Number
    collector = 0
    for i in 2:length(depths)
        if sum(depths[i]) > sum(depths[i-1])
            collector += 1
        end
    end
    return collector
end

grouping_size = 3

depths |> get_depths |> (x) -> window(x, grouping_size) |> windowed_depth_increases


text:
include("ABCDQR.jl")

include("ABCDQR_simula_test.jl")
include("ABCDQR_kfilter_test.jl")
include("ABCDQR_kfilter_s_test.jl")
include("ABCDQR_em_starting_test.jl")
include("ABCDQR_em_test.jl")
include("ABCDQR_em_s_test.jl")


text: include("Arduino.jl")


text:
# generic APIs

## dense vector descriptor

mutable struct CuDenseVectorDescriptor
    handle::cusparseDnVecDescr_t

    function CuDenseVectorDescriptor(v::CuVector{T}) where {T}
        vec_ref = Ref{cusparseDnVecDescr_t}()
        cusparseCreateDnVec(vec_ref, length(v), v, T)
        obj = new(vec_ref[])
        finalizer(cusparseDestroyDnVec, obj)
        obj
    end
end

Base.unsafe_convert(::Type{cusparseDnVecDescr_t}, desc::CuDenseVectorDescriptor) = desc.handle

## sparse matrix descriptor

mutable struct CuSparseMatrixDescriptor
    handle::cusparseSpMatDescr_t
end

function CuSparseMatrixDescriptor(A::CuSparseMatrixCSR{T}) where {T}
    desc_ref = Ref{cusparseSpMatDescr_t}()
    cusparseCreateCsr(
        desc_ref,
        A.dims..., length(A.nzVal),
        A.rowPtr, A.colVal, A.nzVal,
        eltype(A.rowPtr), eltype(A.colVal), 'O', eltype(A.nzVal)
    )
    obj = CuSparseMatrixDescriptor(desc_ref[])
    finalizer(cusparseDestroySpMat, obj)
    return obj
end

Base.unsafe_convert(::Type{cusparseSpMatDescr_t}, desc::CuSparseMatrixDescriptor) = desc.handle

## SpMV

function mv!(
    transa::SparseChar,
    alpha::T,
    A::CuSparseMatrixCSR{T},
    X::CuVector{T},
    beta::T,
    Y::CuVector{T},
    index::SparseChar
) where {T}

    m,n = size(A)

    if transa == 'N'
        chkmvdims(X,n,Y,m)
    end
    if transa == 'T' || transa == 'C'
        chkmvdims(X,m,Y,n)
    end

    cusparseSpMV(
        handle(),
        transa,
        [alpha],
        CuSparseMatrixDescriptor(A),
        CuDenseVectorDescriptor(X),
        [beta],
        CuDenseVectorDescriptor(Y),
        T,
        CUSPARSE_MV_ALG_DEFAULT,
        CU_NULL
    )

    Y
end


text:
using AWS.EC2
using AWS

include("config.jl")

instances = ec2_instances_by_tag("ClusterName", ec2_clustername)
hosts = ec2_hostnames(instances)

println("Running hosts in cluster $ec2_clustername")
println(hosts)


text:
# Autogenerated wrapper script for Pango_jll for x86_64-unknown-freebsd
export libpango, libpangocairo, libpangoft

using FriBidi_jll
using FreeType2_jll
using Glib_jll
using Fontconfig_jll
using HarfBuzz_jll
using Cairo_jll
JLLWrappers.@generate_wrapper_header("Pango")
JLLWrappers.@declare_library_product(libpango, "libpango-1.0.so.0")
JLLWrappers.@declare_library_product(libpangocairo, "libpangocairo-1.0.so.0")
JLLWrappers.@declare_library_product(libpangoft, "libpangoft2-1.0.so.0")
function __init__()
    JLLWrappers.@generate_init_header(FriBidi_jll, FreeType2_jll, Glib_jll, Fontconfig_jll, HarfBuzz_jll, Cairo_jll)
    JLLWrappers.@init_library_product(
        libpango,
        "lib/libpango-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpangocairo,
        "lib/libpangocairo-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpangoft,
        "lib/libpangoft2-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
base_vals = Base.ImmutableDict(DNA_A=>0,DNA_C=>1,DNA_G=>2,DNA_T=>3)

function get_mer_idx(mer,k=5)
    idx = 0
    @turbo for i in 1:k
        idx = 4*idx + base_vals[mer[i]]
    end
    return idx +1
end

function count_mers_4_dist(seq,k = 5)
    counts = zeros(Int,4^k)
    for mer in each(DNAMer{k}, seq)
        current_mer = mer.fw
        current_mer_idx = get_mer_idx(current_mer)
        counts[current_mer_idx] += 1
    end
    return counts
end
function count_mers_4_dist!(seq,mer_array,i,k = 5)
    for mer in each(DNAMer{k}, seq)
        current_mer = mer.fw
        current_mer_idx = get_mer_idx(current_mer)
        mer_array[current_mer_idx,i] += 1
    end
end
function kord_distance(n₁,n₂,k = 5)
    l1,l2 = length(n₁), length(n₂)
    if l1 == l2
        dist = 0.0
        step = k-1
        for i in 1:l1-step
            dist += n₁[i:i+step] == n₂[i:i+step]
        end
        return 1-dist/(l1-step)
    else
        return -1.0
    end
end

kmer_distance(n₁,n₂,L1,L2,k=5) = 1-sum(min.(n₁, n₂))/(min(L1,L2) - k + 1)
function kmer_distance(derep,k = 5)
    seqs = derep.sequence
    n = length(seqs)
    mer_array = zeros(Int,4^k,n)
    for i in 1:n
        count_mers_4_dist!(seqs[i],mer_array,i)
    end
    dist_mat = Array{Float64}(undef,n,n)
    dist_mat = zeros(n,n)
    for i in 1:n
        for j in 1:i
            dist_mat[i,j] = kmer_distance(mer_array[:,i],mer_array[:,j],length(seqs[i]),length(seqs[j]),k)
        end
    end
    return dist_mat
end
function kmer_count(derep,k = 5)
    seqs = derep.sequence
    n = length(seqs)
    mer_array = zeros(Int,4^k,n)
    for i in 1:n
        count_mers_4_dist!(seqs[i],mer_array,i)
    end

    return mer_array
end


text:
using QuantEcon
using PyPlot

function get_simultaneous_P(epsilon, N, p)
    P = Array(Float64, (N+1, N+1))
    for i in 1:(N+1)
        if i - 1 < p * N
            P[i, 1:end] = [binomial(N, j)*(epsilon)^j*(1-epsilon)^(N-j) for j in 0:N]
            P[i, end] = 1-sum(P[i, 1:end-1])
        elseif i - 1 == p * N
            P[i, 1:end] = [binomial(N, j)*(1/2)^N for j in 0:N]
            P[i, end] = 1-sum(P[i, 1:end-1])
        else
            P[i, 1:end] = [binomial(N, j)*(epsilon)^(N-j)*(1-epsilon)^j for j in 0:N]
            P[i, end] = 1-sum(P[i, 1:end-1])
        end
    end
    return P
end

function get_sequential_P(epsilon, N, p)
    P = zeros(Float64, (N+1, N+1))
    for i in 1:N+1
        if i == 1
            P[i, i + 1] = (N - (i-1))/N*epsilon/2
            if (i-1) > p * (N-1)
                P[i, i + 1] += 1 - epsilon
            elseif (i-1) == p * (N-1)
                P[i, i + 1] += 1/2
            end
            P[i, i] = 1 - P[i, i + 1]
        elseif i == N+1
            P[i, i-1] = (i-1)/N*epsilon/2
            if i - 2 == p * (N - 1)
                P[i, i-1] = 1/2
            elseif i - 2 < p * (N - 1)
                P[i, i-1] = 1 - epsilon
            end
            P[i, i] = 1 - P[i, i-1]
        else
            P[i, i+1] = (N - (i-1))/N*epsilon/2
            if 1 > p * (N - (i-1))
                P[i, i+1] += 1 - epsilon
            elseif 1 == p * (N - (i-1))
                P[i, i+1] += 1/2
            end
            P[i, i-1] = (i-1)/N*epsilon/2
            if i - 2 == p * (N - 1)
                P[i, i-1] = 1/2
            elseif i - 2 < p * (N - 1)
                P[i, i-1] = 1 - epsilon
            end
            P[i, i] = 1 - P[i, i-1] - P[i, i + 1]
        end
    end
    return P
end

type KMR
    mc
    N
    epsilon
    Xs
end

get_p_2x2(M) = (M[1, 1, 1] - M[2, 1, 1])/(M[1, 1, 1] - M[1, 2, 1] - M[2, 1, 1] + M[2, 2, 1])
KMR(M, N::Int, epsilon::Float64; simul=true) = KMR(simul ? MarkovChain(get_simultaneous_P(epsilon, N, get_p_2x2(M))) : MarkovChain(get_sequential_P(epsilon, N, get_p_2x2(M))), N, epsilon, None)

function simulate!(kmr::KMR, T; init=0)
    kmr.Xs = simulate(kmr.mc, T, init+1) .- 1
end

plot_sample_path(kmr::KMR) = plot(kmr.Xs)#, label="epsilon = $kmr.epsion")

function plot_stationary_dist(kmr::KMR)
    fig, ax = subplots()
    #ax[:hist](mc_compute_stationary(kmr.mc))
    ax[:bar](collect(0:kmr.N), mc_compute_stationary(kmr.mc))
end

function plot_empirical_dist(kmr::KMR)
    fig, ax = subplots()
    ax[:hist](kmr.Xs, normed=true)
end

function plot_sample_path(kmr::KMR)
    plot(kmr.Xs)
end

"""
M = Array(Int, (2, 2, 2))
M[1, 1, :] = [4, 4]
M[1, 2, :] = [0, 3]
M[2, 1, :] = [3, 0]
M[2, 2, :] = [2, 2]

N = 20
T = 100
epsilon = 1/5

kmr = KMR(M, N, epsilon)
simulate!(kmr, T)
#plot_stationary_dist(kmr)
plot_empirical_dist(kmr)
#plot_sample_path(kmr)

for e in linspace(1/2, 0, 7)
    kmr = KMR(M, N, e)
    simulate!(kmr, T)
    #plot_sample_path(kmr)
    #plot_empirical_dist(kmr)
    plot_stationary_dist(kmr)
end
#println(Xs)
#plot(Xs)
#println(P)
#println(mean(X .== 2))
"""


text:
module Minimise_Energy

using JuLIP
using JuLIP.Potentials
using Main.Fast_Hertz
using LinearAlgebra
#using Pressure
import JuLIP: energy, forces, hessian_pos

export min_Energy

function min_Energy(at::AbstractAtoms)
    max_iter_num = 1000

    α = get_data(at, :Alpha)
    k = 1.0

    r = at.M
    rmin = minimum(r)
    rmax = maximum(r)
    r0 = 2*rmax
    r1 = 2*rmin

    V = Fast_Hertz.Fast_Hertz_Potential(get_data(at, :Length), r0, r1, α)

    set_calculator!(at, V)

    myresult = minimise!(at, precond = :id, method = :lbfgs, gtol=1e-13, verbose=2)
    is_converged = myresult.f_converged || myresult.g_converged || myresult.x_converged
    optmsg = myresult

    iter_num = myresult.iterations
    f_calls = myresult.f_calls
    g_calls = myresult.g_calls

    while  (myresult.iterations == max_iter_num) & !(is_converged)
        myresult = minimise!(at, precond = :id, method = :lbfgs, gtol=1e-13, verbose=2)
        iter_num += myresult.iterations
        f_calls += myresult.f_calls
        g_calls += myresult.g_calls
        @show iter_num
    end
    optmsg = myresult
    optmsg.iterations = iter_num
    optmsg.f_calls = f_calls
    optmsg.g_calls = g_calls

    @show optmsg.iterations
    set_data!(at, :Optmsg, optmsg)

    Energy = energy(V,at)
    set_data!(at, :Energy, Energy)

    # p = Pressure.cal_Pressure(at)
    # set_data!(at, :Stress, p)

    return at

end

end


text:
module DemoPackageBXW
print_greeting() = printstyled("Hello, world!"; color=:purple)
# Write your package code here.

end


text:
import Base: broadcast
import Base.Broadcast.instantiate
using Base.Broadcast: BroadcastStyle, Broadcasted, combine_axes, check_broadcast_axes, instantiate, flatten, broadcastable, _broadcast_getindex, newindex
using SparseArrays

struct AssocStyle <: Broadcast.BroadcastStyle end

Broadcast.BroadcastStyle(::Type{<:Assoc}) = AssocStyle()
Broadcast.BroadcastStyle(::Type{<:Assoc}, ::Type{<:Assoc}) = AssocStyle()

Base.broadcastable(A::Assoc) = A

function Base.similar(bc::Broadcast.Broadcasted{AssocStyle}, ::Type{ElType}) where ElType

    # Scan the inputs for the ArraAssocyAndChar:
    #print(flatten(bc).args)

    ABrow,ABcol = rowcolunion(bc.args...)

    #=
    A,B = flatten(bc).args
    #First, create the row and col of the intersection
    ABrow = sortedintersect(A.row,B.row)
    ABcol = sortedintersect(A.col,B.col)
    #Filling the sparse matrix with

    rowMapping = searchsortedmapping(ABrow,A.row)
    colMapping = searchsortedmapping(ABcol,A.col)
    AA = A.A[rowMapping,colMapping]
    #AA = round(Int64,AA)

    BB = spzeros(size(ABrow,1), size(ABcol,1))
    rowMapping = searchsortedmapping(ABrow,B.row)
    colMapping = searchsortedmapping(ABcol,B.col)
    BB = B.A[rowMapping,colMapping]
    #BB = round(Int64,BB)
    =#
    AA = spzeros(size(ABrow,1), size(ABcol,1))
    # Use the char field of A to create the output
    Assoc(ABrow,ABcol,promote([1.0],ABrow)[1], similar(AA, Base.axes(AA)))
end

#=
function rowcolintersect(A::Assoc,B::Assoc)
    ABrow = sortedintersect(A.row,B.row)
    ABcol = sortedintersect(A.col,B.col)

    return ABrow,ABcol

end
=#

function rowcolunion(A::Assoc,B::Assoc)
    ABrow = sortedunion(A.row,B.row)
    ABcol = sortedunion(A.col,B.col)

    return ABrow,ABcol

end

function rowcolunion(A::Assoc,ABrow,ABcol)
    ABrow = sortedunion(A.row, ABrow)
    ABcol = sortedunion(A.col, ABcol)

    return ABrow,ABcol

end
rowcolunion(A::Assoc,bc::Broadcast.Broadcasted{AssocStyle}) = rowcolunion(A,rowcolunion((bc).args...)...)
rowcolunion(bc::Broadcast.Broadcasted{AssocStyle},B) = rowcolunion(combinedims(bc).args[2],B)

function combinedims(bc::Broadcast.Broadcasted{AssocStyle})
    ABrow,ABcol = rowcolunion(bc.args...)

    return combinedims(bc,ABrow,ABcol)
end

function combinedims(bc::Broadcast.Broadcasted{AssocStyle},ABrow,ABcol)

    if isa(bc.args[1],Broadcast.Broadcasted{AssocStyle})
        A = combinedims(bc.args[1],ABrow,ABcol)
    else

        newAA =  spzeros(length(ABrow),length(ABcol))
        Arow = searchsortedmapping(bc.args[1].row,ABrow)
        Acol = searchsortedmapping(bc.args[1].col,ABcol)
        newAA[Arow,Acol] = bc.args[1].A
        A = Assoc(ABrow,ABcol, Array{Union{AbstractString,Number}}([1.0]), newAA)

    end

    if isa(bc.args[2],Broadcast.Broadcasted{AssocStyle})
        B = combinedims(bc.args[2])
    else

        newBA =  spzeros(length(ABrow),length(ABcol))
        Brow = searchsortedmapping(bc.args[2].row,ABrow)
        Bcol = searchsortedmapping(bc.args[2].col,ABcol)
        newBA[Brow,Bcol] = bc.args[2].A
        B = Assoc(ABrow,ABcol, Array{Union{AbstractString,Number}}([1.0]), newBA)
    end

    return Broadcast.broadcasted(bc.f,A,B)
end

function Base.Broadcast.instantiate(bc::Broadcast.Broadcasted{AssocStyle})
    if bc.axes isa Nothing # Not done via dispatch to make it easier to extend instantiate(::Broadcasted{Style})
        bc = combinedims(bc)
        axes = combine_axes(bc.args...)
    else
        axes = bc.axes
        check_broadcast_axes(axes, bc.args...)
    end
    return Broadcasted{AssocStyle}(bc.f, bc.args, axes)
end

Base.@propagate_inbounds Base.Broadcast._broadcast_getindex(A::Assoc, I) = A.A[Base.Broadcast.newindex(A, I)]

function Base.setindex!(A::Assoc, val, inds)
    if ~isnan(val) && val != 0
        A.A[inds[1],inds[2]] = val
    end
end


text:
using BayesOpt
using Base.Test

# simple quadratic function
function testfunc(n, x, grad, fdata)
    f = 0.0
    let x = pointer_to_array(x, n)
        for i in 1:n
            f += (x[i] - .53)^2
        end
    end
    return f::Cdouble
end

# low-level interfaces (not exported)
let
    n = 10
    fdata = C_NULL
    lb = ones(n) * 0
    ub = ones(n) * 1
    xmin = similar(lb)
    fmin = similar(lb)
    params = BayesOpt.initialize_parameters_to_default()
    @test isa(params, BayesOpt.bopt_params)
    params.verbose_level = 0
    r = BayesOpt.bayes_optimization(n, testfunc, fdata, lb, ub, xmin, fmin, params)
    @test r == 0
    @test vecnorm(xmin - 0.53) < 0.1
end

# high-level interfaces (exported)
if VERSION > v"0.5-"
    let
        n = 10
        x0 = rand(n)
        lb = ones(n) * 0
        ub = ones(n) * 1
        fmin, xmin = bayesopt(
            x -> vecnorm(x - .53),
            x0,
            lb,
            ub,
            n_iterations=100,
            verbose_level=0)
        @test fmin < 0.1
        @test vecnorm(xmin - 0.53) < 0.1
    end
end


text:
using Distributed
using ThreadPools
using JLD2
using LinearMaps
using Printf

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

function inadmissible_error(eqind, eqn, point, val)
    vars = tuple(eqn.vsyms...)
    args = tuple((eval(Expr(:(=), sym, val)) for (sym, val) in zip(eqn.vsyms, point[eqn.vinds]))...)
    if typeof(val) <: AbstractArray
        @error "Singular gradient in equation $eqind: $eqn" vars args val
        error("Singular gradient in equation $eqind: $eqn")
    else
        @error "Inadmissible point in equation $eqind: $eqn" vars args val
        error("Inadmissible point in equation $eqind: $eqn")
    end
end

"""
    global_SS_R!(res, point, equations)

Compute the residual vector of the given set of equations at the given point.
The `equations` argument can be a container or an iterator whose `eltype` is `SteadyStateEquations`.
The residual `res` is updated in place.
"""
function global_SS_R!(resid::AbstractVector{Float64}, point::AbstractVector{Float64}, eqns::EqnIter) where EqnIter
    if ! (eltype(eqns) <: SteadyStateEquation)
        error("Expected a set of steady state equations, not $(EqnIter).")
    end
    for (eqind, eqn) in enumerate(eqns)
        rr = try eqn.eval_resid(point[eqn.vinds]) catch; NaN64 end
        if isnan(rr) || isinf(rr)
            inadmissible_error(eqind, eqn, point, rr)
            # vars = tuple(eqn.vsyms...)
            # args = tuple((eval(Expr(:(=), sym, val)) for (sym, val) in zip(eqn.vsyms, point[eqn.vinds]))...)
            # @error "Inadmissible point in equation $eqind: $eqn" vars args rr
            # error("Inadmissible point in equation $eqind: $eqn")
        end
        resid[eqind] = rr
    end
    return nothing
end
@assert precompile(global_SS_R!, (Vector{Float64}, Vector{Float64}, Vector{SteadyStateEquation}))

"""
    global_SS_R!(res, point, model::Model)

When a model is given, we compute the residual of the entire steady state system.
"""
@inline global_SS_R!(resid::AbstractVector{Float64}, point::AbstractVector{Float64}, model::Model) = global_SS_R!(resid, point, ModelBaseEcon.alleqns(model.sstate))
@assert precompile(global_SS_R!, (Vector{Float64}, Vector{Float64}, Model))

"""
    global_SS_RJ(point, equations)

Compute the residual vector and the Jacobian matrix of the given set of equations at the given point.
The return value is a tuple with two values.
The `equations` argument can be a container or an iterator whose `eltype` is `SteadyStateEquations`.
"""
function global_SS_RJ(point::AbstractVector{Float64}, eqns::EqnIter) where EqnIter
    if ! (eltype(eqns) <: SteadyStateEquation)
        error("Expected a set of steady state equations, not $(EqnIter).")
    end
    neqns = length(eqns)
    nvars = length(point)
    R = zeros(neqns)
    J = zeros(neqns, nvars)
    for (eqind, eqn) in enumerate(eqns)
        rr, jj = try eqn.eval_RJ(point[eqn.vinds]) catch; NaN64, fill(NaN64, size(eqn.vinds)) end
        if isnan(rr) || isinf(rr) || any(@. isnan(jj) | isinf(jj))
            inadmissible_error(eqind, eqn, point, rr)
            # args = tuple((eval(Expr(:(=), sym, val)) for (sym, val) in zip(eqn.vsyms, point[eqn.vinds]))...)
            # vars = tuple(eqn.vsyms...)
            # @error "Inadmissible point in equation $eqind: $eqn" vars args rr
            # error("Inadmissible point in equation $eqind: $eqn")
        end
        if any(@. isnan(jj) | isinf(jj))
            inadmissible_error(eqind, eqn, point, jj)
        #     args = tuple((eval(Expr(:(=), sym, val)) for (sym, val) in zip(eqn.vsyms, point[eqn.vinds]))...)
        #     vars = tuple(eqn.vsyms...)
        #     @error "Singular gradient in equation $eqind: $eqn" vars args rr
        #     error("Singular gradient in equation $eqind: $eqn")
        end
        R[eqind] = rr
        J[eqind,eqn.vinds] .= jj
    end
    return R, J
end
@assert precompile(global_SS_RJ, (Vector{Float64}, Vector{SteadyStateEquation}))

"""
    global_SS_RJ(point, model::Model)

When a model is given, we compute the residual of the entire steady state system.
"""
@inline global_SS_RJ(point::Vector{Float64}, model::Model) = global_SS_RJ(point, ModelBaseEcon.alleqns(model.sstate))
@assert precompile(global_SS_RJ, (Vector{Float64}, Model))