text:
if !isdefined(Main, :Groebner)
    import Groebner
end

include("biomodels/parser.jl")

import Singular
import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

BenchmarkTools.DEFAULT_PARAMETERS.samples = 3

function benchmark_system_my(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    @btime gb = Groebner.groebner($system)
    # gb = Groebner.groebner(system, reduced=false)
end

function benchmark_system_singular(system)
    R = AbstractAlgebra.parent(system[1])
    n = AbstractAlgebra.nvars(R)
    ground_s = Singular.QQ
    R_s, _ = Singular.PolynomialRing(ground_s, ["x$i" for i in 1:n], ordering=:degrevlex)

    system_s = map(
        f -> AbstractAlgebra.change_base_ring(
                    ground_s,
                    AbstractAlgebra.map_coefficients(c -> ground_s(numerator(c), denominator(c)), f),
                    parent=R_s),
        system)

    ideal_s = Singular.Ideal(R_s, system_s)

    Singular.std(Singular.Ideal(R_s, [system_s[1]]))

    @btime Singular.std($ideal_s)
end

function run_biomodels_benchmarks()
    systems = read_BIOMDs(5:20)

    println()
    for (name, system) in systems
        println("$name")
        println("my")
        benchmark_system_my(system)
        println("singular")
        benchmark_system_singular(system)
    end
end

run_biomodels_benchmarks()

#=

=#


text:
using PyPlot, MRIReco, FFTW, LinearAlgebra, BenchmarkTools

FFTW.set_num_threads(1);BLAS.set_num_threads(1)

threads = [1,2,4,8]
times = zeros(length(threads))

for t=1:length(times)
  @show threads[t]
  times[t] = parse(Float64,read(`julia-1.5 --threads=$(threads[t]) reconstruction.jl`, String))
end

plot(threads,times)

ext: jl


text:
using Wakame
using Documenter

DocMeta.setdocmeta!(Wakame, :DocTestSetup, :(using Wakame); recursive=true)

makedocs(;
    modules=[Wakame],
    authors="Bernard Brenyah",
    repo="https://github.com/PyDataBlog/Wakame.jl/blob/{commit}{path}#{line}",
    sitename="Wakame.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://PyDataBlog.github.io/Wakame.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/PyDataBlog/Wakame.jl",
    devbranch="main",
)


text:
# DataFrames
FILE = joinpath(DSTESTDIR, "randoms_small.csv")
DF = readtable(FILE)
if typeof(DF[:hiredate]) <: NullableVector
    DF[:hiredate] = NullableArray(Date[isnull(x) ? Date() : Date(get(x)) for x in DF[:hiredate]], [isnull(x) for x in DF[:hiredate]])
    DF[:lastclockin] = NullableArray(DateTime[isnull(x) ? DateTime() : DateTime(get(x)) for x in DF[:lastclockin]], [isnull(x) for x in DF[:lastclockin]])
else
    for i = 1:5
        T = eltype(DF.columns[i])
        DF.columns[i] = NullableArray(T[isna(x) ? (T <: String ? "" : zero(T)) : x for x in DF.columns[i]], [isna(x) for x in DF.columns[i]])
    end
    DF.columns[6] = NullableArray(Date[isna(x) ? Date() : Date(x) for x in DF[:hiredate]], [isna(x) for x in DF[:hiredate]])
    DF.columns[7] = NullableArray(DateTime[isna(x) ? DateTime() : DateTime(x) for x in DF[:lastclockin]], [isna(x) for x in DF[:lastclockin]])
end
DF2 = deepcopy(DF)
dfsource = Tester("DataFrame", x->x, false, DataFrame, (:DF,), scalartransforms, vectortransforms, x->x, x->nothing)
dfsink = Tester("DataFrame", x->x, false, DataFrame, (:DF2,), scalartransforms, vectortransforms, x->x, x->nothing)
function DataFrames.DataFrame(sym::Symbol; append::Bool=false)
    return @eval $sym
end
function DataFrames.DataFrame(sch::Data.Schema, ::Type{Data.Field}, append::Bool, ref::Vector{UInt8}, sym::Symbol)
    return DataFrame(DataFrame(sym), sch, Data.Field, append, ref)
end
function DataFrame(sink, sch::Data.Schema, ::Type{Data.Field}, append::Bool, ref::Vector{UInt8})
    rows, cols = size(sch)
    newsize = max(0, rows) + (append ? size(sink, 1) : 0)
    # need to make sure we don't break a NullableVector{WeakRefString{UInt8}} when appending
    if append
        for (i, T) in enumerate(Data.types(sch))
            if T <: Nullable{WeakRefString{UInt8}}
                sink.columns[i] = NullableArray(String[string(get(x, "")) for x in sink.columns[i]])
                sch.types[i] = Nullable{String}
            end
        end
    else
        for (i, T) in enumerate(Data.types(sch))
            if T != eltype(sink.columns[i])
                sink.columns[i] = NullableArray(eltype(T), newsize)
            end
        end
    end
    newsize != size(sink, 1) && foreach(x->resize!(x, newsize), sink.columns)
    if !append
        for (i, T) in enumerate(eltypes(sink))
            if T <: Nullable{WeakRefString{UInt8}}
                sink.columns[i] = NullableArray{WeakRefString{UInt8}, 1}(Array{WeakRefString{UInt8}}(newsize), fill(true, newsize), isempty(ref) ? UInt8[] : ref)
            end
        end
    end
    sch.rows = newsize
    return sink
end

# SQLite
dbfile = joinpath(dirname(@__FILE__), "randoms.sqlite")
dbfile2 = joinpath(dirname(@__FILE__),"test2.sqlite")
db = SQLite.DB(dbfile)
SQLite.load(db, "randoms_small", DF)
cp(dbfile, dbfile2; remove_destination=true)
db2 = SQLite.DB(dbfile2)
SQLite.createtable!(db2, "randoms2_small", Data.schema(SQLite.Source(db2, "select * from randoms_small")))
sqlitesource = Tester("SQLite.Source", SQLite.query, true, SQLite.Source, (db2, "select * from randoms_small"), scalartransforms, vectortransforms, x->x, ()->nothing)
sqlitesink = Tester("SQLite.Sink", SQLite.load, true, SQLite.Sink, (db2, "randoms2_small"), scalartransforms, vectortransforms, x->SQLite.query(db2, "select * from randoms2_small"), (x,y)->nothing)

DataStreamsIntegrationTests.teststream([dfsource, sqlitesource], [dfsink, sqlitesink]; rows=99)
finalize(db)
finalize(db2)
db = db2 = nothing;
gc(); gc();
rm(dbfile)
rm(dbfile2)

max_stars_repo_name: beffaCo/updateSqlite

max_issues_repo_name: beffaCo/updateSqlite

max_forks_repo_name: beffaCo/updateSqlite


text:
using ParetoSmoothedImportanceSampling
using StanSample, StatsFuns
using Printf
using JSON
using Test

ProjDir = @__DIR__

include(joinpath(ProjDir, "cvit.jl"))

# Data
data = JSON.parsefile(joinpath(ProjDir, "wells.data.json"))
y = Float64.(data["switched"])
x = Float64[data["arsenic"]  data["dist"]]
n, m = size(x)

# Model
model_str = read(open(joinpath(ProjDir, "arsenic_logistic.stan")), String)
sm1 = SampleModel("arsenic_logistic", model_str)

data1 = (p = m, N = n, y = Int.(y), x = x)
# Fit the model in Stan
rc1 = stan_sample(sm1; data=data1)
if success(rc1)
    nt1 = read_samples(sm1)

    # Compute LOO and standard error
    log_lik = nt1.log_lik'
    loo, loos, pk = psisloo(log_lik)
    elpd_loo = sum(loos)
    se_elpd_loo = std(loos) * sqrt(n)
    @printf(">> elpd_loo = %.1f, SE(elpd_loo) = %.1f\n", elpd_loo, se_elpd_loo)

    # Check the shape parameter k of the generalized Pareto distribution
    if all(pk .< 0.5)
        println("All Pareto k estimates OK (k < 0.5)")
    else
        pkn1 = sum((pk .>= 0.5) & (pk .< 1))
        pkn2 = sum(pk .>= 1)
        @printf(">> %d (%.0f%%) PSIS Pareto k estimates between 0.5 and 1\n", pkn1, pkn1/n*100)
        @printf(">> %d (%.0f%%) PSIS Pareto k estimates greater than 1\n", pkn2, pkn2/n*100)
    end
end
println()

# Fit a second model, using log(arsenic) instead of arsenic
x2 = Float64[log.(data["arsenic"])  data["dist"]]

# Model
data2 = (p = m, N = n, y = Int.(y), x = x2)
# Fit the model in Stan
rc2 = stan_sample(sm1; data=data2)

if success(rc2)
    nt2 = read_samples(sm1)
    # Compute LOO and standard error
    log_lik = nt2.log_lik'
    loo2, loos2, pk2 = psisloo(log_lik)
    elpd_loo = sum(loos2)
    se_elpd_loo = std(loos2) * sqrt(n)
    @printf(">> elpd_loo = %.1f, SE(elpd_loo) = %.1f\n", elpd_loo, se_elpd_loo)

    # Check the shape parameter k of the generalized Pareto distribution
    if all(pk .< 0.5)
        println("All Pareto k estimates OK (k < 0.5)")
    else
        pkn1 = sum((pk .>= 0.5) & (pk .< 1))
        pkn2 = sum(pk .>= 1)
        @printf(">> %d (%.0f%%) PSIS Pareto k estimates between 0.5 and 1\n", pkn1, pkn1/n*100)
        @printf(">> %d (%.0f%%) PSIS Pareto k estimates greater than 1\n", pkn2, pkn2/n*100)
    end
end

if success(rc1) && success(rc2)
    ## Compare the models
    loodiff = loos - loos2
    @printf("elpd_diff = %.1f, SE(elpd_diff) = %.1f\n",sum(loodiff), std(loodiff) * sqrt(n))
end
println()

## k-fold-CV
# k-fold-CV should be used if several khats>0.5
# in this case it is not needed, but provided as an example

model_str = read(open(joinpath(ProjDir, "arsenic_logistic_t.stan")), String)
sm3 = SampleModel("arsenic_logistic_t", model_str);

cvitr, cvitst = cvit(n, 10, true)
kfcvs = similar(loos)
for cvi in 1:3
    standatacv = (p = m, N = length(cvitr[cvi]), Nt = length(cvitst[cvi]),
        x = x[cvitr[cvi],:], y = Int.(y[cvitr[cvi]]),
        xt = x[cvitst[cvi],:], yt = Int.(y[cvitst[cvi]]))

    # Fit the model in Stan
    rc3 = stan_sample(sm3; data=standatacv)
    if success(rc3)
        nt3 = read_samples(sm3)
        # Compute LOO and standard error
        log_likt = nt3.log_likt'
        local n_sam, n_obs = size(log_likt)
        kfcvs[cvitst[cvi]] .=
            reshape(logsumexp(log_likt .- log(n_sam), dims=1), n_obs)
    end
end


text:
##############################################################################
#
# RedistributionFunction.jl
#
# Part of cvortex.jl
# Get vortex redistribution methods.
#
# Copyright 2019 HJA Bird
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
##############################################################################

"""
A vortex redistribution function.

The method by which vortex particles are redistibuted in lagrangian vortex
particles onto a grid to ensure long term stability.

# Examples
```
using CVortex
my_redistribution_kernel = singular_regularisation()
my_redistribution_kernel = winckelmans_regularisation()
my_redistribution_kernel = planetary_regularisation()
my_redistribution_kernel = gaussian_regularisation()
```
"""
struct RedistributionFunction
	func:: Ptr{Cvoid}			# Actually float(*g_3D)(float rho)
	critical_dist :: Float32	# Char[32]
end

#= Functions to to get VortexFunc structures =#

"""
Implements 0th order redistribution method.
The use of this redistribution is not advised.
lambda3_redistribution or m4p_redistribution are advised.
returns RedistributionFunction struct.
"""
function lambda0_redistribution()
	ret = ccall(("cvtx_RedistFunc_lambda0", libcvortex),
        RedistributionFunction, ())
	return ret;
end
"""
Implements 1st order redistribution method.
This method is dissipative.
returns RedistributionFunction struct.
"""
function lambda1_redistribution()
	ret = ccall(("cvtx_RedistFunc_lambda1", libcvortex),
        RedistributionFunction, ())
	return ret;
end
"""
Implements 2nd order redistribution method.
This redistribution is discontinuous, and consequently numerically
less useful than the lambda3_redistribution or m4p_redistribution.
returns RedistributionFunction struct.
"""
function lambda2_redistribution()
	ret = ccall(("cvtx_RedistFunc_lambda2", libcvortex),
        RedistributionFunction, ())
	return ret;
end
"""
Implements 3th order redistribution method.
returns RedistributionFunction struct.
"""
function lambda3_redistribution()
	ret = ccall(("cvtx_RedistFunc_lambda3", libcvortex),
        RedistributionFunction, ())
	return ret;
end
"""
Implements M_4' method of vorticity redistribution.
returns RedistributionFunction struct.
"""
function m4p_redistribution()
	ret = ccall(("cvtx_RedistFunc_m4p", libcvortex),
        RedistributionFunction, ())
	return ret;
end


text:
"""
Distribute the weights over n (about) equal-sized bins by swapping between bins with high total weight and low total weight.
"""
function distribute_weights_swapping(weights,n;nswaps = length(weights)^2)
    # we first 'normally' distribute the weights (first k to first thread, second k to second thread etc.)
    # Then we try to improve it a number of times by swapping weights between sets with a high and low total weight
    step = div(length(weights),n)+1 # first number of sets have this size, second part has this size-1
    nstep = n-(step*n-length(weights)) # this is the number of sets with size 'step'
    sets = vcat([collect((i-1)*step+1:i*step) for i=1:nstep],
            [collect(nstep*step+(i-1)*(step-1)+1:nstep*step+i*(step-1)) for i=1:n-nstep])
    set_weights = [sum(weights[sets[i]]) for i=1:length(sets)] #the total weights

    # edge case: less than 1 weight for every core
    if length(weights) <= n
        return sets,set_weights,[weights[s] for s in sets]
    end

    index_set = 1
    index_el = 1
    for k=1:nswaps # nswaps is actually not the number of swaps, but the number of tries
        # The current (large) set and the current element of the set
        max_set = sort([(set_weights[i],i) for i=1:length(set_weights)],rev=true)[index_set][2]
        max_el = sets[max_set][sort([(weights[sets[max_set]][i],i) for i=1:length(weights[sets[max_set]])],rev=true)[index_el][2]]
        # The small set and the small element of the set
        min_set = argmin(set_weights)
        min_el = sets[min_set][argmin(weights[sets[min_set]])]
        #Determine whether the swap works: the min_set should not increase by too much, and the max should of course also not increase
        if set_weights[min_set]+weights[max_el]-weights[min_el] < set_weights[max_set] &&
            set_weights[max_set]-weights[max_el]+weights[min_el] < set_weights[max_set]
            # swapping distributes the weights more equally, so we do it

            # Make the new large set
            sets[max_set] = [i for i in sets[max_set] if i!= max_el]
            push!(sets[max_set],min_el)
            set_weights[max_set] += weights[min_el] - weights[max_el]

            #make the new small set
            sets[min_set] = [i for i in sets[min_set] if i!= min_el]
            push!(sets[min_set],max_el)
            set_weights[min_set] += weights[max_el] - weights[min_el]
            # We swapped, so we start looking again from the extreme sets
            index_el = 1
            index_set = 1
        elseif index_el < step-1 #some sets have 'step' elements, but others have 'step-1' elements. We don't look at which set we have but just be safe
            # We didn't swap, so we want to look at different elements in the next iteration
            index_el+=1
        elseif index_el == step-1 && index_set < n-1 # try all pairs with the smallest before quitting
            # We didn't swap anything in this (large) set, so we try the next set
            index_set+=1
            index_el = 1
        else
            # nothing changed, so we can as well stop
            break
        end
    end
    return sets,set_weights,[weights[s] for s in sets]
end

struct ThreadingInfo
    j_order::Vector{Int} # the order of the clusters for threading
    jl_order::Vector{Tuple{Int,Int}} # the order of the (j,l) blocks for threading
end

function ThreadingInfo(sdp::ClusteredLowRankSDP,n=Threads.nthreads())
    # First we determine the blocksizes of Y. This determines the distribution over the threads
    subblocksizes = [[0 for l in eachindex(sdp.A[j])] for j in eachindex(sdp.A)]
    for j in eachindex(sdp.A)
        for l in eachindex(sdp.A[j])
            for r=1:size(sdp.A[j][l],1),s=1:size(sdp.A[j][l],2)
                for p in keys(sdp.A[j][l][r,s])
                    if subblocksizes[j][l] == 0
                        subblocksizes[j][l] = size(sdp.A[j][l][r,s][p],1)
                    elseif subblocksizes[j][l] != size(sdp.A[j][l][r,s][p],1)
                        s1 = subblocksizes[j][l]
                        s2 = size(sdp.A[j][l][r,s][p],1)
                        error("The subblocks (j,l,(r,s)) must have the same size for every r,s. $s1 != $s2")
                    end
                end
            end
        end
    end
    Y_blocksizes = [[size(sdp.A[j][l],1)*subblocksizes[j][l] for l in eachindex(sdp.A[j])] for j in eachindex(sdp.A)]

    # Now we find a good order for the (j,l) pairs:
    # the only relevant parameter here is the blocksize Y_blocksizes[j][l]
    # Computing the step length takes a cholesky ~n^3 . Eigenvalues should also take ~n^3
    jl_order = [(j,l) for j in eachindex(sdp.A) for l in eachindex(sdp.A[j])]
    jl_pair_weights = [Y_blocksizes[j][l]^3 for (j,l) in jl_order]
    setsjl, weightsjl, weight_distjl = distribute_weights_swapping(jl_pair_weights, n, nswaps=length(jl_pair_weights)^2)
    # to get these sets on the cores, we need to put the longer sets first
    sort!(setsjl, by=length, rev=true)
    jl_order = jl_order[vcat(setsjl...)]

    # Determine a good order for the clusters j:
    # We base the weights for the clusters on the Cholesky decomposition.
    # In principle we could instead use the time for the solving (P[j]^2*N). Not sure what is better, maybe P[j]^2*N/2 + 1/3*P[j]^3
    j_order_weights = [size(sdp.c[j],1)^3 for j in eachindex(sdp.c)] # [1//2 * size(sdp.c[j],1)^2*size(sdp.b,1) + 1//3 * size(sdp.c[j].1)^3 for j in eachindex(sdp.c)]
    setsj, weightsj, weight_distj = distribute_weights_swapping(j_order_weights, n, nswaps = length(j_order_weights)^2)
    sort!(setsj, by=length, rev=true)
    j_order = vcat(setsj...)
    return ThreadingInfo(j_order,jl_order)
end


text:
#gnum-division.jl
#division algorithms between unums, ubounds and gnums.

#=

doc"""
  `div!(::Unum, ::Gnum)` divides a unum BY a gnum.  the gnum will be altered as
  a result.  This function may also alter the unum by appending g-flags into the
  flags section.  Returns a reference to the gnum to allow chaining.
"""
function div!{ESS,FSS}(a::Unum{ESS,FSS}, b::Gnum{ESS,FSS})
  #make sure our ignore_sides parameters are set.
  clear_ignore_sides!(b)
  #go ahead and set "g-flags" on the value a.
  set_g_flags!(a)

  #ascertain the result sign.  Note that all division algorithms are
  #sign-agnostic and return positive-definite results.
  result_sign::UInt16 = @signof(a) $ @signof(b.lower)

  #perform hard override checks on the multiply.
  __check_div_hard_override!(a, b)

  #flip sides?

  #test sidedness and trampoline to the appropriate procedure.
  #is_onesided(b) ? div_onesided!(a, b) : div_twosided!(a, b)

  (result_sign != 0) && additive_inverse!(b)
end
=#
#=
doc"""
  `mul!(::Ubound, ::Gnum)` multiplies a ubound INTO a gnum.
  the gnum will be altered as a result.  This function may also alter the ubound
  by appending g-flags into the flags section of the constituent unums.  Returns
  a reference to the gnum to allow chaining.
"""
function mul!{ESS,FSS}(a::Ubound{ESS,FSS}, b::Gnum{ESS,FSS})
  throw(ArgumentError("not implemented yet."))
end

doc"""
  `mul!(a::Gnum, b::Gnum)` multiplies a gnum `a` INTO a gnum `b`.
  the *second* gnum will be altered as a result.  Following at&t assembly syntax,
  the result is on the right hand side (this is the equivalent of b \*= a ).
  Returns a reference to the gnum to allow chaining.
"""
function mul!{ESS,FSS}(a::Gnum{ESS,FSS}, b::Gnum{ESS,FSS})
  throw(ArgumentError("not implemented yet."))
end

################################################################################
## onesided/twosided multiplies.

doc"""
  `mul_onesided!(::Unum, ::Gnum)` performs a multiplication of a unum a into
  a unum b, where the starting gnum is one-sided.  This function does not check
  that the passed parameters match this description.
"""
function mul_onesided!{ESS,FSS}(a::Unum{ESS,FSS}, b::Gnum{ESS,FSS})
  #all procedures in onesided multiply presume that a is positive.

  if is_zero(b, LOWER_UNUM)
    ignore_side!(b, LOWER_UNUM)
  elseif is_inf(b, LOWER_UNUM)
    b.lower.flags &= ~UNUM_SIGN_MASK
    ignore_side!(b, LOWER_UNUM)
  elseif should_calculate(b, LOWER_UNUM) && is_unit(b.lower)
    copy_unum!(a, b.lower)
    #overwrite the sign; this will be re-xored later.
    b.lower.flags &= ~UNUM_SIGN_MASK
    ignore_side!(b, LOWER_UNUM)
  elseif is_mmr(a)
    should_calculate(b, LOWER_UNUM) && mmr_onesided_mult!(b)
  elseif is_sss(a)
    should_calculate(b, LOWER_UNUM) && sss_onesided_mult!(b)
  end

  #check for the gnum value being mmr or sss.
  if should_calculate(b, LOWER_UNUM)
    if is_mmr(b, LOWER_UNUM)
      copy_unum!(a, b.lower)
      mmr_onesided_mult!(b)
    elseif is_sss(b, LOWER_UNUM)
      copy_unum!(a, b.lower)
      sss_onesided_mult!(b)
    end
  end

  #we have exhausted all of the special cases, so just do the multiplication.
  should_calculate(b, LOWER_UNUM) && (__multiply!(a, b, LOWER_UNUM, UPPER_UNUM))
  #swap parity if a was negative.

  #return the gnum.
  b
end

doc"""
  `mul_twosided!(::Unum, ::Gnum)` performs a multiplication of a unum a into
  a unum b, where the starting gnum is one-sided.  This function does not check
  that the passed parameters match this description.
"""
function mul_twosided!{ESS,FSS}(a::Unum{ESS,FSS}, b::Gnum{ESS,FSS})
  throw(ArgumentError("twosided mults not implemented yet"))
end
=#
#=
doc"""
  `__transits_zero(::Gnum)` checks to see if a Gnum has a negative-definite left
  and a positive-definite right.  This is useful for division because such a gnum
  will throw an error.
"""
function __transits_zero{ESS,FSS}(b::Gnum{ESS,FSS})
  #onesided gnums cannot transit zero
  return is_twosided(b) & is_negative(b.lower) & is_positive(b.upper) & !(is_zero(b, LOWER_UNUM)) & !(is_zero(b, UPPER_UNUM))
end

################################################################################
## shortcut division, for special values.

doc"""
  `__check_div_hard_override!(::Unum, ::Gnum)` checks to see if the terms in
  the divisor or the dividend trigger a "hard override".  These values are,
  in cascading precedence:

  nan (dividend, divisor) -> outputs nan.
  value that transits zero (divisor) -> outputs nan.
"""
function __check_div_hard_override!{ESS,FSS}(a::Utype{ESS,FSS}, b::Gnum{ESS,FSS})
  ############################################
  # analysis of division with NaNs.
  # if our dividend is nan, then set the product to nan.
  is_nan(a) && (@scratch_this_operation!(b); return)
  is_nan(b) && (ignore_both_sides!(b); return)

  ############################################
  # check for values that transit zero
  __transits_zero(b) && (@scratch_this_operation!(b); return)
end

#checking for values in the dividend which trigger override.
function __check_div_dividend_override!{ESS,FSS}(a::Unum{ESS,FSS}, b::Gnum{ESS,FSS})
  ############################################
  # analysis of division with infs
  # inf / anything becomes inf, except inf, which becomes NaN.

  if (is_g_inf(a))
    #infinity can't multiply times just zero, either
    should_calculate(b, LOWER_UNUM) && is_inf(b, LOWER_UNUM) && (@scratch_this_operation!(b))
    should_calculate(b, UPPER_UNUM) && is_inf(b, UPPER_UNUM) && (@scratch_this_operation!(b))
    #set the lower side to infinity, then ignore further calculations, and make it onesided
    inf!(b, z16, LOWER_UNUM); ignore_side!(b, LOWER_UNUM); set_onesided!(b)
  end

  ############################################
  # analysis of division with zeros
  # zero divided by anything becomes zero.
  # n.b. this does not check twosided unums, because that has already been checked.
  if is_g_zero(a)
    zero!(b, LOWER_UNUM)
    set_onesided!(b)
    ignore_side!(b, LOWER_UNUM)
  end
end

function __check_div_divisor_override!{ESS,FSS}(a::Unum{ESS,FSS}, b::Gnum{ESS,FSS})
  #anything / inf is zero, but make sure the zero gets placed in the lower unum.
  if should_calculate(b, LOWER_UNUM) && is_inf(b, LOWER_UNUM)
    zero!(b, LOWER_UNUM)
    ignore_side!(b, LOWER_UNUM); set_onesided!(b)
  end
  if should_calculate(b, UPPER_UNUM) && is_inf(b, UPPER_UNUM)
    #move the value of the lower unum into the upper segment and then do the calculation.
    #note tthat using this implies that the upper unum is active, s o we don't need to
    #set this to onesided.
    copy_unum_with_gflags!(b.lower, b.upper)
    zero!(b, LOWER_UNUM)
    ignore_side!(b, LOWER_UNUM)
  end
end

doc"""
  `sss_onesided_mult!(::Gnum)` multiplies a onesided Gnum times an sss.
  the lower value of the Gnum contains the value to be multiplied.  The upper
  value will contain the result value.
"""
function sss_onesided_mult!{ESS,FSS}(b::Gnum{ESS,FSS})
  #if the magnitude is more than one, then we multiply the value of b times
  #smallsubnormal.
  if is_magnitude_less_than_one(b.lower)
    sss!(b, z16, LOWER_UNUM)
    ignore_side!(b, LOWER_UNUM)
  else
    #we're going to have a two-valued solution, so alias the inner and outer values.
    set_twosided!(b)
    #first, make the upper part small_exact.
    small_exact!(b.upper, z16)
    __outward_exact!(b.lower)
    #do the multiplication
    __multiply!(b.lower, b, UPPER_UNUM, DISCARD_SECONDARY)
    #check to see if we need to retrace.
    is_exact(b.upper) && __inward_ulp!(b.upper)
    #set the lower unum to be small subnormal.
    sss!(b, z16, LOWER_UNUM)
    #finished with calculations.
    ignore_both_sides!(b)
  end
end

doc""" `mmr_onesided_mult!(::Gnum)` multiplies a onesided Gnum times an mmr."""
function mmr_onesided_mult!{ESS,FSS}(b::Gnum{ESS,FSS})
  if is_magnitude_less_than_one(b.lower)
    #we're going to have a two-valued solution, so alias the inner and outer values.
    set_twosided!(b)
    #first, make the buffer big_exact.
    big_exact!(b.buffer, z16)
    #multiply the buffer times the value in the lower unum, overwriting it.
    __multiply!(b.buffer, b, LOWER_UNUM, DISCARD_SECONDARY)
    #check to see if we need to retrace.
    is_exact(b.lower) && make_ulp!(b.lower)
    #make the upper side mmr.
    mmr!(b, z16, UPPER_UNUM)
    ignore_both_sides!(b)
  else
    #then everything stays mmr.
    mmr!(b, z16, LOWER_UNUM)
    ignore_side!(b, LOWER_UNUM)
  end
end
=#
#=
################################################################################
## actual algorithmic multiplication

doc"""
  `Unums.__multiply!(::Unum, ::Gnum, ::Type{Val}, ::Type{Val})`
  performs the algorithmic multiplication of an unum into one side of the gnum.
  Because a unum muliplication times another unum must always be on the same
  side of zero, this procedure ignores sign, always returning a positive value.

  The *primary result* is the result of the exact multiplication between the two
  values; the *secondary result* is the result of multiplying the ulp values.

  The first value parameter contains the side that contains the multiplicand:
  this may be UPPER_GNUM, LOWER_GNUM, or BUFFER.

  The second value parameter contains a directive that determines what to do
  when the value is inexact.  If this is UPPER_GNUM, LOWER_GNUM, or BUFFER, then
  the secondary result is stored in the respective destination.  This destination
  cannot be identical to the first value parameter.

  You may also specify DISCARD_PRIMARY or DISCARD_SECONDARY as the directives,
  which will discard the primary or secondary results.
"""
@gen_code function __multiply!{ESS,FSS,side,directive}(a::Unum{ESS,FSS}, b::Gnum{ESS,FSS}, ::Type{Val{side}}, ::Type{Val{directive}})
  max_exp = max_exponent(ESS)
  sml_exp = min_exponent(ESS, FSS)
  min_exp = min_exponent(ESS)
  mesize = max_esize(ESS)
  mfsize = max_fsize(FSS)

  (side == directive) && throw(ArgumentError("multiply algorithm cannot send an ulp result to the same side as source"))

  @code quote
    a_exp::Int64 = decode_exp(a)
    b_exp::Int64 = decode_exp(b.$side)
    a_dev::UInt64 = is_exp_zero(a) * o64
    b_dev::UInt64 = is_exp_zero(b.$side) * o64

    #make the result positive
    b.scratchpad.flags &= ~UNUM_SIGN_MASK

    scratchpad_exp::Int64 = a_exp + b_exp + a_dev + b_dev

    #preliminary comparisons that will stop us from performing unnecessary steps.
    (scratchpad_exp > $max_exp + 1) && (@preserve_sflags b mmr!(b, z16, Val{side}); return)
    (scratchpad_exp < $sml_exp - 2) && (@preserve_sflags b sss!(b, z16, Val{side}); return)
  end

  if FSS < 6
    #fss < 6 can perform a multiplication within the existing 64-bit integer width.
    @code :(b.scratchpad.fraction = __chunk_mult_small(a.fraction, b.$side.fraction))
  elseif FSS == 6
    #fss = 6 requires access to the scratchpad array.
    @code :(nan!(b); return)
  else
    @code :(nan!(b); return)
  end

  @code quote

    #fraction multiplication:  where va is "virtual bit" that is 1 for normals or 0 for subnormals.
    #  va + a
    #  vb + b
    #  va*vb + vb * a + va * b + ab
    #
    carry::UInt64 = o64 - (a_dev | b_dev) #if either a or b are subnormal, the carry is zero.
    #note that adding the fractions back in has a cross-multiplied effect.
    (b_dev == 0) && (carry = __carried_add_frac!(carry, a, b.scratchpad))
    (a_dev == 0) && (carry = __carried_add_frac!(carry, b.$side, b.scratchpad))

    #next, do carry analysis. carry can be 0, 1, 2, or 3 at the most.
    if (carry == 0)
      #shift over the fraction as far as we need to.
      shift::Int64 = leading_zeros(b.scratchpad.fraction) + 1
      #check if we're zero.
      is_frac_zero(b.scratchpad) && (sss!(b, z16, Val{side}); return)
      __leftshift_frac!(b.scratchpad, shift)
      #adjust the exponent.
      scratchpad_exp -= shift
      (scratchpad_exp < $sml_exp) && (sss!(b, z16, Val{side}); return)
      if (scratchpad_exp < $min_exp)
        #calculate the leftshift.
        shift = $min_exp - scratchpad_exp
        __rightshift_frac_with_underflow_check!(b.scratchpad, shift)
        set_frac_bit!(b.scratchpad, shift)
        b.scratchpad.esize = $mesize
        b.scratchpad.exponent = 0
      else
        (b.scratchpad.esize, b.scratchpad.exponent) = encode_exp(scratchpad_exp)
        carry = 1
      end
    elseif (carry == 1)
      if scratchpad_exp > $max_exp
        @preserve_sflags b (mmr!(b, z16, Val{side}); return)
      end
      (b.scratchpad.esize, b.scratchpad.exponent) = encode_exp(scratchpad_exp)
    else
      (scratchpad_exp += 1)
      if scratchpad_exp > $max_exp
        @preserve_sflags b (mmr!(b, z16, Val{side}); return)
      else
        #shift things over by one since we went up in size.
        __rightshift_frac_with_underflow_check!(b.scratchpad, 1)
        #carry from the carry over into the fraction.
        (carry == 3) && set_frac_top!(b.scratchpad)
        #re-encode the exponent.
        (b.scratchpad.esize, b.scratchpad.exponent) = encode_exp(scratchpad_exp)
      end
      #set the carry to one
      carry = 1
    end
  end

  #if we're instructed to discard the secondary value, then simply copy over
  #the values, and we're done.
  (directive == :discard_secondary) && (@code :(copy_unum_with_gflags!(b.scratchpad, b.$side); return); @goto __mult_end)
  #assign the destination for the secondary varibale.
  sdest = (directive == :discard_primary) ? side : directive
  #this code ge
  retrace_code1 = (directive == :discard_primary) ? :(nothing) : :(carry = __tandem_copy_add_frac!(carry, b.$side, b.scratchpad, a.fsize))
  retrace_code2 = (directive == :discard_primary) ? :(nothing) : :(copy_unum!(b.scratchpad, b.$side))
  needs_twosided = (sdest == :upper) ? (:(set_twosided!(b))) : :(nothing)
  #if we're writing to the lower side, then be prepared to collapse an mmr into a onesided.
  mmr_collapse = (side == :lower) && !(directive == :discard_primary) ? :(set_onesided!(b)) : :(nothing)

  @code quote
    #if both sides were exact, the result is fairly simple.  Just copy the scratchpad results
    #on over to the destination side, we don't have to do anything further.

    is_exact(a) && is_exact(b.$side) && (copy_unum_with_gflags!(b.scratchpad, b.$side); return)
    $needs_twosided

    #the scratchpad needs to be exact, but as a temporary shim we need to make
    #it an ulp to assess if the result is mmr.
    make_ulp!(b.scratchpad)
    is_mmr(b.scratchpad) && ($mmr_collapse; mmr!(b, z16, Val{side}); return)
    make_exact!(b.scratchpad)

#    println("a: ", bits(a))
#    println("1) scr: ", bits(b.scratchpad), " carry: ", carry)

    #from here on out we engage the heuristic algorithm.  The 'direct' algorithm
    #would be to do two multiplications, but it is possible to save on that
    #calculation.  The algorithm is as follows:
    # we're going to do the multiplication X1 * X2
    # X1 = (2^P1)(F1) (+) (2^(P1-U1))
    # X2 = (2^P2)(F2) (+) (2^(P2-U2))
    # X1 * X2 = 2^(P1P2)(F1F2 (+) (2^-U1)F2 + (2^-U2)F1 + 2^(-U1-U2))
    # where x (+) a ≡ x + (0, a)

    #copy the scratchpad back, and add into the scratchpad b.$side.fraction >> a.fsize
    (is_ulp(a)) ? $retrace_code1 : $retrace_code2

#    println("2) scr: ", bits(b.scratchpad), " carry: ", carry)
    #add into the scratchpad a.fraction >> b.$side.fsize
    (is_ulp(b.$side)) && (carry = __add_frac_with_shift!(carry, a, b.scratchpad, b.$side.fsize))
    #check to make sure the far value isn't exact; move it inward if so.

    #tentatively make the fsize the maximum
    b.$side.fsize = $mfsize

    #resolve the scratchpad's carry
    __carry_resolve!(carry, b.scratchpad, Val{true})

    is_exact(b.scratchpad) && __inward_ulp!(b.scratchpad)

    #move the scratchpad back to the destination
    copy_unum!(b.scratchpad, b.$sdest)
  end

  @label __mult_end
end
=#

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
function visualizeDensityMesh!(vc,
                               fgl::FactorGraph,
                               lbl::Symbol;
                               levels=3,
                               meshid::Int=2  )

  pl1 = marginal(getVertKDE(fgl,lbl),[1;2;3])

  gg = (x, a=0.0) -> evaluateDualTree(pl1, [x[1] x[2] x[3]]')[1]-a  #([x[1];x[2];x[3]]')'

  x = getKDEMax(pl1)
  maxval = gg(x)

  vv = getKDERange(pl1)
  lower_bound = Vec(vec(vv[:,1])...)
  upper_bound = Vec(vec(vv[:,2])...)

  levels = linspace(0.0,maxval,levels+2)

  # MD = []
  for val in levels[2:(end-1)]
    meshdata = GeometryData(contour_mesh(x -> gg(x,val), lower_bound, upper_bound))
    meshdata.color = RGBA( val/(1.5*maxval),0.0,1.0,val/(1.5*maxval))
    # push!(MD, meshdata)
    setgeometry!(vc[:meshes][lbl][Symbol("lev$(val)")], meshdata)
  end
  # mptr = Any(MD)
  # vc.meshes[lbl] = mptr
  # Visualizer(mptr, meshid) # meshdata
  nothing
end

function drawgt!(vc,
                 sym::Symbol,
                 gtval::Tuple{Symbol, Vector{Float64}};
                 session::AbstractString="NA"  )
  #
  if gtval[1] == :XYZ
    drawpoint!(vc, sym, tf=Translation(gtval[2][1],gtval[2][2],gtval[2][3]),
          session=session,
          color=RGBA(1.0,0,0,0.5),
          collection=:gt_landm  )
  elseif gtval[1] == :XYZqWXYZ
    drawpose!(vc, sym,
          tf = Translation(gtval[2][1],gtval[2][2],gtval[2][3]) ∘
               LinearMap(CoordinateTransformations.Quat(gtval[2][4],gtval[2][5],gtval[2][6],gtval[2][7])),
          session=session,
          collection=:gt_poses  )
  else
    warn("unknown ground truth drawing type $(gtval[1])")
  end

  nothing
end

# TODO -- maybe we need RemoteFactorGraph type
function visualizeallposes!(vc,
                            fgl::FactorGraph;
                            drawlandms::Bool=true,
                            drawtype::Symbol=:max,
                            gt::Dict{Symbol, Tuple{Symbol,Vector{Float64}}}=Dict{Symbol, Tuple{Symbol,Vector{Float64}}}(),
                            api::DataLayerAPI=localapi  )
  #
  session = fgl.sessionname
  topoint = gettopoint(drawtype)

  dotwo = false
  dothree = false
  po,ll = ls(fgl)
  if length(po) > 0
    sym = po[1]
    X = getVal(fgl, sym, api=api )
    dotwo, dothree = getdotwothree(sym, X)
  end

  # TODO -- move calls higher in abstraction to be more similar to drawdbdirector()
  for p in po
    vert = getVert(fgl, p, api=api )
    drawpose!(vc, vert, topoint, dotwo, dothree, session=session)
    if haskey(gt, p)
      drawgt!(vc, p, gt[p], session=session)
    end
  end
  if drawlandms
    for l in ll
      den = getVertKDE(fgl, l, api=api)
      pointval = topoint(den)
      drawpoint!(vc, l, tf=Translation(pointval[1:3]...), session=session)
      if haskey(gt, l)
        drawgt!(vc, l, gt[l], session=session)
      end
    end
  end

  nothing
end

function drawposepoints!(vis,
                         vert::Graphs.ExVertex;
                         session::AbstractString="NA"  )
  #
  vsym = Symbol(vert.label)
  X = getVal(vert)

  dotwo, dothree = getdotwothree(vsym, X)
  makefromX = (X::Array{Float64,2}, i::Int) -> X[1:3,i]
  if dotwo
    makefromX = (X::Array{Float64,2}, i::Int) -> Float64[X[1:2,i];0.0]
  end

  XX = Vector{Vector{Float64}}()
  for i in 1:size(X,2)
    push!(XX, makefromX(X,i))
  end
  pointcloud = PointCloud(XX)
  if string(vsym)[1] == 'l'
    pointcloud.channels[:rgb] = [RGB(1.0, 1.0, 0) for i in 1:length(XX)]
  elseif string(vsym)[1] == 'x'
    pointcloud.channels[:rgb] = [colorwheel(vert.index) for i in 1:length(XX)]
  end
  setgeometry!(vis[Symbol(session)][:posepts][vsym], pointcloud)
  nothing
end

function drawposepoints!(vis,
      fgl::FactorGraph,
      sym::Symbol;
      session::AbstractString="NA",
      api::DataLayerAPI=dlapi  )
  #
  vert = getVert(fgl, sym, api=api)
  drawposepoints!(vis, vert, session=session, api=localapi) # definitely use localapi
  nothing
end

function deletemeshes!(vc)
  delete!(vc[:meshes])
end

"""
    $(SIGNATURES)

Draw a line segment between to vertices.
"""
function drawLineBetween!(vis,
                          session::AbstractString,
                          fr::Graphs.ExVertex,
                          to::Graphs.ExVertex;
                          scale=0.01,
                          name::Symbol=:edges,
                          subname::Union{Nothing,Symbol}=nothing,
                          color=RGBA(0,1.0,0,0.5)   )
  #
  dotwo, dothree = getdotwothree(Symbol(fr.label), getVal(fr))

  xipt = zeros(3); xjpt = zeros(3);
  if dothree
    xi = marginal(getVertKDE( fr ),[1;2;3] )
    xj = marginal(getVertKDE( to ),[1;2;3] )
    xipt[1:3] = getKDEMax(xi)
    xjpt[1:3] = getKDEMax(xj)
  elseif dotwo
    xi = marginal(getVertKDE( fr ),[1;2] )
    xj = marginal(getVertKDE( to ),[1;2] )
    xipt[1:2] = getKDEMax(xi)
    xjpt[1:2] = getKDEMax(xj)
  end

  lbl = Symbol(string(fr.label,to.label))
  place = vis[Symbol(session)][name][lbl]
  if subname != nothing
    place = vis[Symbol(session)][name][subname][lbl]
  end
  drawLine!(place, xipt, xjpt, color=color, scale=scale )
  nothing
end

"""
    $(SIGNATURES)

Draw a line segment between to nodes in the factor graph.
"""
function drawLineBetween!(vis,
                          fgl::FactorGraph,
                          fr::Symbol,
                          to::Symbol;
                          scale=0.01,
                          name::Symbol=:edges,
                          subname::Union{Nothing,Symbol}=nothing,
                          color=RGBA(0,1.0,0,0.5),
                          api::DataLayerAPI=dlapi  )
  #
  v1 = getVert(fgl, fr, api=api)
  v2 = getVert(fgl, to, api=api)

  drawLineBetween!(vis,fgl.sessionname, v1,v2,scale=scale,name=name,subname=subname,color=color   )
end

"""
    $(SIGNATURES)

Assume odometry chain and draw edges between subsequent poses. Use keyword arguments to change colors, etc.
"""
function drawAllOdometryEdges!(vis,
                               fgl::FactorGraph;
                               scale=0.01,
                               name::Symbol=:edges,
                               color=RGBA(0,1.0,0,0.5),
                               api::DataLayerAPI=dlapi  )
  #
  xx, ll = ls(fgl)

  for i in 1:(length(xx)-1)
    drawLineBetween!(vis, fgl, xx[i],xx[i+1], api=api , color=color, scale=scale, name=name )
  end

  nothing
end

function drawAllBinaryFactorEdges!(vis,
                                   fgl::FactorGraph;
                                   scale=0.01,
                                   api::DataLayerAPI=dlapi  )
  #
  sloth = findAllBinaryFactors(fgl, api=api)

  for (teeth, toe) in sloth
    color = pointToColor(toe[3])
    drawLineBetween!(vis, fgl, toe[1], toe[2], subname=toe[3], scale=scale, color=color)
  end
  nothing
end

"""
    $(SIGNATURES)

Update all triads listed in poseswithdepth[Symbol(vert.label)] with wTb. Prevents cycles in
remote tree viewer of Visualizer.
"""
function updateparallelposes!(vis,
                              vert::Graphs.ExVertex,
                              poseswithdepth::Dict;
                              wTb::CoordinateTransformations.AbstractAffineMap=
                                    Translation(0,0,0.0) ∘ LinearMap(
                                    CoordinateTransformations.Quat(1.0, 0, 0, 0))    )
  #

  if haskey(poseswithdepth, Symbol(vert.label))
    for cp in poseswithdepth[Symbol(vert.label)]
      settransform!(cp, wTb)
    end
  end
  nothing
end

function drawdbdirector(;addrdict::NothingUnion{Dict{AbstractString, AbstractString}}=nothing)
  # Uncomment out for command line operation
  cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict,drawdepth=true, drawedges=true)
  session = addrdict["session"]

  poseswithdepth = Dict()
  # poseswithdepth[:x1] = 0 # skip this pose -- there is no big data before ICRA

  vis = startDefaultVisualization()
  sleep(1.0)

  param = Dict()
  try
    param = robotsetup(cloudGraph, session)
  catch
    warn("No robot parameters found for the session, continuing with basic visualization only")
  end

  drawloop = Bool[true]
  println("Starting draw loop...")
  while drawloop[1]
    drawdbsession(vis, cloudGraph, addrdict, param, poseswithdepth) #,  db[collection]
    println(".")
    sleep(0.5)
  end

  println("Finishing askdrawdirectordb.")
  nothing
end

#syntax support lost in 0.5, but see conversation at
# function (dmdl::DrawObject)(vc::VisualizationContainer)
# issue  https://github.com/JuliaLang/julia/issues/14919


text:
module PlotTests
using ControlSystems, Plots
using VisualRegressionTests, ControlExamplePlots
using Test
gr()
default(show=false)

@testset "test_plots" begin
funcs, refs, eps = getexamples()
# Make it easier to pass tests on different systems
# Set to a factor 2 of common errors
eps = [0.15, 0.015, 0.1, 0.01, 0.01, 0.02, 0.01, 0.15, 0.15, 0.01, 0.01]
res = genplots(funcs, refs, eps=eps, popup=false)

##Explicit enumeration for simpler debugging

#"bode.png"
@test  res[1] |> success
#"nyquist.png"
@test  res[2] |> success
#"sigma.png"
@test  res[3] |> success
#"nichols.png"
@test  res[4] |> success
#"step.png"
@test  res[5] |> success
#"impulse.png"
@test  res[6] |> success
#"lsim.png"
@test  res[7] |> success
#"margin.png"
@test  res[8] |> success
#"gangoffour.png"
@test  res[9] |> success
#"pzmap.png"
@test  res[10] |> success
#"rlocus.png"
@test  res[11] |> success
end

end


text:
using SDPSymmetryReduction
using Documenter
using Literate
using Test

ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "100"

DocMeta.setdocmeta!(SDPSymmetryReduction, :DocTestSetup, :(using SDPSymmetryReduction); recursive=true)

## Use Literate.jl to generate examples (functions modified from https://github.com/jump-dev/JuMP.jl/blob/master/docs/make.jl)

function _file_list(full_dir, relative_dir, extension)
    return map(
        file -> joinpath(relative_dir, file),
        filter(file -> endswith(file, extension), sort(readdir(full_dir))),
    )
end

"""
    _include_sandbox(filename)
Include the `filename` in a temporary module that acts as a sandbox. (Ensuring
no constants or functions leak into other files.)
"""
function _include_sandbox(filename)
    mod = @eval module $(gensym()) end
    return Base.include(mod, filename)
end

function _literate_directory(dir)
    rm.(_file_list(dir, dir, ".md"))
    for filename in _file_list(dir, dir, ".jl")
        # `include` the file to test it before `#src` lines are removed. It is
        # in a testset to isolate local variables between files.
        Test.@testset "$(filename)" begin
            _include_sandbox(filename)
        end
        Literate.markdown(
            filename,
            dir;
            documenter = true,
            credit = true,
        )
    end
    return nothing
end

_literate_directory.(joinpath(@__DIR__, "src", "examples"))

## Generate docs

makedocs(;
    modules=[SDPSymmetryReduction],
    authors="Daniel Brosch <daniel.brosch@outlook.com> and contributors",
    repo="https://github.com/DanielBrosch/SDPSymmetryReduction.jl/blob/{commit}{path}#{line}",
    sitename="SDPSymmetryReduction.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DanielBrosch.com/SDPSymmetryReduction.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => [
            "examples/ErdosRenyiThetaFunction.md",
            "examples/QuadraticAssignmentProblems.md",
            "examples/ReduceAndSolveJuMP.md"
        ]
    ],
)
##
deploydocs(;
    repo="github.com/DanielBrosch/SDPSymmetryReduction.jl",
    devbranch="main",
)


text:
"""
    gen_prn()

Generate a pseudorandom number.

# Notes

Uses a linear congruential generator (LCG) with [POSIX parameters](https://en.wikipedia.org/wiki/Linear_congruential_generator#Parameters_in_common_use):

``X_n = 25214903917 X_{n-1} + 11 \\quad mod \\quad 2^{48}``

# Examples

```julia-repl
julia> gen_prn()
 156750217634815

julia> gen_prn()
 63914890472862
```
"""
function gen_prn()
    seed = get_seed()
	prn = mod(A * seed + C, MOD)
    set_seed(prn)
    return prn
end

"""
    get_std_uniform(shape=1; seed=nothing)

Generate a `shape` element array of random variables from a standard Uniform(0,1) distribution. Optionally you can set a specific seed.

# Examples

```julia-repl
julia> get_std_uniform()
1-element Vector{Float64}:
 0.42443098343863284

julia> get_std_uniform(seed=43)
1-element Vector{Float64}:
 0.09636209187468836

julia> get_std_uniform(5)
5-element Vector{Float64}:
 0.6584669595802204
 0.33437978955868886
 0.509019330923099
 0.12156905126458639
 0.917393216014684

```
"""
function get_std_uniform(shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    seed_setter(seed)
    U = zeros(shape)  # preallocate array
    U .= gen_prn.()  # vectorize assignment for efficiency
    U = U./MOD
	return U
end

"""
    uniform_rng(a, b, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Uniform(`a`, `b`) distribution. Optionally you can set a specific seed.

# Notes

The Uniform distribution is given by:

``f(x, a, b) = \\frac{1}{b-a} \\quad \\text{for} \\quad a ≤ x ≤ b``

# Examples

```julia-repl
julia> uniform_rng(1, 6)
1-element Vector{Float64}:
 2.638331960912094

julia> uniform_rng(1, 6, seed=42)
1-element Vector{Float64}:
 2.6333962626438314

julia> uniform_rng(0, 1, (4,4))
4×4 Matrix{Float64}:
 0.640603   0.757195  0.325722  0.645452
 0.955188   0.155203  0.953206  0.0046541
 0.0923526  0.490721  0.451705  0.516445
 0.661619   0.527063  0.212847  0.832298

```

# References

Walk, C. Handbook on statistical distributions for experimentalists. 2007.
"""
function uniform_rng(a::Real=0, b::Real=1, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    U = get_std_uniform(shape, seed=seed)
    X = a .+ (b-a) .* U
    return X
end


text:
include("08_two_site_operator_custom.jl")

println("
#######################################################
# Tutorial 9                                          #
#                                                     #
# 2-site energy minimization via gradient             #
# optimization                                        #
#######################################################
")

# Used for computing gradients "automatically".
# Julia's source-to-source reverse mode AD library.
using Zygote

# Energy function to minimize with gradient descent.
# Depends implicitly on the Hamiltonian.
function E(ψ)
  ψHψ = inner(ψ', H, ψ)

  # Alternative:
  # Hψ = apply(H, ψ)
  # ψHψ = inner(ψ, Hψ)

  ψψ = inner(ψ, ψ)
  return ψHψ / ψψ
end

# Extremely simple gradient descent.
# Better to use a library like `OptimKit.jl` in general.
function minimize(f, ∂f, x; nsteps, γ)
  for n in 1:nsteps
    println("n = ", n, ", f_x = ", f(x))

    # Gradient descent step with step size `γ`
    x = x - γ * ∂f(x)
  end
  return x
end

# Use the cat state as a starting state
ψ₀ = (Zp1 * Zp2 + Zm1 * Zm2) / √2

# Minimize the energy!
∂E(ψ) = gradient(E, ψ)[1]
ψ = minimize(E, ∂E, ψ₀; nsteps=10, γ=0.1)

@show E(ψ₀), norm(∂E(ψ₀))
@show E(ψ), norm(∂E(ψ))


text:
# This file is a part of TextClassification.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt
export MicroTC_Config, MicroTC_ConfigSpace, TextConfigSpace

struct MicroTC_Config{T_, C_}
    textconfig::TextConfig
    textmodel::T_
    cls::C_
end

function MicroTC_Config(;
        textconfig=TextConfig(),
        textmodel=EntModelConfig(),
        cls=LiblinearConfig(1.0, 0.1)
    )

    MicroTC_Config(textconfig, textmodel, cls)
end

function Base.show(io::IO, config::MicroTC_Config)
    print(io, "{MicroTC_Config ")
    show(io, config.textconfig)
    print(io, " ")
    show(io, config.textmodel)
    print(io, " ")
    show(io, typeof(config.cls))
    print(io, " }")
end

struct MicroTC_ConfigSpace <: AbstractSolutionSpace
    textconfig::TextConfigSpace
    textmodel
    cls
end

Base.eltype(::MicroTC_ConfigSpace) = MicroTC_Config
Base.hash(c::MicroTC_Config) = hash(repr(c))
Base.isequal(a::MicroTC_Config, b::MicroTC_Config) = repr(a) == repr(b)

function MicroTC_ConfigSpace(;
        kernel=[k_(CosineDistance()) for k_ in [DirectKernel, ReluKernel]],
        centerselection=[
            TextCentroidSelection(),
            #MedoidSelection(dist=CosineDistance()),
            #KnnCentroidSelection(sel1=CentroidSelection(), sel2=CentroidSelection(), dist=CosineDistance())
        ],
        ncenters=[3, 7, 11],
        textmodel=[EntModelConfigSpace(), VectorModelConfigSpace()],
        cls=[
            KncConfigSpace(centerselection=centerselection, kernel=kernel),
            # KncPerClassConfigSpace{0.3}(centerselection=centerselection, kernel=kernel, ncenters=ncenters),
            LiblinearConfigSpace(),
            KnnClassifierConfigSpace(k=1:2:11, keeptop=0.5:0.1:1.0)
        ],
        textconfig::TextConfigSpace = TextConfigSpace()
    )

    MicroTC_ConfigSpace(textconfig, textmodel, cls)
end

function Base.rand(space::MicroTC_ConfigSpace)
    T = space.textmodel isa AbstractArray ? rand(space.textmodel) : space.textmodel
    C = space.cls isa AbstractArray ? rand(space.cls) : space.cls

    MicroTC_Config(;
        textconfig=rand(space.textconfig),
        textmodel=rand(T),
        cls=rand(C)
    )
end

function mutate(space::MicroTC_ConfigSpace, c::MicroTC_Config, iter)
    textconfig = mutate(space.textconfig, c.textconfig, iter)
    textmodel = mutate(space.textmodel, c.textmodel, iter)
    cls = mutate(space.cls, c.cls, iter)
    MicroTC_Config(textconfig=textconfig, textmodel=textmodel, cls=cls)
end

function combine_select(a::MicroTC_Config, L::AbstractVector)
    textconfig = combine_select(a.textconfig, [b.first.textconfig => b.second for b in L])
    textmodel = combine_select(a.textmodel, [b.first.textmodel => b.second for b in L])
    cls = combine_select(a.cls, [b.first.cls => b.second for b in L])
    MicroTC_Config(textconfig=textconfig, textmodel=textmodel, cls=cls)
end


text:
function update_mask(line, one_mask, x_mask)
    one_mask = 0
    x_mask = 0
    for idx in 1:36
        if line[7 + idx] == '1'
            one_mask |= 1 << (36 - idx)
        elseif line[7 + idx] == 'X'
            x_mask |= 1 << (36 - idx)
        end
    end
    return (one_mask, x_mask)
end

function store(memory::Dict{Int,Int}, base_addr::Int, value::Int, x_mask::Int)
    addr = base_addr & ~x_mask
    memory[addr] = value
    max_addr = base_addr | x_mask
    while addr < max_addr
        for idx in 36:-1:1
            c = (1 << (36 - idx))
            if x_mask & c != 0
                if addr & c == 0
                    zero_mask = (c - 1) & addr  & x_mask
                    addr = (addr | c) & ~zero_mask
                    break
                end
            end
        end
        memory[addr] = value
    end
end

function run(s)
    # Your code here
    one_mask::Int = 0
    x_mask::Int = 0
    memory = Dict{Int,Int}()
    addr::Int = 0
    value::Int = 0
    for line in readlines(IOBuffer(s))
        if startswith(line, "mask")
            (one_mask, x_mask) = update_mask(line, one_mask, x_mask)
        else
            parts = split(line, " = ")
            addr = parse(Int, SubString(parts[1], 5, lastindex(parts[1]) - 1))
            addr |= one_mask
            value = parse(Int, parts[2])
            store(memory, addr, value, x_mask)
        end
    end
    return sum(values(memory))
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
# Demo of multivariate normal HMM using Viterbi algorithm and Baum-Welch.
include("models.jl")

module demo_mvn
using PyPlot
using Distributions
using MVNdiagHMM

# settings
m = 10
n = 1000
d = 5
tolerance = 1e-3

srand(1)

# parameters
p = round(rand(Dirichlet(ones(m))),2)
p[m] = 1-sum(p[1:m-1])
log_pi = log(p)
T = zeros(m,m)
alpha = 3.0
for r = 1:m
    T[r,:] = floor(rand(Dirichlet((alpha/m)*ones(m)))*100)/100
    T[r,m] = 1-sum(T[r,1:m-1])
end
log_T = log(T)
mu = round(4*randn(m,d),2)
sig = round(sqrt(rand(InverseGamma(1,2),(m,d))),2)
phi = [MvNormal(vec(mu[s,:]),vec(sig[s,:])) for s=1:m]

# simulate data
x,z0 = MVNdiagHMM.generate(n,log_pi,log_T,phi)

# compute optimal z
z = MVNdiagHMM.viterbi(x,log_pi,log_T,phi)

# compute naive z using each obs separately
z_naive = [indmax([MVNdiagHMM.log_q(x[i],phi[s]) for s=1:m]) for i=1:n]

# write to file
#writedlm("pi_T.txt",[p'; zeros(m)'; T])
#writedlm("mu_sig.txt",[mu; zeros(d)'; sig])
#writedlm("x.txt",[x[i][j] for i=1:n, j=1:d])
#writedlm("z_true.txt",z0)

# estimate parameters
log_pi_est,log_T_est,phi_est,log_m_est = MVNdiagHMM.estimate(x,m,tolerance)

# compute log marginal likelihood under the true parameters for comparison
G,log_m = MVNdiagHMM.forward(x,log_pi,log_T,phi)

# display results
println("\nViterbi percent correct:")
println(mean(z.==z0))
println("\nNaive percent correct:")
println(mean(z_naive.==z0))
println("\nLog marginal likelihood under true parameters:")
println(log_m)
println("\nLog marginal likelihood under estimated parameters:")
println(log_m_est)

# plots
x1 = Float64[xi[1] for xi in x]
m10 = Float64[mu[s,1] for s in z0]
m1p = Float64[mu[s,1] for s in z]
figure(1); clf(); hold(true)
plot(1:n,x1,"k.")
plot(1:n,m10,"b-")

figure(2); clf(); hold(true)
plot(1:n,z0,"b-",linewidth=2)
plot(1:n,z,"r--",linewidth=2)
plot(1:n,z_naive,"g:",linewidth=2)
xlim(0,100)

end # module

max_stars_repo_name: jwmi/HMM

max_issues_repo_name: jwmi/HMM

max_forks_repo_name: jwmi/HMM


text:
using Test
using AdvancedHMC
using AdvancedHMC: DualValue, PhasePoint
using CUDA

@testset "AdvancedHMC GPU" begin
    include("common.jl")

    n_chains = 1000
    n_samples = 1000
    dim = 5

    T = Float32
    m, s, θ₀ = zeros(T, dim), ones(T, dim), rand(T, dim, n_chains)
    m, s, θ₀ = CuArray(m), CuArray(s), CuArray(θ₀)

    target = Gaussian(m, s)
    metric = UnitEuclideanMetric(T, size(θ₀))
    ℓπ, ∇ℓπ = get_ℓπ(target), get_∇ℓπ(target)
    hamiltonian = Hamiltonian(metric, ℓπ, ∇ℓπ)
    integrator = Leapfrog(one(T) / 5)
    proposal = HMCKernel(Trajectory{EndPointTS}(integrator, FixedNSteps(5)))

    samples, stats = sample(hamiltonian, proposal, θ₀, n_samples)
end

@testset "PhasePoint GPU" begin
    for T in [Float32, Float64]
        init_z1() = PhasePoint(
            CuArray([T(NaN) T(NaN)]),
            CuArray([T(NaN) T(NaN)]),
            DualValue(CuArray(zeros(T, 2)), CuArray(zeros(T, 1, 2))),
            DualValue(CuArray(zeros(T, 2)), CuArray(zeros(T, 1, 2)))
        )
        init_z2() = PhasePoint(
            CuArray([T(Inf) T(Inf)]),
            CuArray([T(Inf) T(Inf)]),
            DualValue(CuArray(zeros(T, 2)), CuArray(zeros(T, 1, 2))),
            DualValue(CuArray(zeros(T, 2)), CuArray(zeros(T, 1, 2)))
        )

        @test_logs (:warn, "The current proposal will be rejected due to numerical error(s).") init_z1()
        @test_logs (:warn, "The current proposal will be rejected due to numerical error(s).") init_z2()

        z1 = init_z1()
        z2 = init_z2()

        @test z1.ℓπ.value == z2.ℓπ.value
        @test z1.ℓκ.value == z2.ℓκ.value
    end
end


text:
function detect_leave(x,tipo)
    if tipo == :has_left
        v = x .!= lag(x)
        v[1] = false
    elseif tipo == :will_leave
        v = lead(x) .!= x
    end
    return v
end

function count_bout(rewards,side)
    sidechanges = detect_leave(side,:has_left)
    v = Vector{Int64}(undef,length(rewards))
    c = 1
    for i =1:length(rewards)
        if  rewards[i] || sidechanges[i]
            c = c+1
        end
        v[i] = c
    end
    return v
end

function process_bouts(df::DataFrames.AbstractDataFrame)
    #df[!,:sidechange] = df.Side .!= lag(df.Side)
    df[:,:Leave] = lead(df.Side) .!= df.Side
    #df.sidechange[1] = false
    df[:,:Bout] = count_bout(df.Reward,df.Side)
    #df[!,:Omissions_plus_one].=0
    dayly_vars_list = [:Session, :MouseID, :Day, :Daily_Session, :Box,:Protocol,:Stim_Day];
    bout_table = by(df, :Bout) do dd
        dt = DataFrame(
        Omissions_plus_one = size(dd,1),
        Leave = dd[end,:Leave],
        Start_Reward = dd[1,:Reward],
        Side = dd[1,:Side],
        Trial_duration = (dd[end,:PokeOut]-dd[1,:PokeIn]),
        Pre_Interpoke = size(dd,1) > 1 ? maximum(skipmissing(dd[!,:Pre_Interpoke])) : missing,
        Post_Interpoke = size(dd,1) > 1 ? maximum(skipmissing(dd[!,:Post_Interpoke])) : missing,
        Wall = dd[1,:Wall],
        Correct_start = dd[1,:Correct],
        Correct_leave = !dd[end,:Correct],
        Block = dd[1,:Block],
        Streak_within_Block = dd[1,:Streak_within_Block],
        Streak = dd[1,:Streak],
        ReverseStreak = dd[1,:ReverseStreak],
        Protocol_Day = dd[1,:Protocol_Day],
        Exp_Day = dd[1,:Exp_Day]
        )
        for s in dayly_vars_list
            if s in names(df)
                dt[!,s] .= df[1, s]
            end
        end
        return dt
    end
    return bout_table
end


text:
using Test, ScientificTypes, ScientificTypesBase, Random
using Tables, CategoricalArrays, DataFrames
using ColorTypes
using Dates

# using CSV # dropped until julia release new LTS as issue for 1.0
import Distributions

const Dist = Distributions
const Arr  = AbstractArray
const CArr = CategoricalArray
const Cat  = CategoricalValue
const Vec  = AbstractVector

include("convention/coerce.jl")
include("convention/scitype.jl")

include("autotype.jl")
include("coerce.jl")
include("schema.jl")


text:
import Base64
export dcc_send_file, dcc_send_string, dcc_send_bytes
"""
    dbc_send_file(path::AbstractString, filename = nothing; type = nothing)

Convert a file into the format expected by the Download component.

# Arguments
- `path` - path to the file to be sent
- `filename` - name of the file, if not provided the original filename is used
- `type` - type of the file (optional, passed to Blob in the javascript layer)
"""
function dcc_send_file(path, filename = nothing; type = nothing)
    filename = isnothing(filename) ? basename(path) : filename
    return dcc_send_bytes(read(path), filename, type = type)
end

"""
    dcc_send_bytes(src::AbstractVector{UInt8}, filename; type = nothing)
    dcc_send_bytes(writer::Function, data, filename; type = nothing)

Convert vector of bytes into the format expected by the Download component.
`writer` function must have signature `(io::IO, data)`

# Examples

Sending binary content
```julia
file_data = read("path/to/file")
callback!(app, Output("download", "data"), Input("download-btn", "n_clicks"), prevent_initial_call = true) do n_clicks
    return dcc_send_bytes(file_data, "filename.fl")
end
```

Sending `DataFrame` in `Arrow` format
```julia
using DataFrames, Arrow
...
df = DataFrame(...)
callback!(app, Output("download", "data"), Input("download-btn", "n_clicks"), prevent_initial_call = true) do n_clicks
    return dcc_send_bytes(Arrow.write, df, "df.arr")
end
```
"""
function dcc_send_bytes(src::AbstractVector{UInt8}, filename; type = nothing)

    return Dict(
        :content => Base64.base64encode(src),
        :filename => filename,
        :type => type,
        :base64 => true
    )
end

function dcc_send_bytes(writer::Function, data, filename; type = nothing)
    io = IOBuffer()
    writer(io, data)
    return dcc_send_bytes(take!(io), filename, type = type)
end

"""
    dcc_send_data(src::AbstractString, filename; type = nothing)
    dcc_send_data(writer::Function, data, filename; type = nothing)

Convert string into the format expected by the Download component.
`writer` function must have signature `(io::IO, data)`

# Examples

Sending string content
```julia
text_data = "this is the test"
callback!(app, Output("download", "data"), Input("download-btn", "n_clicks"), prevent_initial_call = true) do n_clicks
    return dcc_send_string(text_data, "text.txt")
end
```

Sending `DataFrame` in `CSV` format
```julia
using DataFrames, CSV
...
df = DataFrame(...)
callback!(app, Output("download", "data"), Input("download-btn", "n_clicks"), prevent_initial_call = true) do n_clicks
    return dcc_send_string(CSV.write, df, "df.csv")
end
```
"""
function dcc_send_string(src::AbstractString, filename; type = nothing)

    return Dict(
        :content => src,
        :filename => filename,
        :type => type,
        :base64 => false
    )
end

function dcc_send_string(writer::Function, data, filename; type = nothing)
    io = IOBuffer()
    writer(io, data)
    return dcc_send_string(String(take!(io)), filename, type = type)
end


text:
# Checks PELT runs

println("Running PELT tests...")

n = 1000;        # Number of samples
λ = 100;         # Frequencey of changes

########################
# Normal mean segments #
########################
μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)
seg_costs = NormalMeanSegment(sample);
PELT(seg_costs, n);
@PELT sample Normal(?,σ)

#######################
# Normal var segments #
#######################
μ, σ = 1.0, Uniform(2.0, 15.0)
sample, cps = @changepoint_sampler n λ Normal(μ, σ)
seg_costs = NormalVarSegment(sample, μ)
PELT(seg_costs, n)
@PELT sample Normal(μ, ?)

############################
# Exponential changepoints #
############################
μ = Uniform(0.0, 10.0)
sample, cps = @changepoint_sampler n λ Exponential(μ)
seg_costs = ExponentialSegment(sample)
PELT(seg_costs, n)
@PELT sample Exponential(?)

########################
# Poisson changepoints #
########################
μ = Uniform(0.0, 10.0)
sample, cps = @changepoint_sampler n λ Poisson(μ)
seg_costs = PoissonSegment(sample)
PELT(seg_costs, n)
@PELT sample Poisson(?)


text:
# Model tests, to ensure that the loaded models are running
# with a decent accuracy.

# Tests files have been taken from ONNX.
# You need ONNX.jl and ProtoBuf.jl installed to run these tests.

using ONNX
using ProtoBuf
using Keras

function read_value(filename)
    return ONNX.get_array(readproto(open(filename), ONNX.Proto.TensorProto()))
end

# Test loop
loss(x, y) = sum((x .- y).^2)
netloss = 0
for x=0:2
    input = read_value("input_$x.pb")
    output = Keras.load("vgg19.json", "vgg19.h5", input)
    op_expected = read_value("output_$x.pb")
    netloss += loss(output, op_expected)
end

print("Average Net loss is ", netloss)


text:
using LinearAlgebra
using Zygote
using ForwardDiff
using ReverseDiff
using StaticArrays
using Base.Iterators

## Finite differencing
central_2nd(φp, φ, φm)     = φm - 2φ + φp
Δ²_central(φp, φ, φm, Δx)  = central_2nd(φp, φ, φm) / Δx^2
Δ_central(φp, φm, Δx)      = (φp - φm) / 2Δx
Δ_forward(φ1, φ2, Δt)      = (φ2 - φ1) / Δt

## Stencil generation
neighbours(CI, deg) = CI .+ CartesianIndices(ntuple(x -> -deg:deg, length(CI)))

stencil_depth(CI, depth) = [ CI + CartesianIndex(ntuple(x -> ifelse(x == i, depth, 0), length(CI))) for i in 1:length(CI) ]

stencils(CI, deg, poke = false) = ifelse(poke, [ stencil_depth(CI, -(deg + 1)); neighbours(CI, deg)[:]; stencil_depth(CI, deg + 1) ], neighbours(CI, deg)[:])

stencil_1(CI) = stencils(CI, 0, true)

## FTCS scheme
ftcs_stencil(φ, index) = @view φ[stencil_1(index)[2:end]]

function substitute(φ, index, val)
    if index ∈ CartesianIndices(φ)
        φ[index]
    else
        val
    end
end

function ftcs_stencil(φ, index, φ_boundaries)
    if any(index.I .<= 1) || any(index.I .>= size(φ)) # There have to be more conditions...
        φ_vec = substitute.(Ref(φ), stencil_1(index), φ_boundaries)
    else
        φ_vec = @view φ[stencil_1(index)[2:end]]
    end
end

## Difference operators specific to FTCS stencil

# First derivatives
Δt(φ, dt) = Δ_forward(φ[4], φ[3], dt)
Δx(φ, dx) = Δ_central(φ[5], φ[1], dx)
Δy(φ, dy) = Δ_central(φ[6], φ[2], dy)

# Second derivatives
Δ²x(φ, dx) = Δ²_central(φ[5], φ[3], φ[1], dx)
Δ²y(φ, dy) = Δ²_central(φ[6], φ[3], φ[2], dy)

# Governing equations in residual form
vorticity_residual(ω, ψ, ν, dt, dx, dy) = Δt(ω, dt) + Δy(ψ, dy) * Δx(ω, dx) - Δx(ψ, dx) * Δy(ω, dy) - ν * (Δ²x(ω, dx) + Δ²y(ω, dy))
streamfunction_residual(ω, ψ, dx, dy)   = Δ²x(ψ, dx) + Δ²y(ψ, dy) + ω[3]

# Compute residual equations
function compute_residuals!(R, p, sizer, ω_bound, ψ_bound, ν, dt, dx, dy)
    ωs = @view p[1:prod(sizer)]
    ψs = @view p[prod(sizer)+1:end]

    ωs = reshape(ωs, sizer)
    ψs = reshape(ψs, sizer)

    # Get indices and bounds
    inds = CartesianIndices(ωs)
    R    = reshape(R, (sizer..., 2))
    # ∂R∂p = reshape(∂R∂p, ((sizer..., sizer...)..., 4))

    # Iterate over grid
    for ind in inds
        # Get neighbours of a cell
        ω = ftcs_stencil(ωs, ind, ω_bound)
        ψ = ftcs_stencil(ψs, ind, ψ_bound)

        # Evaluate residuals
        R[ind,1] = vorticity_residual(ω, ψ, ν, dt, dx, dy) # Vorticity equation
        R[ind,2] = streamfunction_residual(ω, ψ, dx, dy)   # Streamfunction equation
    end

    R[:]
    # , reshape(∂R∂p, (2 .* (prod(size(ωs)), prod((reverse ∘ size)(ωs)))...))
end

newton_update(p, δp, α = 1.0) = p - α * δp

# Pre-allocated version
newton_update!(p, δp, α = 1.0) = @. p = p - α * δp

## Initial setup
dt, dx, dy  = 0.1, 0.2, 0.1
ts, xs, ys  = 0:dt:1, -1:dx:1, -1:dy:1
grid        = product(xs, ys, ts)
α           = 1.0
ν           = 0.05
U           = 2.0
ω_bound     = 1.0
ψ_bound     = 0.1

ωs = ω_bound * ones(size(grid))
ψs = ψ_bound * ones(size(grid))

p = [ ωs[:]; ψs[:] ]

R    = similar(p)
∂R∂p = zeros(2 .* (prod(size(ωs)), prod((reverse ∘ size)(ωs)))...)

num_iters = 3
compute_residuals!(R, p) = compute_residuals!(R, p, size(ωs), ω_bound, ψ_bound, ν, dt, dx, dy)
compute_grad_fwd(R, p) = ForwardDiff.jacobian(compute_residuals!, R, p)

##
# using NLsolve
# prob = nlsolve(compute_residuals!, p)

## NEWTON ITERATIONS
#=================================================================#

function newton_solver!(R, p, num_iters = 3, α = 1.0)
    # Array to store errors
    ε = zeros(num_iters);

    # Newton iteration loop
    for i in 1:num_iters
        # Compute residuals
        R    = compute_residuals!(R, p)

        # Compute Jacobian
        ∂R∂p = compute_grad_fwd(R, p)

        # Compute Newton step
        Δp   = ∂R∂p \ -R

        # Update state with relaxation factor
        p   .= newton_update!(p, Δp, α)

        # Error processing
        ε[i] = maximum(abs.(Δp))
        println("Newton step error L² norm: $(ε[i])")
    end
    R, p, ε
end

##
R, p, ε = newton_solver!(R, p, num_iters, α)

## Plotting
using GLMakie

# function grad_vorticity_residual(ω, ψ, ν, dt, dx, dy)

# end

# grad_streamfunc_residual(ω, ψ, dx, dy) = gradient(streamfunction_residual, ω, ψ, dx, dy)

# function grad_residual_update!(∂R∂p, grad, inds)
#     j, k, l = inds.I
#     a, b, c, d, e, f = size(∂R∂p)
#              ∂R∂p[inds, j  , k  , l  ] = grad[3]
#     if k > 1 ∂R∂p[inds, j  , k-1, l  ] = grad[1] end
#     if l > 1 ∂R∂p[inds, j  , k  , l-1] = grad[2] end
#     if j < a ∂R∂p[inds, j+1, k  , l  ] = grad[4] end
#     if k < b ∂R∂p[inds, j  , k+1, l  ] = grad[5] end
#     if l < c ∂R∂p[inds, j  , k  , l+1] = grad[6] end
# end

# function update_jacobian!(∂R∂p, ind, ω, ψ, ν, dt, dx, dy)
#     # Evaluate Jacobian
#     grad_ω = grad_vorticity_residual(ω, ψ, ν, dt, dx, dy)
#     grad_ψ = grad_streamfunc_residual(ω, ψ, dx, dy)

#     grad_residual_update!(∂R∂p[:,:,:,:,:,:,1], grad_ω[1], ind)
#     grad_residual_update!(∂R∂p[:,:,:,:,:,:,2], grad_ω[2], ind)
#     grad_residual_update!(∂R∂p[:,:,:,:,:,:,3], grad_ψ[1], ind)
#     grad_residual_update!(∂R∂p[:,:,:,:,:,:,4], grad_ψ[2], ind)

#     nothing
# end

max_stars_repo_name: GodotMisogi/AdjointCFD

max_issues_repo_name: GodotMisogi/AdjointCFD

max_forks_repo_name: GodotMisogi/AdjointCFD


text:
@testset "NullArray" begin
  a1 = NullArray{Float64}(1)
  a2 = NullArray{Float64}(1, 2)
  a3 = NullArray{Float64}(1, 2, 3)
  a4 = NullArray{Float64}(1, 2, 3, 4)
  a5 = NullArray{Float64}(1, 2, 3, 4, 5)
  a6 = NullArray{Float64}(1, 2, 3, 4, 5, 6)
  a7 = NullArray{Float64}(1, 2, 3, 4, 5, 6, 7)

  @test size(a1) == (1,)
  @test size(a2) == (1, 2)
  @test size(a3) == (1, 2, 3)
  @test size(a4) == (1, 2, 3, 4)
  @test size(a5) == (1, 2, 3, 4, 5)
  @test size(a6) == (1, 2, 3, 4, 5, 6)
  @test size(a7) == (1, 2, 3, 4, 5, 6, 7)

  arrays = [a1, a2, a3, a4, a5, a6, a7]
  len = 1
  for i=1:length(arrays)
    len *= i
    @test length(arrays[i]) == len
    @test_throws ErrorException arrays[i][1]
    @test_throws ErrorException arrays[i][1] = 1
  end

end

@testset "EmptyArray" begin
  a1 = EmptyArray{Float64}(1)
  a2 = EmptyArray{Float64}(1, 2)
  a3 = EmptyArray{Float64}(1, 2, 3)
  a4 = EmptyArray{Float64}(1, 2, 3, 4)
  a5 = EmptyArray{Float64}(1, 2, 3, 4, 5)
  a6 = EmptyArray{Float64}(1, 2, 3, 4, 5, 6)
  a7 = EmptyArray{Float64}(1, 2, 3, 4, 5, 6, 7)

  @test size(a1) == (1,)
  @test size(a2) == (1, 2)
  @test size(a3) == (1, 2, 3)
  @test size(a4) == (1, 2, 3, 4)
  @test size(a5) == (1, 2, 3, 4, 5)
  @test size(a6) == (1, 2, 3, 4, 5, 6)
  @test size(a7) == (1, 2, 3, 4, 5, 6, 7)

  arrays = [a1, a2, a3, a4, a5, a6, a7]
  len = 1
  v = 0
  for i=1:length(arrays)
    a = arrays[i]
    len *= i
    @test length(a) == len

    # make sure these don't error
    v += a[i]
    a[i] = 1
    a[1] += 1
  end

end


text:
# utext.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files utext.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
"""
module utext end

macro libtext(s)    ; _libicu(s, iculib,     "utext_")    ; end

type UText
    p::Ptr{Void}
    s

    function UText(str::ByteStr)
        err = UErrorCode[0]
        v = Vector{UInt8}(str)
        p = ccall(@libtext(openUTF8), Ptr{Void},
                  (Ptr{Void}, Ptr{UInt8}, Int64, Ptr{UErrorCode}),
                  C_NULL, v, sizeof(v), err)
        @assert SUCCESS(err[1])
        # Retain pointer to data so that it won't be GCed
        self = new(p, v)
        finalizer(self, close)
        self
    end

    function UText(str::UniStr)
        err = UErrorCode[0]
        v = Vector{UInt16}(str)
        p = ccall(@libtext(openUChars), Ptr{Void},
                  (Ptr{Void}, Ptr{UChar}, Int64, Ptr{UErrorCode}),
                  C_NULL, v, length(v)-1, err)
        @assert SUCCESS(err[1])
        # Retain pointer to data so that it won't be GCed
        self = new(p, v)
        finalizer(self, close)
        self
    end
end

close(t::UText) =
    t.p == C_NULL ||
        (ccall(@libtext(close), Void, (Ptr{Void},), t.p) ; t.p = C_NULL ; t.s = Void())


text:
@testset "Trajectories" begin

trajectories = ["Spiral", "Radial", "Cartesian", "EPI", "OneLine", "SpiralVarDens",
                "Cartesian3D", "StackOfStars", "Kooshball"]

for trajName in trajectories
  for T in [Float32, Float64]

    tr = trajectory(T, trajName, 8, 9)
    @test trajName == string(tr)
    #println(kspaceNodes(tr)) # TODO: make real tests
  end
end

end


text:
struct ELM{T <: Number}
    n_features::Int
    n_outputs::Int
    hidden_layer::MutableHiddenLayer{T}

    function ELM{T}(n_features::Int,
                    n_outputs::Int) where {T <: Number}
        hidden_layer = MutableHiddenLayer(T)
        new{T}(n_features, n_outputs, hidden_layer)
    end

    function ELM(n_features::Int,
                 n_outputs::Int)
        ELM{Float64}(n_features, n_outputs)
    end
end

function add_neurons!(elm::T₁,
                      n_neurons::Int,
                      activation_function::T₂) where {T₁ <: ELM,
                                                      T₂ <: ActivationFunction}
    param(::HiddenLayer{T}) where {T} = T

    @assert begin
        isnothing(elm.hidden_layer.HH) &&
        isnothing(elm.hidden_layer.TH)
    end "Call clear! on ELM before adding neurons."

    neurons = Neurons(param(elm.hidden_layer), n_neurons,
                      elm.n_features, activation_function)
    elm.hidden_layer.n_neurons += n_neurons
    push!(elm.hidden_layer.neurons, neurons)
    elm
end

function add_data!(elm::T₁,
                   samples::T₂,
                   targets::T₃,
                   sample_weights::T₄;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractMatrix,
                                                  T₃ <: AbstractMatrix,
                                                  T₄ <: AbstractVector}
    @assert elm.hidden_layer.n_neurons > 0 "Add neurons to ELM before adding data."

    if !is_initialized(elm.hidden_layer)
        initialize!(elm.hidden_layer, elm.n_outputs)
    end

    X = samples
    T = targets
    ψ = sample_weights
    Dₓ, Nₓ = size(X)
    Qₜ, Nₜ = size(T)

    @assert Dₓ == elm.n_features "Input dimensionality mismatch."
    @assert Qₜ == elm.n_outputs "Output dimensionality mismatch."
    @assert Nₓ == Nₜ == length(ψ) "Sample count mismatch."

    N = Nₓ

    for batch in partition_range(1:N, batch_size)
        X₀ = X[:,batch]
        T₀ = T[:,batch]
        ψ₀ = ψ[batch]
        add_batch!(elm, X₀, T₀, ψ₀)
    end

    elm
end

function add_data!(elm::T₁,
                   samples::T₂,
                   targets::T₃;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractMatrix,
                                                  T₃ <: AbstractMatrix}
    N = last(size(samples))
    sample_weights = [1.0 for n in 1:N]
    add_data!(elm, samples, targets, sample_weights, batch_size = batch_size)
end

function add_data!(elm::T₁,
                   samples::T₂,
                   targets::T₃,
                   sample_weights::T₄;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractMatrix,
                                                  T₃ <: AbstractVector,
                                                  T₄ <: AbstractVector}
    add_data!(elm, samples, reshape(targets, 1, :), sample_weights,
              batch_size = batch_size)
end

function add_data!(elm::T₁,
                   samples::T₂,
                   targets::T₃;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractMatrix,
                                                  T₃ <: AbstractVector}
    add_data!(elm, samples, reshape(targets, 1, :), batch_size = batch_size)
end

function add_data!(elm::T₁,
                   sample::T₂,
                   targets::T₃,
                   sample_weights::T₄;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractVector,
                                                  T₃ <: AbstractMatrix,
                                                  T₄ <: AbstractVector}
    add_data!(elm, reshape(sample, 1, :), targets, sample_weights,
              batch_size = batch_size)
end

function add_data!(elm::T₁,
                   sample::T₂,
                   targets::T₃;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractVector,
                                                  T₃ <: AbstractMatrix}
    add_data!(elm, reshape(sample, 1, :), targets, batch_size = batch_size)
end

function add_data!(elm::T₁,
                   sample::T₂,
                   targets::T₃,
                   sample_weights::T₄;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractVector,
                                                  T₃ <: AbstractVector,
                                                  T₄ <: AbstractVector}
    add_data!(elm, reshape(sample, 1, :), reshape(targets, 1, :),
              sample_weights, batch_size = batch_size)
end

function add_data!(elm::T₁,
                   sample::T₂,
                   targets::T₃;
                   batch_size::Int = 1000) where {T₁ <: ELM,
                                                  T₂ <: AbstractVector,
                                                  T₃ <: AbstractVector}
    add_data!(elm, reshape(sample, 1, :). reshape(targets, 1, :),
              batch_size = batch_size)
end

function add_batch!(elm::T₁,
                    samples::T₂,
                    targets::T₃,
                    sample_weights::T₄) where {T₁ <: ELM,
                                               T₂ <: AbstractMatrix,
                                               T₃ <: AbstractMatrix,
                                               T₄ <: AbstractVector}
    X = samples
    T = targets
    Ψ = LinearAlgebra.Diagonal(sample_weights)

    # Get H, the nonlinear representation of X.
    H = project(elm.hidden_layer, X)

    # Adjust the weights of the samples and targets to compensate for any
    # imbalance in the frequency with which categories appear in the dataset.
    H = H * Ψ
    T = T * Ψ

    # Increment `elm.HH` and `elm.TH`, two covariance matrices which preserve
    # the intermediate state of the ELM, before it is solved.
    elm.hidden_layer.HH .+= (H * H')
    elm.hidden_layer.TH .+= (T * H')

    elm
end

function clear!(elm::T) where {T <: ELM}
    elm.hidden_layer.HH = nothing
    elm.hidden_layer.TH = nothing
    elm.hidden_layer.initialized = false
    elm
end

# Find the optimal output weights for the model by multiplying `elm.TH` by the
# pseudoinverse of `elm.HH`, and return the solved model as an `SLFN` object.
function solve(elm::T) where {T <: ELM}
    D = elm.n_features
    Q = elm.n_outputs
    hidden_layer = make_immutable(elm.hidden_layer)
    output_weights = elm.hidden_layer.TH * LinearAlgebra.pinv(elm.hidden_layer.HH)
    SLFN(D, Q, hidden_layer, output_weights)
end


text:
immutable BigFloatA <: AbstractFloat
    prec::Clong
    sign::Cint
    exp::Clong
    signif::Vector{Limb}
end

precision(x::BigFloatA) = x.prec

function new_data(prec)
    k = sizeof(Limb)*8
    p = div(prec + (k-1), k)
    Array(Limb,p)
end

function new_mpfr_data()
    prec = precision(BigFloat)
    s = new_data(prec)
    Mpfr_t(prec,0,0,unsafe_convert(Ptr{Limb}, s)), s
end

BigFloatA(m::Mpfr_t, s::Vector{Limb}) =
    BigFloatA(m.prec, m.sign, m.exp, s)

convert(::Type{Mpfr_t}, x::BigFloatA) =
    Mpfr_t(x.prec, x.sign, x.exp, unsafe_convert(Ptr{Limb},x.signif))

for (fJ, fC) in ((:si,:Clong), (:ui,:Culong), (:d,:Float64))
    @eval begin
        function convert(::Type{BigFloatA}, x::($fC))
            m, s = new_mpfr_data()
            ccall(($(string(:mpfr_set_,fJ)), :libmpfr), Int32, (Ptr{Mpfr_t}, ($fC), Int32), &m, x, ROUNDING_MODE[end])
            return BigFloatA(m, s)
        end
    end
end

# Basic arithmetic without promotion
for (fJ, fC) in ((:+,:add), (:*,:mul))
    @eval begin
        # BigFloatA
        function ($fJ)(x::BigFloatA, y::BigFloatA)
            m, s = new_mpfr_data()
            ccall(($(string(:mpfr_,fC)),:libmpfr), Int32, (Ptr{Mpfr_t}, Ptr{Mpfr_t}, Ptr{Mpfr_t}, Int32), &m, &x, &y, ROUNDING_MODE[end])
            return BigFloatA(m, s)
        end

    end
end

zero(::Type{BigFloatA}) = convert(BigFloatA,0)

function string(x::BigFloatA)
    # In general, the number of decimal places needed to read back the number exactly
    # is, excluding the most significant, ceil(log(10, 2^precision(x)))
    k = ceil(Int32, precision(x) * 0.3010299956639812)
    lng = k + Int32(8) # Add space for the sign, the most significand digit, the dot and the exponent
    buf = Array(UInt8, lng + 1)
    # format strings are guaranteed to contain no NUL, so we don't use Cstring
    lng = ccall((:mpfr_snprintf,:libmpfr), Int32, (Ptr{UInt8}, Culong, Ptr{UInt8}, Ptr{Mpfr_t}), buf, lng + 1, "%.Re", &x)
    if lng < k + 5 # print at least k decimal places
        lng = ccall((:mpfr_sprintf,:libmpfr), Int32, (Ptr{UInt8}, Ptr{UInt8}, Ptr{Mpfr_t}), buf, "%.$(k)Re", &x)
    elseif lng > k + 8
        buf = Array(UInt8, lng + 1)
        lng = ccall((:mpfr_snprintf,:libmpfr), Int32, (Ptr{UInt8}, Culong, Ptr{UInt8}, Ptr{Mpfr_t}), buf, lng + 1, "%.Re", &x)
    end
    return bytestring(pointer(buf), (1 <= x < 10 || -10 < x <= -1 || x == 0) ? lng - 4 : lng)
end


text:
"""
    bin_sparse_matrix!(M::AbstractArray, M_id::AbstractSparseMatrixCSC, max_num_bins::Integer)

Compute a Binning pattern in `M_id` from the sparsity pattern in `M_id`.

`M_id` is `SparseMatrixCSC{Int64, Int64}` with the same shape as `M`,
when handed over it contains only 0 and 1. After usage it contains an Integer
corresponding to a bin or 0. According to the Binning Pattern,
similar valued entries can be considered the same.

`M` is the matrix which entries are used for binning. `max_num_bins` is the maximum
number of bins and must be non-negative, if `max_num_bins=0`, no binning is
performed and only the sparsity pattern in `M_id` is used then.

See also: [`sparsify`](@ref), [`sparsity_pattern`](@ref).

# Examples
```julia-repl
julia> bin_sparse_matrix!([4 1 4.01; 0.1 17.1 17; 0.2 4 29], sparse([1 0 1; 0 1 1; 0 0 1]), 200)
3×3 SparseMatrixCSC{Int64, Int64} with 5 stored entries:
 1  ⋅  1
 ⋅  2  2
 ⋅  ⋅  3
```
"""
function bin_sparse_matrix!(M::AbstractArray{T}, M_id::AbstractSparseMatrixCSC, max_num_bins::Integer) where {T}
    max_num_bins >= 0 || throw(DomainError(max_num_bins, "max_num_bin must be non-negative"))
    size(M, 1) == size(M_id, 1) && size(M, 2) == size(M_id, 2) || throw(ArgumentError("M must be the same size as M_id"))

    m = size(M_id, 1)
    n = size(M_id, 2)
    M_i, M_j = findnz(M_id)
    M_pat_nnz = count(!iszero, M_id)

    separated_at = zero(T)

    if max_num_bins == 0
        M_id[:, :] = sparse(M_i, M_j, 1:M_pat_nnz, m, n)
    else
        nz_M = zeros(T, M_pat_nnz)

        for k in 1:M_pat_nnz
            nz_M[k] = M[M_i[k], M_j[k]]
        end
        min_left, max_left, min_right, max_right = separated_min_max(nz_M, separated_at)
        bin_ids = binned_with_separated_min_max(nz_M, max_num_bins, min_left, max_left, min_right, max_right, separated_at)
        new_ids = bin_mapping(bin_ids)
        M_id[:, :] = sparse(M_i, M_j, new_ids, m, n)
    end
    return M_id
end

function bin_sparse_matrix!(M::AbstractArray{<:Int}, M_id::AbstractSparseMatrixCSC, max_num_bins::Integer)
    return bin_sparse_matrix!(float(M), M_id, max_num_bins)
end

_extrema(v) = isempty(v) ? (Inf, -Inf) : extrema(v)

function separated_min_max(v::AbstractVector{T}, separated_at::T) where {T<:Real}
    min_left, max_left = _extrema(v[v .<= separated_at])
    min_right, max_right = _extrema(v[v .>= separated_at])

    return min_left, max_left, min_right, max_right
end

function binned_with_separated_min_max(v::AbstractVector{T}, max_num_bins::Integer,
    min_left::Real, max_left::Real, min_right::Real, max_right::Real, separated_at::Real) where{T}
    n = length(v)

    if min_left <= max_left
        all(min_left .<= v) || throw(ArgumentError("some element of v is less than min_left"))
        min_left <= separated_at || throw(ArgumentError("incompatible min_left and separated_at"))
    end
    if min_right <= max_right
        all(v .<= max_right) || throw(ArgumentError("some element of v is greater than max_right"))
        separated_at <= max_right || throw(ArgumentError("incompatible separated_at and max_right"))
    end

    if max_left == min_right
        loc_max_num_left_right_bins = max_num_bins - 1
    else
        loc_max_num_left_right_bins = max_num_bins
    end

    if separated_at > min_left
        left_dist = separated_at - min_left
    else
        left_dist = 0
    end

    if max_right > separated_at
        right_dist = max_right - separated_at
    else
        right_dist = 0
    end

    total_dist = left_dist + right_dist

    if total_dist != 0
        max_n_left_bins = floor(Int, (loc_max_num_left_right_bins * left_dist) / total_dist)
        max_n_right_bins = floor(Int, (loc_max_num_left_right_bins * right_dist) / total_dist)
    else
        max_n_left_bins = Inf
        max_n_right_bins = Inf
    end
    @assert max_n_left_bins + max_n_right_bins <= loc_max_num_left_right_bins

    inv_h_l = max_n_left_bins / left_dist
    inv_h_r = max_n_right_bins / right_dist

    fuzz::T = 1E2  # magic constant from matlab

    left_tol = fuzz * eps(T) * left_dist
    right_tol = fuzz * eps(T) * right_dist

    bin_ids = -ones(Int, n)

    middle_id = max_n_left_bins + max_n_right_bins

    for i in 1:n
        if v[i] == separated_at
            b = middle_id
        elseif v[i] < separated_at
            if separated_at - v[i] <= left_tol
                b = middle_id
            else
                b = floor(Int, (v[i] - min_left) * inv_h_l)
            end
        else  # if (separated_at < v[i])
            if v[i] - separated_at <= right_tol
                b = middle_id
            else
                b = floor(Int, (max_right - v[i]) * inv_h_r) + max_n_left_bins
            end
        end
        bin_ids[i] = b
    end
    return bin_ids
end

function bin_mapping(bin_ids::AbstractVector{<:Integer})

    min_id, max_id = _extrema(bin_ids)
    impossible_id = typemin(Int)

    impossible_id != min_id || throw(DomainError("input array must not contain the minimum int"))

    num_bins = max_id - min_id + 1
    work_array = impossible_id * ones(Int, num_bins)
    n_unique_vals = 0

    n = length(bin_ids)

    for i in 1:n
        id = work_array[bin_ids[i] - min_id + 1]

        if id == impossible_id
            n_unique_vals = n_unique_vals + 1
            work_array[bin_ids[i] - min_id + 1] = n_unique_vals
        end
    end
    new_bin_ids = copy(bin_ids)
    for i in 1:n
        new_bin_ids[i] = work_array[bin_ids[i] - min_id + 1]
    end
    return new_bin_ids
end


text:
type U <: Component
	x
	units
end
Physical=U
*(p1::Physical,p2::Physical)=simplify(Physical(p1.x*p2.x,p1.units*p2.units))
+(p1::Physical,p2::Physical)=simplify(Physical(p1.x+p2.x,p1.units+p2.units))
-(p1::Physical,p2::Physical)=simplify(Physical(p1.x-p2.x,p1.units-p2.units))
/(p1::Physical,p2::Physical)=simplify(Physical(p1.x/p2.x,p1.units/p2.units))
*(p::Physical,n::Number)=simplify(Physical(p.x*n,p.units))
function simplify(p::Physical)
	u=simplify(p.units)
	if isa(u,Array) && length(u)==1 && length(u[1])>1 && isa(u[1][1],Number)
		shift!(u[1])
	end
	return U(simplify(p.x),u)
end


text:
for op in (:exp, :expm1, :log, :log1p,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth)
    @eval begin

        function ($op)(a::Float65)
            state = getstate(a)
            value = Float65( ($op)(reflect(a.fp)) )
            if state
               value = setstate(value)
            end
            value
         end

    end
end


text:
using miniWeather
using Test

@testset "test mass and energy conservation" begin
    nx = 100
    nz = 50
    sim_time = 400
    out_freq = 400
    weather_type = "thermal"

    model, grid = init(nx, nz, sim_time, weather_type);
    mass0, te0 = reductions(model, grid)
    run!(model, grid, out_freq, "output.nc")
    mass, te = reductions(model, grid)

    @test (mass - mass0) / mass0 < 1.0e-9
    @test (te - te0) / te0 < 4.5e-5

end


text:
"""
```
struct VarName{sym}
    indexing  ::    String
end
```

A variable identifier. Every variable has a symbol `sym` and `indices `indexing`.
The Julia variable in the model corresponding to `sym` can refer to a single value or
to a hierarchical array structure of univariate, multivariate or matrix variables. `indexing` stores the indices that can access the random variable from the Julia
variable.

Examples:

- `x[1] ~ Normal()` will generate a `VarName` with `sym == :x` and `indexing == "[1]"`.
- `x[:,1] ~ MvNormal(zeros(2))` will generate a `VarName` with `sym == :x` and
 `indexing == "[Colon(),1]"`.
- `x[:,1][2] ~ Normal()` will generate a `VarName` with `sym == :x` and
 `indexing == "[Colon(),1][2]"`.
"""
struct VarName{sym}
    indexing::String
end

"""
    @varname(var)

A macro that returns an instance of `VarName` given the symbol or expression of a Julia variable, e.g. `@varname x[1,2][1+5][45][3]` returns `VarName{:x}("[1,2][6][45][3]")`.
"""
macro varname(expr::Union{Expr, Symbol})
    expr |> varname |> esc
end
function varname(expr)
    ex = deepcopy(expr)
    (ex isa Symbol) && return quote
        DynamicPPL.VarName{$(QuoteNode(ex))}("")
    end
    (ex.head == :ref) || throw("VarName: Mis-formed variable name $(expr)!")
    inds = :(())
    while ex.head == :ref
        if length(ex.args) >= 2
            strs = map(x -> :($x === (:) ? "Colon()" : string($x)), ex.args[2:end])
            pushfirst!(inds.args, :("[" * join($(Expr(:vect, strs...)), ",") * "]"))
        end
        ex = ex.args[1]
        isa(ex, Symbol) && return quote
            DynamicPPL.VarName{$(QuoteNode(ex))}(foldl(*, $inds, init = ""))
        end
    end
    throw("VarName: Mis-formed variable name $(expr)!")
end

macro vsym(expr::Union{Expr, Symbol})
    expr |> vsym
end

"""
    vsym(expr::Union{Expr, Symbol})

Returns the variable symbol given the input variable expression `expr`. For example, if the input `expr = :(x[1])`, the output is `:x`.
"""
function vsym(expr::Union{Expr, Symbol})
    ex = deepcopy(expr)
    (ex isa Symbol) && return QuoteNode(ex)
    (ex.head == :ref) || throw("VarName: Mis-formed variable name $(expr)!")
    while ex.head == :ref
        ex = ex.args[1]
        isa(ex, Symbol) && return QuoteNode(ex)
    end
    throw("VarName: Mis-formed variable name $(expr)!")
end

"""
    @vinds(expr)

Returns a tuple of tuples of the indices in `expr`. For example, `@vinds x[1,:][2]` returns
`((1, Colon()), (2,))`.
"""
macro vinds(expr::Union{Expr, Symbol})
    expr |> vinds |> esc
end
function vinds(expr::Union{Expr, Symbol})
    ex = deepcopy(expr)
    inds = Expr(:tuple)
    (ex isa Symbol) && return inds
    (ex.head == :ref) || throw("VarName: Mis-formed variable name $(expr)!")
    while ex.head == :ref
        pushfirst!(inds.args, Expr(:tuple, ex.args[2:end]...))
        ex = ex.args[1]
        isa(ex, Symbol) && return inds
    end
    throw("VarName: Mis-formed variable name $(expr)!")
end

"""
    split_var_str(var_str, inds_as = Vector)

This function splits a variable string, e.g. `"x[1:3,1:2][3,2]"` to the variable's symbol `"x"` and the indexing `"[1:3,1:2][3,2]"`. If `inds_as = String`, the indices are returned as a string, e.g. `"[1:3,1:2][3,2]"`. If `inds_as = Vector`, the indices are returned as a vector of vectors of strings, e.g. `[["1:3", "1:2"], ["3", "2"]]`.
"""
function split_var_str(var_str, inds_as = Vector)
    ind = findfirst(c -> c == '[', var_str)
    if inds_as === String
        if ind === nothing
            return var_str, ""
        else
            return var_str[1:ind-1], var_str[ind:end]
        end
    end
    @assert inds_as === Vector
    inds = Vector{String}[]
    if ind === nothing
        return var_str, inds
    end
    sym = var_str[1:ind-1]
    ind = length(sym)
    while ind < length(var_str)
        ind += 1
        @assert var_str[ind] == '['
        push!(inds, String[])
        while var_str[ind] != ']'
            ind += 1
            if var_str[ind] == '['
                ind2 = findnext(c -> c == ']', var_str, ind)
                push!(inds[end], strip(var_str[ind:ind2]))
                ind = ind2+1
            else
                ind2 = findnext(c -> c == ',' || c == ']', var_str, ind)
                push!(inds[end], strip(var_str[ind:ind2-1]))
                ind = ind2
            end
        end
    end
    return sym, inds
end

@generated function inargnames(::VarName{s}, ::Model{_F, argnames}) where {s, argnames, _F}
    return s in argnames
end

@generated function inmissings(::VarName{s}, ::Model{_F, _a, _T, missings}) where {s, missings, _F, _a, _T}
    return s in missings
end