text:
# ## Stock market data
#
# Let's load the usual packages and the data

using MLJ, RDatasets, ScientificTypes,
      DataFrames, Statistics, StatsBase

smarket = dataset("ISLR", "Smarket")
@show size(smarket)
@show names(smarket)

# Let's get a description too

describe(smarket, :mean, :std, :eltype)

# The target variable is `:Direction`:

y = smarket.Direction
X = select(smarket, Not(:Direction));

# We can compute all the pairwise correlations; we use `Matrix` so that the dataframe entries are considered as one matrix of numbers (otherwise `cor` won't work):

cm = X |> Matrix |> cor
round.(cm, sigdigits=1)

# Let's see what the `:Volume` feature looks like:

using PyPlot
figure(figsize=(8,6))
plot(X.Volume)
xlabel("Tick number", fontsize=14)
ylabel("Volume", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

savefig("assets/ISL-volume.svg") # hide

# ![volume](/assets/ISL-volume.svg)

# ### Logistic Regression
#
# We will now try to train models; the target `:Direction` has two classes: `Up` and `Down`; it needs to be interpreted as a Multiclass object first:

yc = coerce(y, Multiclass)
unique(yc)

# Let's now try fitting a simple logistic classifier (aka logistic regression) not using `:Year` and `:Today`:

@load LogisticClassifier pkg=MLJLinearModels
X2 = select(X, Not([:Year, :Today]))
clf = machine(LogisticClassifier(), X2, y)

# Let's fit it to the data and try to reproduce the output:

fit!(clf)
ŷ = predict(clf, X2)
cross_entropy(ŷ, y) |> mean

# Note that here the `ŷ` are _scores_; in order to recover the class, we could use the mode and compare the misclassification rate:

ŷ = predict_mode(clf, X2)
misclassification_rate(ŷ, y)

# Well that's not fantastic...
#
# Let's visualise how we're doing building a confusion matrix manually,
# first is predicted, second is truth:

TN = down_down = sum(ŷ .== y .== "Down")
FN = down_up = sum(ŷ .!= y .== "Up")
FP = up_down = sum(ŷ .!= y .== "Down")
TP = up_up = sum(ŷ .== y .== "Up")

conf_mat = [down_down down_up; up_down up_up]

# We can then compute the accuracy or precision easily for instance:

acc = (TN + TP) / length(y)
prec = TP /  (TP + FP)
rec  = TP / (TP + FN)
@show round(acc, sigdigits=3)
@show round(prec, sigdigits=3)
@show round(rec, sigdigits=3)

# Let's now train on the data before 2005 and use it to predict on the rest.
# Let's find the row indices for which the condition holds

train = 1:findlast(X.Year .< 2005);
test = last(train)+1:length(y)

# We can now just re-fit the machine that we've already defined just on those rows and predict on the test:

fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
mcr = misclassification_rate(ŷ, y[test])
accuracy = 1 - mcr

# Let's retrain a machine using only `:Lag1` and `:Lag2`:

X3 = select(X2, [:Lag1, :Lag2])
clf = machine(LogisticClassifier(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
mean(ŷ .== y[test])

# Interesting... it has higher accuracy than the model with more features! This could be investigated further by increasing the regularisation parameter but we'll leave it aside for  now.
#
# We can use a trained machine to predict on new data:

Xnew = (Lag1 = [1.2, 1.5], Lag2 = [1.1, -0.8])
@show ŷ = predict(clf, Xnew)

# Note: when specifying data, we used a simple `NamedTuple`; we could also have defined a dataframe or any other compatible tabular container.
# Note also that we retrieved the raw predictions here i.e.: a score for each class; we could have used `predict_mode` or indeed

mode.(ŷ)

# ### LDA
#
# Let's do a similar thing but with a LDA model this time:

@load BayesianLDA pkg=MultivariateStats

clf = machine(BayesianLDA(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

# Note: `BayesianLDA` is LDA using a multivariate normal model for each class with a default prior inferred from the proportions for each class in the training data.
# You can also use the bare `LDA` model which does not make these assumptions and allows using a different metric in the transformed space, see the docs for details.

@load LDA pkg=MultivariateStats
using Distances

clf = machine(LDA(dist=CosineDist()), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

# ### QDA
#
# Bayesian QDA is available via ScikitLearn:

@load BayesianQDA pkg=ScikitLearn

# Using it is done in much the same way as before:

clf = machine(BayesianQDA(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

# ### KNN
#
# Multiple packages offer KNN, we go via NearestNeighbors:

@load KNNClassifier pkg=NearestNeighbors

knnc = KNNClassifier(K=1)
clf = machine(knnc, X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
@show mean(ŷ .== y[test])

# Let's try with three neighbors

knnc.K = 3
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
@show mean(ŷ .== y[test])

# ## Caravan insurance data
#
# The caravan dataset is part of ISLR as well:

caravan  = dataset("ISLR", "Caravan")
@show size(caravan)

# The target variable is `Purchase`, effectively  a categorical

purchase = caravan.Purchase
vals     = unique(purchase)

# Let's see how many of each we have

nl1 = sum(purchase .== vals[1])
nl2 = sum(purchase .== vals[2])
println("#$(vals[1]) ", nl1)
println("#$(vals[2]) ", nl2)

# that's quite unbalanced.
#
# Apart from the target, all other variables are numbers; we can standardize the data:

y, X = unpack(caravan, ==(:Purchase), col->true)

std = machine(Standardizer(), X)
fit!(std)
Xs = transform(std, X)

var(Xs[:,1])

# **Note**: in MLJ, it is recommended to work with pipelines / networks when possible and not do "step-by-step" transformation and fitting of the data as this is more error prone. We do it here to stick to the ISL tutorial.
#
# We split the data in the first 1000 rows for testing and the rest for training:

test = 1:1000
train = last(test)+1:nrows(Xs);

# Let's now fit a KNN model and check the misclassification rate

clf = machine(KNNClassifier(K=3), Xs, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

misclassification_rate(ŷ, y[test])

# that looks good but recall the problem is very unbalanced

mean(y[test] .!= "No")

# Let's fit a logistic classifier to this problem

clf = machine(LogisticClassifier(), Xs, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

misclassification_rate(ŷ, y[test])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
import Base: size, eachindex, getindex, view, ==, eltype, convert, zero, getproperty
import Distances: pairwise
import Zygote: literal_getproperty, accum
export ColVecs

"""
    ColVecs{T, TX<:AbstractMatrix}

A lightweight box for an `AbstractMatrix` to make it behave like a vector of vectors.
"""
struct ColVecs{T, TX<:AbstractMatrix{T}} <: AbstractVector{Vector{T}}
    X::TX
    ColVecs(X::TX) where {T, TX<:AbstractMatrix{T}} = new{T, TX}(X)
end

@adjoint function ColVecs(X::AbstractMatrix)
    back(Δ::NamedTuple) = (Δ.X,)
    back(Δ::AbstractMatrix) = (Δ,)
    function back(Δ::AbstractVector{<:AbstractVector{<:Real}})
        throw(error("In slow method"))
    end
    return ColVecs(X), back
end

==(D1::ColVecs, D2::ColVecs) = D1.X == D2.X
size(D::ColVecs) = (size(D.X, 2),)
getindex(D::ColVecs, n::Int) = D.X[:, n]
getindex(D::ColVecs, n::CartesianIndex{1}) = getindex(D, n[1])
getindex(D::ColVecs, n) = ColVecs(D.X[:, n])
view(D::ColVecs, n::Int) = view(D.X, :, n)
view(D::ColVecs, n) = ColVecs(view(D.X, :, n))
eltype(D::ColVecs{T}) where T = Vector{T}
zero(D::ColVecs) = ColVecs(zero(D.X))

################################ Fancy block data set type #################################

"""
    BlockData{T, TV<:AbstractVector{T}, TX<:AbstractVector{TV}} <: AbstractVector{T}

A strictly ordered collection of `AbstractVector`s, representing a ragged array of data.
"""
struct BlockData{T, V<:AbstractVector{<:T}} <: AbstractVector{T}
    X::Vector{V}
end

BlockData(X::Vector{AbstractVector}) = BlockData{Any, AbstractVector}(X)
==(D1::BlockData, D2::BlockData) = D1.X == D2.X
size(D::BlockData) = (sum(length, D.X),)

blocks(X::BlockData) = X.X
function getindex(D::BlockData, n::Int)
    b = 1
    while n > length(D.X[b])
        n -= length(D.X[b])
        b += 1
    end
    return D.X[b][n]
end
function getindex(D::BlockData, n::BlockVector{<:Integer})
    @assert eachindex(D) == n
    return D
end
view(D::BlockData, b::Int, n) = view(D.X[b], n)
eltype(D::BlockData{T}) where T = T
function eachindex(D::BlockData)
    lengths = map(length, blocks(D))
    return BlockArray(1:sum(lengths), lengths)
end

ext: jl


text:
const _maxConstName = 33
const _maxConstValue = 1024

"""
   constants(eph)

   Retrieve the constants stored in the ephemeris associated to handle eph as a dictionary

"""
function constants(eph::Ephem)
    res = Dict{Symbol,Any}()
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    NC::Int = ccall((:calceph_getconstantcount , libcalceph), Cint,
    (Ptr{Cvoid},),eph.data)
    value = Ref{Cdouble}(0.0)
    name = Vector{UInt8}(undef,_maxConstName)
    for i=1:NC
       numberOfValues = ccall((:calceph_getconstantindex , libcalceph), Cint,
                  (Ptr{Cvoid},Cint,Ptr{UInt8},Ref{Cdouble}),
                  eph.data, i ,name ,value)
       if (numberOfValues==1)
          res[Symbol(strip(unsafe_string(pointer(name))))] = value[]
       elseif (numberOfValues>1)
          values = Array{Float64,1}(undef,numberOfValues)
          stat = ccall((:calceph_getconstantvd , libcalceph), Cint,
                  (Ptr{Cvoid},Ptr{UInt8},Ptr{Cdouble}, Cint),
                  eph.data, name ,values, numberOfValues)
          if (stat>0)
             res[Symbol(strip(unsafe_string(pointer(name))))] = values
          end
       else
          numberOfValues = ccall((:calceph_getconstantvs , libcalceph), Cint,
                  (Ptr{Cvoid},Ptr{UInt8},Ptr{Ptr{Char}}, Cint),
                  eph.data, name ,C_NULL, 0)
          if (numberOfValues>0)
             storage = Array{UInt8}(undef,_maxConstValue,numberOfValues)
             stat = ccall((:calceph_getconstantvs , libcalceph), Cint,
             (Ptr{Cvoid},Ptr{UInt8},Ptr{UInt8}, Cint),
             eph.data, name ,storage, numberOfValues)
             if (stat>0)
                 values = [ strip(unsafe_string(pointer(storage,i))) for i in 1:_maxConstValue:length(storage) ]
                 if (numberOfValues==1)
                    res[Symbol(strip(unsafe_string(pointer(name))))] = values[1]
                 else
                    res[Symbol(strip(unsafe_string(pointer(name))))] = values
                 end
             end
          end
       end
    end

    return res
end


text:
"""
    getlocalizations()

Extract `Localization`s from the data file for the given channel.
"""
function getlocalizations(alllocalizations::Vector{Localization}, channelname, startframe, nframes,
    starttrimframes, endtrimframes, radius, t_off, mergeradius)

    lowerlimit = startframe + starttrimframes
    upperlimit = startframe + nframes - endtrimframes - 1

    localizations = filter(l -> l.channel == channelname, alllocalizations)
end

"""
    getmolecules()

Group `Localization`s into `Molecule`s, using temporal and spatial proximity.
"""
function getmolecules(alllocalizations::Vector{Localization}, channelname, startframe, nframes,
    starttrimframes, endtrimframes, radius, t_off, mergeradius)

    lowerlimit = startframe + starttrimframes
    upperlimit = startframe + nframes - endtrimframes - 1

    localizations = filter(l -> l.channel == channelname, alllocalizations)
    localizationstrimmed = filter(l -> within(l.frame, lowerlimit, upperlimit), localizations)
    molecules = groupby_localmax_temporallimit(localizationstrimmed, radius, t_off)
    moleculesmerged = merge_close_molecules(molecules, mergeradius)

    moleculesmerged, localizationstrimmed
end


text:
module Evolutionary
    using Random, LinearAlgebra, Statistics
    import NLSolversBase: AbstractObjective, NonDifferentiable, ConstraintBounds,
           value, value!, nconstraints_x
    import Base: length, push!, show, @kwdef, minimum, copy, identity
    import UnPack.@unpack
    import NLSolversBase: NonDifferentiable, f_calls
    export AbstractStrategy, strategy, mutationwrapper,
           IsotropicStrategy, AnisotropicStrategy, NoStrategy,
           # ES mutations
           gaussian, cauchy,
           # GA mutations
           flip, bitinversion, domainrange, inversion, insertion, swap2, scramble, shifting,
           # ES recombinations
           average, marriage,
           # GA recombinations
           singlepoint, twopoint, uniform,
           discrete, waverage, intermediate, line,
           pmx, ox1, cx, ox2, pos,
           # GA selections
           ranklinear, uniformranking, roulette, rouletteinv, sus, susinv, tournament, truncation,
           # DE selections
           random, permutation, randomoffset, best,
           # DE recombinations
           uniformbin, exponential,
           # Optimization methods
           ES, CMAES, GA, DE,
           es, cmaes, ga,
           NonDifferentiable

    # optimize API
    include("api/types.jl")
    include("api/results.jl")
    include("api/utilities.jl")
    include("api/optimize.jl")

    # ES & GA recombination functions
    include("recombinations.jl")

    # ES & GA mutation functions
    include("mutations.jl")

    # GA selection functions
    include("selections.jl")

    # Evolution Strategy
    include("es.jl")
    include("cmaes.jl")

    # Genetic Algorithms
    include("ga.jl")

    # Differential Evolution
    include("de.jl")

    # deprecations
    @deprecate isotropic(recombinant, strategy) gaussian(recombinant, strategy)
    @deprecate anisotropic(recombinant, strategy) gaussian(recombinant, strategy)
    @deprecate isotropicSigma(strategy) gaussian(strategy)
    @deprecate anisotropicSigma(strategy) gaussian(strategy)
    @deprecate averageSigma(strategy) average(strategy)
    @deprecate strategy(kwargs...) IsotropicStrategy(N) false
end


text:
push!(Base.LOAD_PATH, joinpath(dirname(@__FILE__), "..", ".."))

using Test
import ProxSDP
import MathOptInterface
import LinearAlgebra
import Random
import SparseArrays
import DelimitedFiles

const MOI = MathOptInterface
const MOIB = MOI.Bridges
const MOIU = MOI.Utilities

const optimizer_bridged = MOI.instantiate(
    ()->ProxSDP.Optimizer(
        tol_gap = 1e-6, tol_feasibility= 1e-6,
        # max_iter = 100_000,
        time_limit = 3., #seconds FAST
        warn_on_limit = true,
        # log_verbose = true, log_freq = 100000
        ),
    with_bridge_type = Float64)

@testset "SolverName" begin
    @test MOI.get(optimizer_bridged, MOI.SolverName()) == "ProxSDP"
end

@testset "SolverVersion" begin
    ver = readlines(joinpath(@__DIR__, "..", "Project.toml"))[4][12:16]
    @test MOI.get(optimizer_bridged, MOI.SolverVersion()) == ver
end

function test_runtests()

    config = MOI.Test.Config(
        atol = 1e-4,
        rtol = 1e-3,
        exclude = Any[
            MOI.ConstraintBasisStatus,
            MOI.VariableBasisStatus,
            MOI.ConstraintName,
            MOI.VariableName,
            MOI.ObjectiveBound,
            MOI.ScalarFunctionConstantNotZero, # ignored by UniversalFallback
        ],
    )

    opt = ProxSDP.Optimizer(
        tol_gap = 1e-6, tol_feasibility= 1e-6,
        # max_iter = 100_000,
        time_limit = 1.0, #seconds FAST
        warn_on_limit = true,
        # log_verbose = true, log_freq = 100000
        )
    model = MOI.instantiate(()->opt, with_bridge_type = Float64)
    MOI.set(model, MOI.Silent(), true)

    MOI.Test.runtests(
        model,
        config,
        exclude = String[
            # Not ProxSDP fail.
            # see: https://github.com/jump-dev/MathOptInterface.jl/issues/1665
            "test_model_UpperBoundAlreadySet",
            "test_model_LowerBoundAlreadySet",
            # poorly scaled problem (solved bellow with higher accuracy)
            "test_linear_add_constraints",
            # time limit hit
            "test_linear_INFEASIBLE",
            "test_solve_DualStatus_INFEASIBILITY_CERTIFICATE_VariableIndex_LessThan_max",
        ],
    )

    MOI.set(model, MOI.RawOptimizerAttribute("time_limit"), 5.0)
    MOI.empty!(model)
    MOI.Test.test_solve_DualStatus_INFEASIBILITY_CERTIFICATE_VariableIndex_LessThan_max(
        model,
        config,
    )
    MOI.empty!(model)
    MOI.Test.test_linear_INFEASIBLE(
        model,
        config,
    )

    MOI.set(model, MOI.RawOptimizerAttribute("tol_primal"), 1e-7)
    MOI.set(model, MOI.RawOptimizerAttribute("tol_dual"), 1e-7)
    MOI.set(model, MOI.RawOptimizerAttribute("tol_gap"), 1e-7)
    MOI.set(model, MOI.RawOptimizerAttribute("tol_feasibility"), 1e-7)

    MOI.set(model, MOI.Silent(), true)
    MOI.empty!(model)
    MOI.Test.test_linear_add_constraints(
        model,
        config,
    )

    @test MOI.get(model, ProxSDP.PDHGIterations()) >= 0

    return
end

@testset "MOI Unit" begin
    test_runtests()
end

@testset "ProxSDP MOI Units tests" begin
    include("moi_proxsdp_unit.jl")
end

@testset "MIMO Sizes" begin
    include("base_mimo.jl")
    include("moi_mimo.jl")
    for i in 2:5
        @testset "MIMO n = $(i)" begin
            moi_mimo(optimizer_bridged, 123, i, test = true)
        end
    end
end

# hitting time limit
# probably infeasible/unbounded
# @testset "RANDSDP Sizes" begin
#     include("base_randsdp.jl")
#     include("moi_randsdp.jl")
#     for n in 10:11, m in 10:11
#         @testset "RANDSDP n=$n, m=$m" begin
#             moi_randsdp(optimizer, 123, n, m, test = true, atol = 1e-1)
#         end
#     end
# end

# This problems are too large for Travis
opt_low_acc = ProxSDP.Optimizer(
    tol_gap = 1e-3, tol_feasibility = 1e-3,
    # log_verbose = true, log_freq = 1000
    )
cache = MOI.default_cache(opt_low_acc, Float64)
optimizer_low_acc = MOI.Utilities.CachingOptimizer(cache, opt_low_acc)
@testset "SDPLIB Sizes" begin
    datapath = joinpath(dirname(@__FILE__), "data")
    include("base_sdplib.jl")
    include("moi_sdplib.jl")
    @testset "EQPART" begin
        # badly conditioned
        # moi_sdplib(optimizer_low_acc, joinpath(datapath, "gpp124-1.dat-s"), test = true)
        moi_sdplib(optimizer_low_acc, joinpath(datapath, "gpp124-2.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "gpp124-3.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "gpp124-4.dat-s"), test = true)
    end
    @testset "MAX CUT" begin
        moi_sdplib(optimizer_low_acc, joinpath(datapath, "mcp124-1.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "mcp124-2.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "mcp124-3.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "mcp124-4.dat-s"), test = true)
    end
end

@testset "Sensor Localization" begin
    include("base_sensorloc.jl")
    include("moi_sensorloc.jl")
    for n in 5:5:10
        moi_sensorloc(optimizer_bridged, 0, n, test = true)
    end
end

@testset "Unsupported argument" begin
    @test_throws ErrorException optimizer_unsupportedarg = ProxSDP.Optimizer(unsupportedarg = 10)
    # @test_throws ErrorException MOI.optimize!(optimizer_unsupportedarg)
end

include("test_terminationstatus.jl")


text:
"""
    evalbezier(control, t)

Evaluate a Bezier curve using Horner's method.
"""
function evalbezier(control::AbstractArray{T,2}, t::T) where {T<:AbstractFloat}
    # Horner's method
    tᶜ = 1.0 - t
    tᵢ = 1.0
    cᵢ = 1.0
    bᵢ = control[:, 1]*tᶜ
    n = size(control, 2) - 1
    for i in 1:n-1
        tᵢ = tᵢ*t
        cᵢ = (n-i+1)*cᵢ/i
        bᵢ = (bᵢ+ cᵢ*control[:, i+1]*tᵢ)*tᶜ
    end
    bᵢ + control[:, n+1]*tᵢ*t
end

"""
    evalbezier(control, tvec)

Evaluate along at array of time instants.
"""
function evalbezier(control::AbstractArray{T,2}, tvec::AbstractArray{T, 1}) where {T<:AbstractFloat}
    # TODO Just a hack (only really meant for viz)
    p = zeros(size(control, 1), size(tvec, 1))
    for (i,t) in enumerate(tvec)
        p[:,i] = evalbezier(control, t)
    end
    p
end


text:
module TestUtilities

using MLJIteration
using MLJBase
using ..DummyModel
using Test

@testset "get and set iteration paramter" begin
    model = DummyIterativeModel()
    X, y = make_dummy(N=2)
    mach = machine(model, X, y)
    iter = MLJBase.iteration_parameter(model)
    @test MLJIteration.rget(model, iter) == 10
    @test MLJIteration.rget(mach, iter) == 10
    MLJIteration.rset!(model, iter, 42)
    @test MLJIteration.rget(model, iter) == 42
    @test MLJIteration.rget(mach, iter) == 42
    MLJIteration.rset!(mach, iter, 17)
    @test MLJIteration.rget(model, iter) == 17
    @test MLJIteration.rget(mach, iter) == 17
end

end

true


text:
using ArgParse

include("../src/InsarLOS.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--outfile", "-o"
            default = "los_map.h5"
            help = "Name output .h5 file to save solution stack"
        "--full-elevation-dem"
            help = "File with full elevation.dem used to make .geo SLCs"
        "--dbfile"
            help = "Name output .db file to use for satllite params"
        "--dbpath"
            help = "Path to directory contining some relevant .db files"
        "--demrsc"
            default = "dem.rsc"
            help = "Name of .rsc file for igrams- the output LOS file will match this DEM area"
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        @show arg, val
    end

    # Assign the appropriate dataset name based on the order
    InsarLOS.create_los_map(outfile=parsed_args["outfile"],
                            demrsc=parsed_args["demrsc"],
                            dbpath=parsed_args["dbpath"],
                            dbfile=parsed_args["dbfile"],
                            full_elevation_dem=parsed_args["full-elevation-dem"])
end


text:
function obs_env(M,Au,Ad)
    _, FL = leftenv(Au, Ad, M)
    _, FR = rightenv(Au, Ad, M)
    _, FL_n = norm_FL(Au, Ad)
    _, FR_n = norm_FR(Au, Ad)
    M, Au, Ad, FL, FR, FL_n, FR_n
end

"""
    Z(env)

return the up and down partition function of the `env`.
"""
function Z(env)
    M, Au, Ad, FL, FR, FL_n, FR_n = env
    z = ein"(((adf,abc),dgeb),ceh),fgh -> "(FL,Au,M,FR,Ad)[]
    λ = ein"abc,abc -> "(FL,FR)[]
    overlap = ein"(ad,acb),(dce,be) ->"(FL_n,Au,Ad,FR_n)[]/ein"ab,ab ->"(FL_n,FR_n)[]
    return z/λ/overlap
end

"""
    magnetisation(env, model::MT, β)

return the up and down magnetisation of the `model`. Requires that `mag_tensor` are defined for `model`.
"""
function magnetisation(env, model::MT) where {MT <: HamiltonianModel}
    M, Au, Ad, FL, FR, _, _ = env
    Mag = _arraytype(M)(mag_tensor(model))
    mag = ein"(((adf,abc),dgeb),ceh),fgh -> "(FL,Au,Mag,FR,Ad)[]
    λ = ein"(((adf,abc),dgeb),ceh),fgh -> "(FL,Au,M,FR,Ad)[]
    magnet = mag/λ
    @assert imag(magnet) < 1e-6
    return real(magnet)
end

"""
    energy(env, model::MT, β)

return the up and down energy of the `model` as a function of the inverse
temperature `β` and the environment bonddimension `D` as calculated with
vumps. Requires that `model_tensor` are defined for `model`.
"""
function energy(env,model::MT) where {MT <: HamiltonianModel}
    M, Au, Ad, FL, FR, _, _ = env
    Ene = _arraytype(M)(energy_tensor(model))
    energy = ein"(((adf,abc),dgeb),ceh),fgh -> "(FL,Au,Ene,FR,Ad)[]
    λ = ein"(((adf,abc),dgeb),ceh),fgh -> "(FL,Au,M,FR,Ad)[]
    E = energy/λ*2 # factor 2 for counting horizontal and vertical links
    @assert imag(E) < 1e-6
    return real(E)
end

function overlap(env)
    _, Au, Ad, _, _, FL_n, FR_n = env
    _, FLu_n = norm_FL(Au, Au)
    _, FRu_n = norm_FR(Au, Au)
    _, FLd_n = norm_FL(Ad, Ad)
    _, FRd_n = norm_FR(Ad, Ad)
    Au /= sqrt(ein"(ad,acb),(dce,be) ->"(FLu_n,Au,Au,FRu_n)[]/ein"ab,ab ->"(FLu_n,FRu_n)[])
    Ad /= sqrt(ein"(ad,acb),(dce,be) ->"(FLd_n,Ad,Ad,FRd_n)[]/ein"ab,ab ->"(FLd_n,FRd_n)[])
    ein"(ad,acb),(dce,be) ->"(FL_n,Au,Ad,FR_n)[]/ein"ab,ab ->"(FL_n,FR_n)[]
end

"""
    magofβ(::Ising,β)

return the analytical result for the magnetisation at inverse temperature
`β` for the 2d classical ising model.
"""
magofβ(model::Ising) = model.β > isingβc ? (1-sinh(2*model.β)^-4)^(1/8) : 0.

"""
    magofdβ(::Ising,β)

return the analytical result for the derivative of magnetisation at inverse temperature
`β` for the 2d classical ising model.
"""
magofdβ(model::Ising) = model.β > isingβc ? (coth(2*model.β)*csch(2*model.β)^4)/(1-csch(2*model.β)^4)^(7/8) : 0.

"""
    eneofβ(::Ising,β)

return some the numerical integrations of analytical result for the energy at inverse temperature
`β` for the 2d classical ising model.
"""
function eneofβ(model::Ising)
    β = model.β
    if β == 0.0
        return 0
    elseif β == 0.2
        return -0.42822885693016843
    elseif β == 0.4
        return -1.1060792706185651
    elseif β == 0.6
        return -1.909085845408498
    elseif β == 0.8
        return -1.9848514445364174
    elseif β == 1.0
        return -1.997159425837225
    end
end

"""
    Zofβ(::Ising,β)

return some the numerical integrations of analytical result for the partition function at inverse temperature
`β` for the 2d classical ising model.
"""
function Zofβ(model::Ising)
    β = model.β
    if β == 0.0
        return 2.0
    elseif β == 0.2
        return 2.08450374046259
    elseif β == 0.4
        return 2.4093664345022363
    elseif β == 0.6
        return 3.3539286863974582
    elseif β == 0.8
        return 4.96201030069517
    elseif β == 1.0
        return 7.3916307004743125
    end
end


text:
using BinaryBuilder

name = "licensecheck"
version = v"0.3.1"

sources = [
    GitSource("https://github.com/google/licensecheck",
              "5aa300fd3333d0c6592148249397338023cafcce"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
install_license ${WORKSPACE}/srcdir/licensecheck/LICENSE
cd $WORKSPACE/srcdir/licensecheck/

if [[ "${target}" == *-darwin* ]]; then
    # Create official Apple-blessed `xcrun`, needed by the CGO linker
    cat > /usr/bin/xcrun << EOF
#!/bin/sh
exec "\${@}"
EOF
    chmod +x /usr/bin/xcrun
fi

mkdir clib
cp $WORKSPACE/srcdir/main.go clib/main.go
mkdir -p ${libdir}
CGO_ENABLED=1 go build -buildmode=c-shared -o ${libdir}/licensecheck.${dlext} clib/main.go
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("licensecheck", :licensecheck),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :go], preferred_gcc_version=v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
nt = 1000
seed = 99
m = ABCDQR_simula(nt,seed)
y = m["y"]
u = m["u"]

#######################
print("Testing ABCDQR_em_starting: ")
wi = collect(range(0.1, length=2, stop=3))
Ai,Bi,Ci,Di,Qi,Ri,Si,m1i,P1i = ABCDQR_em_starting(y,u,wi,1,true)
if sum(Ai.^2) > 1e-6
    println("OK")
end


text:
"""
    Fibonacci{T} <: UnivariateAlgorithm

The Fibonacci search algorithm optimize a univariate function by iteratively shrinking the
initial bracket around the minimum. The ratio is computed from the Fibonacci sequence.
"""
struct Fibonacci{T} <: UnivariateAlgorithm
    ϵ::T

    function Fibonacci{T}(ϵ) where {T}
        0 < ϵ < 1 ? new{T}(ϵ) : error("`ϵ` must be in the interval (0, 1)")
    end
end

"""
    Fibonacci([ϵ=0.01])

Initialize Fibonacci search algorithm with optional argument `ϵ` in the interval (0, 1)
to control how close the last iterate is to the previous one. Refer to Kochenderfer and
Wheeler's "Algorithms for Optimization" for more information.
"""
Fibonacci(ϵ=0.01) = Fibonacci{typeof(ϵ)}(ϵ)

function _optimize(f, lower::T, upper::T, alg::Fibonacci;
        reltol, abstol, maxiter) where {T}
    ϵ = T(alg.ϵ)
    s = (1 - √5) / (1 + √5)
    p = 1 / (φ * (1 - s^(maxiter + 1)) / (1 - s^maxiter))
    x = T(p * lower + (1 - p) * upper)
    yx = f(x)
    converged = false
    iter = 1
    while true
        x_tol = reltol * abs(x) + abstol
        if abs(upper - lower) < 2x_tol
            converged = true
            break
        end
        iter == maxiter && break
        iter += 1
        if iter == maxiter
            new_x = T(ϵ * upper + (1 - ϵ) * x)
        else
            new_x = T(p * upper + (1 - p) * lower)
        end
        new_yx = f(new_x)
        if yx > new_yx
            lower, x, yx = x, new_x, new_yx
        else
            lower, upper = new_x, lower
        end
        p = 1 / (φ * (1 - s^(maxiter - iter + 2)) / (1 - s^(maxiter - iter + 1)))
    end
    return Solution(converged, iter, x, yx)
end


text:
include("agent.jl")
include("random_policy.jl")
include("q_based_policy.jl")


text:
export make_bonds
export make_paths
export bond2path
export collect_bonds_paths

function make_bonds(atom_group)
    deleteat!(collect(Iterators.product(first(atom_group), atom_group)), 1)
end

function make_paths(atoms3, indexed_atoms)
    a0, a1, a2, a3 = indexed_atoms[atoms3]
    ps1, pe1, ps2, pe2, ps3, pe3 = fill(0, 6)
    p0, p1, p2, p3 = Tuple(atoms3)
    n = 1
    while true
        b1, b2, b3 = broadcast(x -> Line(a0, x), [a1, a2, a3])
        o12, o13, o23, o21, o31, o32 = map(orientation, [b1, b1, b2, b2, b3, b3], [a2, a3, a3, a1, a1, a2])
        if     o12 * o13 == -1
            if o12 == -1
               ps1, pe1, ps2, pe2, ps3, pe3 = p3, p1, p1, p2, p2, p3
            else
               ps1, pe1, ps2, pe2, ps3, pe3 = p3, p2, p2, p1, p1, p3
            end
            break
        elseif o23 * o21 == -1
            if o23 == -1
               ps1, pe1, ps2, pe2, ps3, pe3 = p1, p2, p2, p3, p3, p1
            else
               ps1, pe1, ps2, pe2, ps3, pe3 = p1, p3, p3, p2, p2, p1
            end
            break
        elseif o31 * o32 == -1
            if o31 == -1
               ps1, pe1, ps2, pe2, ps3, pe3 = p2, p3, p3, p1, p1, p2
            else
               ps1, pe1, ps2, pe2, ps3, pe3 = p2, p1, p1, p3, p3, p2
            end
            break
        else
            vibration = 0.01 * n * length2(Line(a1,a2))
            a0 = Point2D(getx(a0) + round(first(randn(1)))*vibration, gety(a0) + round(first(randn(1)))*vibration)
            n += 1
        end
    end
    paths = [(ps1, p0, pe1), (ps2, p0, pe2), (ps3, p0, pe3)]
end

function bond2path(atom_group::Array{Int64,1}, indexed_atoms)
    nn = length(atom_group) - 1
    if     nn == 3
        bonds = make_bonds(atom_group)
        paths = make_paths(atom_group, indexed_atoms)
    elseif nn == 2
        bonds = make_bonds(atom_group)
        path = (atom_group[2], atom_group[1], atom_group[3])
        paths = [path, reverse(path)]
    elseif nn == 1
        bonds = make_bonds(atom_group)
        paths = [(atom_group[2], atom_group[1], atom_group[2])]
    else   nn == 0
        bonds = []
        paths = []
    end
    return bonds, paths
end

function collect_bonds_paths(atoms_groups, indexed_atoms)
    bonds_paths_collection = map(x -> bond2path(x, indexed_atoms), atoms_groups)
    bonds = vcat(unique!(first.(bonds_paths_collection))...)
    paths = vcat(unique!(last.(bonds_paths_collection))...)
    return bonds, paths
end


text:
# Container data structure (structure-of-arrays style) for DG elements
struct ElementContainer1D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 3}                   # [variables, i, elements]
  u_t::Array{Float64, 3}                 # [variables, i, elements]
  u_tmp2::Array{Float64, 3}              # [variables, i, elements]
  u_tmp3::Array{Float64, 3}              # [variables, i, elements]
  inverse_jacobian::Vector{Float64}      # [elements]
  node_coordinates::Array{Float64, 3}    # [orientation, i, elements]
  surface_flux_values::Array{Float64, 3} # [variables, direction, elements]
  cell_ids::Vector{Int}                  # [elements]
end

function ElementContainer1D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, NVARS, n_nodes, capacity)
  u_t = fill(NaN, NVARS, n_nodes, capacity)
  # u_rungakutta is initialized to non-NaN since it is used directly
  u_tmp2 = fill(0.0, NVARS, n_nodes, capacity)
  u_tmp3 = fill(0.0, NVARS, n_nodes, capacity)
  inverse_jacobian = fill(NaN, capacity)
  node_coordinates = fill(NaN, 1, n_nodes, capacity)
  surface_flux_values = fill(NaN, NVARS, 2 * 1, capacity)
  cell_ids = fill(typemin(Int), capacity)

  elements = ElementContainer1D{NVARS, POLYDEG}(u, u_t, u_tmp2, u_tmp3, inverse_jacobian, node_coordinates,
                                                 surface_flux_values, cell_ids)

  return elements
end

# Return number of elements
nelements(elements::ElementContainer1D) = length(elements.cell_ids)

# Container data structure (structure-of-arrays style) for DG interfaces
struct InterfaceContainer1D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 3}      # [leftright, variables, interfaces]
  neighbor_ids::Matrix{Int} # [leftright, interfaces]
  orientations::Vector{Int} # [interfaces]
end

function InterfaceContainer1D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, 2, NVARS, capacity)
  neighbor_ids = fill(typemin(Int), 2, capacity)
  orientations = fill(typemin(Int), capacity)

  interfaces = InterfaceContainer1D{NVARS, POLYDEG}(u, neighbor_ids, orientations)

  return interfaces
end

# Return number of interfaces
ninterfaces(interfaces::InterfaceContainer1D) = length(interfaces.orientations)

# Container data structure (structure-of-arrays style) for DG boundaries
struct BoundaryContainer1D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 3}                # [leftright, variables, boundaries]
  neighbor_ids::Vector{Int}           # [boundaries]
  orientations::Vector{Int}           # [boundaries]
  neighbor_sides::Vector{Int}         # [boundaries]
  node_coordinates::Array{Float64, 2} # [orientation, elements]
end

function BoundaryContainer1D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, 2, NVARS, capacity)
  neighbor_ids = fill(typemin(Int), capacity)
  orientations = fill(typemin(Int), capacity)
  neighbor_sides = fill(typemin(Int), capacity)
  node_coordinates = fill(NaN, 1, capacity)

  boundaries = BoundaryContainer1D{NVARS, POLYDEG}(u, neighbor_ids, orientations, neighbor_sides,
                                                   node_coordinates)

  return boundaries
end

# Return number of boundaries
nboundaries(boundaries::BoundaryContainer1D) = length(boundaries.orientations)


text:
# Dispatch on the day
struct Day{N} end
Day(n::Int) = Day{n}()

function testdata(D::Day{N}) where {N}
    return load(D, joinpath(@__DIR__, "..", "data", @sprintf("test%02d.txt", N)))
end
function data(D::Day{N}) where {N}
    return load(D, joinpath(@__DIR__, "..", "data", @sprintf("data%02d.txt", N)))
end

test1(D::Day) = solve1(D, testdata(D))
test2(D::Day) = solve2(D, testdata(D))

solve1(D::Day) = solve1(D, data(D))
solve2(D::Day) = solve2(D, data(D))


text:
using AxisArrays
using Base.Test

if VERSION < v"0.6.0-dev"
    @test isempty(detect_ambiguities(AxisArrays, Base, Core))
end

include("core.jl")
include("intervals.jl")
include("indexing.jl")
include("sortedvector.jl")
include("search.jl")
include("combine.jl")

include("readme.jl")

nothing


text:
using RCall, JLD

R"
library(gridExtra)
library(cowplot)
library(tidyverse)
library(latex2exp)
library(viridis)
"

DF = load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_3.jld")
DF_PY_apr = DF["DF_PY_apr"]

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
df_r= $DF_PY_apr
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
print(df_r[1:15,])
df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_beta = rgb(df_r$beta_scaled, 0,0.5,1)
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line(alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_py_apr <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.8, linetype='dotted') +
  xlim(0, 100)+ ylim(0,16)+ labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+
   theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15),legend.position='none')
"

R"p_py_apr
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_7_3.pdf', width= 4, height = 4)
plot(p_py_apr)
dev.off()"

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
# AB1 format
# ===========
#
# Reader of the AB1 file format.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioJulia.jl/blob/master/LICENSE.md

module ABIF

import BioGenerics
include("reader.jl")

end # module


text:
export
    Algorithm,
    marginal,
    joint,
    probability,
    mean

"""
    Algorithm

The supertype of all algorithms.

A standard set of queries is defined for algorithms. Any given subtype of `Algorithm`
will implement a subset of these queries.
"""
abstract type Algorithm end


text:
module NMF
    using Compat
    using StatsBase

    import Compat: view
    import Base: sum!

    export nnmf

    include("common.jl")
    include("utils.jl")

    include("initialization.jl")
    include("multupd.jl")
    include("projals.jl")
    include("alspgrad.jl")

    include("interf.jl")

end # module


text:
"""
MurmurHash3 was written by Austin Appleby, and is placed in the public
domain. The author hereby disclaims copyright to this source code.

This version was translated into Julia by Scott P. Jones
It is licensed under the MIT license

Note - The x86 and x64 versions do _not_ produce the same results, as the
algorithms are optimized for their respective platforms. You can still
compile and run any of them on any platform, but your performance with the
non-native version will be less than optimal.
"""
module MurmurHash3
export mmhash128_a, mmhash128_u, mmhash128_c, mmhash32

u8(val)  = val%UInt8
u32(val) = val%UInt32
u64(val) = val%UInt64

@inline rotl(x::Unsigned, r) = (x << r) | (x >>> (sizeof(typeof(x))*8 - r))

@inline xor33(k::UInt64) = xor(k, k >>> 33)

@inline rotl27(k) = rotl(k, 27)
@inline rotl31(k) = rotl(k, 31)
@inline rotl33(k) = rotl(k, 33)

#-----------------------------------------------------------------------------
# Finalization mix - force all bits of a hash block to avalanche

@inline fmix(k::UInt64) = xor33(xor33(xor33(k) * 0xff51afd7ed558ccd) * 0xc4ceb9fe1a85ec53)

const c1 = 0x87c37b91114253d5
const c2 = 0x4cf5ad432745937f

@inline function mhblock(h1, h2, k1, k2)
    h1 = (rotl27(xor(h1, rotl31(k1 * c1) * c2)) + h2) * 5 + 0x52dce729
    h2 = (rotl31(xor(h2, rotl33(k2 * c2) * c1)) + h1) * 5 + 0x38495ab5
    h1, h2
end

@inline function mhbody(nblocks, pnt, h1, h2)
    for i = 1:nblocks
        h1, h2 = mhblock(h1, h2, unsafe_load(pnt), unsafe_load(pnt + 8))
        pnt += 16
    end
    pnt, h1, h2
end

@inline mhtail1(h1, k1) = xor(h1, rotl31(k1 * c1) * c2)
@inline mhtail2(h2, k2) = xor(h2, rotl33(k2 * c2) * c1)

@inline function mhfin(len, h1, h2)
    h1 = xor(h1, u64(len))
    h2 = xor(h2, u64(len))

    h1 += h2
    h2 += h1

    h1 = fmix(h1)
    h2 = fmix(h2)

    h1 += h2
    h1, h1 + h2
end

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

up8(val)  = u32(val) << 8
up16(val) = u32(val) << 16
up24(val) = u32(val) << 24
up32(val) = u64(val) << 32
up40(val) = u64(val) << 40
up48(val) = u64(val) << 48
up56(val) = u64(val) << 56

dn6(val) = u8(val >>> 6)
dn12(val) = u8(val >>> 12)
dn18(val) = u8(val >>> 18)

msk6(val) = u8(val & 0x3f)

# Support functions for UTF-8 handling
@inline get_utf8_2(ch) = (0xc0 | dn6(ch),  0x80 | msk6(ch))
@inline get_utf8_3(ch) = (0xe0 | dn12(ch), 0x80 | msk6(dn6(ch)), 0x80 | msk6(ch))
@inline get_utf8_4(ch) = (0xf0 | dn18(ch), 0x80 | msk6(dn12(ch)),
                          0x80 | msk6(dn6(ch)), 0x80 | msk6(ch))

# Optimized in-place conversion to UTF-8 for hashing compatibly with isequal / String
@inline shift_n(v, n) = u64(v) << (((n & 7)%UInt)<<3)

# if cnt == 0 - 4, bytes must fit in k1
# cnt between 5 - 8, may overflow into k2
# if l == 8 - 12,  bytes must fit in k2
# cnt between 12 - 15, may overflow into k3

mergebytes(b1, b2)         = b1 | up8(b2)
mergebytes(b1, b2, b3)     = b1 | up8(b2) | up16(b3)
mergebytes(b1, b2, b3, b4) = b1 | up8(b2) | up16(b3) | up24(b4)

@inline function add_utf8(cnt, ch, k1::UInt64)
    if ch <= 0x7f
        cnt + 1, k1 | shift_n(ch, cnt)
    elseif ch <= 0x7ff
        b1, b2 = get_utf8_2(ch)
        cnt + 2, k1 | shift_n(mergebytes(b1, b2), cnt)
    elseif ch <= 0xffff
        b1, b2, b3 = get_utf8_3(ch)
        cnt + 3, k1 | shift_n(mergebytes(b1, b2, b3), cnt)
    else
        b1, b2, b3, b4 = get_utf8_4(ch)
        cnt + 4, k1 | shift_n(mergebytes(b1, b2, b3, b4), cnt)
    end
end

@inline function add_utf8_split(cnt, ch, k1::UInt64)
    if ch <= 0x7f
        cnt + 1, k1 | shift_n(ch, cnt), u64(0)
    elseif ch <= 0x7ff
        b1, b2 = get_utf8_2(ch)
        if (cnt & 7) == 7
            cnt + 2, k1 | up56(b1), u64(b2)
        else
            cnt + 2, k1 | shift_n(b1 | up8(b2), cnt), u64(0)
        end
    elseif ch <= 0xffff
        b1, b2, b3 = get_utf8_3(ch)
        if (cnt & 7) == 5
            cnt + 3, k1 | up40(b1) | up48(b2) | up56(b3), u64(0)
        elseif (cnt & 7) == 6
            cnt + 3, k1 | up48(b1) | up56(b2), u64(b3)
        else
            cnt + 3, k1 | up56(b1), u64(b2) | up8(b3)
        end
    else
        # This will always go over, may be 1, 2, 3 bytes in second word
        b1, b2, b3, b4 = get_utf8_4(ch)
        if (cnt & 7) == 5
            cnt + 4, k1 | up40(b1) | up48(b2) | up56(b3), u64(b4)
        elseif (cnt & 7) == 6
            cnt + 4, k1 | up48(b1) | up56(b2), b3 | up8(b4)
        else
            cnt + 4, k1 | up56(b1), u64(b2) | up8(b3) | up16(b4)
        end
    end
end

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

# AbstractString MurmurHash3, converts to UTF-8 on the fly
function mmhash128_8_c(str::AbstractString, seed::UInt32)
    k1 = k2 = u64(0)
    h1 = h2 = u64(seed)
    cnt = len = 0
    @inbounds for ch in str
        if cnt < 5
            cnt, k1 = add_utf8(cnt, u32(ch), k1)
        elseif cnt < 8
            cnt, k1, k2 = add_utf8_split(cnt, u32(ch), k1)
        elseif cnt < 13
            cnt, k2 = add_utf8(cnt, u32(ch), k2)
        else
            cnt, k2, k3 = add_utf8_split(cnt, u32(ch), k2)
            # When k1 and k2 are full, then hash another block
            if cnt > 15
                h1, h2 = mhblock(h1, h2, k1, k2)
                k1 = k3
                k2 = u64(0)
                len += 16
                cnt &= 15
            end
        end
    end
    # We should now have characters in k1 and k2, and total length in len
    if cnt != 0
        h1 = mhtail1(h1, k1)
        cnt > 8 && (h2 = mhtail2(h2, k2))
    end
    mhfin(len + cnt, h1, h2)
end

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

# Note: this is designed to work on the Str/String types, where I know in advance that
# the start of the strings are 8-byte aligned, and it is safe to access a full
# 8-byte chunk always at the end (simply masking off the remaining 1-7 bytes)

@inline mask_load(pnt, left) = unsafe_load(pnt) & ((UInt64(1) << ((left & 7) << 3)) - 0x1)

function mmhash128_8_a(len::Integer, pnt::Ptr, seed::UInt32)
    pnt8, h1, h2 = mhbody(len >>> 4, reinterpret(Ptr{UInt64}, pnt), u64(seed), u64(seed))
    if (left = len & 15) > 0
        h1 = mhtail1(h1, left < 8 ? mask_load(pnt8, left) : unsafe_load(pnt8))
        left > 8 && (h2 = mhtail2(h2, mask_load(pnt8 + 8, left)))
    end
    mhfin(len, h1, h2)
end

function mmhash128_8_a(seed::Integer)
    h1 = fmix(2 * u64(seed))
    h2 = fmix(3 * u64(seed))
    h1 + h2, h1 + 2 * h2
end

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

# Combine bits from k1, k2, k3
@inline shift_mix(shft, k1::T, k2::T, k3::T) where {T<:Union{UInt32,UInt64}} =
    k1 >>> shft | (k2 << (sizeof(T)*8 - shft)),
    k2 >>> shft | (k3 << (sizeof(T)*8 - shft)),
    k3 >>> shft

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

function mmhash128_8_u(len::Integer, unaligned_pnt::Ptr, seed::UInt32)
    # Should optimize handling of short (< 16 byte) unaligned strings
    ulp = reinterpret(UInt, unaligned_pnt)
    pnt = reinterpret(Ptr{UInt64}, ulp & ~u64(7))
    fin = reinterpret(Ptr{UInt64}, (ulp + len + 0x7) & ~u64(7)) - 8
    shft = (ulp & u64(7))<<3
    # println("_mmhash128_8_u($len, $unaligned_pnt, $seed) => $pnt, $fin")
    h1 = h2 = u64(seed)
    k1 = unsafe_load(pnt) # Pick up first 1-7 bytes
    k2 = u64(0)
    while pnt < fin
        k1, k2, k3 = shift_mix(shft, k1, unsafe_load(pnt += 8), unsafe_load(pnt += 8))
        # print(" pnt=$pnt, k1=0x$(outhex(k1)), k2=0x$(outhex(k2))")
        h1, h2 = mhblock(h1, h2, k1, k2)
        # println(" => h1=0x$(outhex(h1)), h2=0x$(outhex(h2))")
        k1 = k3
    end
    # We should now have characters in k1 and k2, and total length in len
    if (len & 15) != 0
        h1 = mhtail1(h1, k1)
        (len & 15) > 8 && (h2 = mhtail2(h2, k2))
    end
    # println(" len=$len, h1=0x$(outhex(h1)), h2=0x$(outhex(h2))")
    mhfin(len, h1, h2)
end

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

@inline xor16(k::UInt32) = xor(k, k >>> 16)
@inline xor13(k::UInt32) = xor(k, k >>> 13)

# I don't think these help the generated code anymore (but they make the code easier to read)
@inline rotl13(k) = rotl(k, 13)
@inline rotl15(k) = rotl(k, 15)
@inline rotl16(k) = rotl(k, 16)
@inline rotl17(k) = rotl(k, 17)
@inline rotl18(k) = rotl(k, 18)
@inline rotl19(k) = rotl(k, 19)

# Constants for mmhash_32
const d1 = 0xcc9e2d51
const d2 = 0x1b873593

@inline fmix(h::UInt32) = xor16(xor13(xor16(h) * 0x85ebca6b) * 0xc2b2ae35)

@inline mhblock(h1, k1) = rotl13(xor(h1, rotl15(k1 * d1) * d2))*5 + 0xe6546b64

@inline function mhbody(nblocks, pnt, h1)
    for i = 1:nblocks
        h1 = mhblock(h1, unsafe_load(pnt))
        pnt += 4
    end
    pnt, h1
end

function mmhash32(len, pnt, seed::UInt32)
    pnt, h1 = mhbody(len >>> 2, reinterpret(Ptr{UInt32}, pnt), seed)
    (len & 3) == 0 || (h1 = xor(h1, rotl15(unsafe_load(pnt)) * d1) * d2)
    fmix(xor(h1, u32(len)))
end

@inline function mhfin(len, h1, h2, h3, h4)
    h1 = xor(h1, u32(len))
    h2 = xor(h2, u32(len))
    h3 = xor(h3, u32(len))
    h4 = xor(h4, u32(len))

    h1 += h2; h1 += h3; h1 += h4; h2 += h1; h3 += h1; h4 += h1

    h1 = fmix(h1)
    h2 = fmix(h2)
    h3 = fmix(h3)
    h4 = fmix(h4)

    h1 += h2; h1 += h3; h1 += h4; h2 += h1; h3 += h1; h4 += h1

    up32(h2) | h1, up32(h4) | h3
end

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

# Calculate MurmurHash for 32-bit platforms

# Constants for mmhash128_4
const e1 = 0x239b961b
const e2 = 0xab0e9789
const e3 = 0x38b34ae5
const e4 = 0xa1e38b93

@inline function mhblock(h1, h2, h3, h4, k1, k2, k3, k4)
    h1 = (rotl19(xor(h1, rotl15(k1 * e1) * e2)) + h2)*5 + 0x561ccd1b
    h2 = (rotl17(xor(h2, rotl16(k2 * e2) * e3)) + h3)*5 + 0x0bcaa747
    h3 = (rotl15(xor(h3, rotl17(k3 * e3) * e4)) + h4)*5 + 0x96cd1c35
    h4 = (rotl13(xor(h4, rotl18(k4 * e4) * e1)) + h1)*5 + 0x32ac3b17
    h1, h2, h3, h4
end

@inline function mhbody(nblocks, pnt, h1, h2, h3, h4)
    for i = 1:nblocks
        h1, h2, h3, h4 =
            mhblock(h1, h2, h3, h4,
                    unsafe_load(pnt), unsafe_load(pnt+4), unsafe_load(pnt+8), unsafe_load(pnt+12))
        pnt += 16
    end
    pnt, h1, h2, h3, h4
end

function mmhash128_4(len, pnt, seed::UInt32)
    pnt, h1, h2, h3, h4 = mhbody(len >>> 4, pnt, seed, seed, seed, seed)
    if (left = len & 15) != 0
        h1  = xor(h1, rotl16(unsafe_load(pnt) * e1) * e2)
        if left > 4
            h2  = xor(h2, rotl16(unsafe_load(pnt+4) * e2) * e3)
            if left > 8
                h3  = xor(h3, rotl17(unsafe_load(pnt+8) * e3) * e4)
                left > 12 && (h4  = xor(h4, rotl18(unsafe_load(pnt+12) * e4) * e1))
            end
        end
    end
    mhfin(len, h1, h2, h3, h4)
end

@inline shift_n_32(v, n) = u32(v) << (((n & 7)%UInt)<<3)

@inline function get_utf8(cnt, ch)
    if ch <= 0x7f
        cnt + 1, u32(ch)
    elseif ch <= 0x7ff
        b1, b2 = get_utf8_2(ch)
        cnt + 2, mergebytes(b1, b2)
    elseif ch <= 0xffff
        b1, b2, b3 = get_utf8_3(ch)
        cnt + 3, mergebytes(b1, b2, b3)
    else
        b1, b2, b3, b4 = get_utf8_4(ch)
        cnt + 4, mergebytes(b1, b2, b3, b4)
    end
end

@inline function add_utf8_split(cnt, ch, k1::UInt32)
    if ch <= 0x7f
        cnt + 1, k1 | shift_n_32(ch, cnt), u32(0)
    elseif ch <= 0x7ff
        b1, b2 = get_utf8_2(ch)
        if (cnt & 7) == 3
            cnt + 2, k1 | up24(b1), u32(b2)
        else
            cnt + 2, k1 | shift_n_32(b1 | up8(b2), cnt), u32(0)
        end
    elseif ch <= 0xffff
        b1, b2, b3 = get_utf8_3(ch)
        if (cnt & 7) == 1
            cnt + 3, k1 | up8(b1) | up16(b2) | up24(b3), u32(0)
        elseif (cnt & 7) == 2
            cnt + 3, mergebytes(k1, b1, b2), u32(b3)
        else
            cnt + 3, k1 | up24(b1), b2 | up8(b3)
        end
    else
        # This will always go over, may be 1, 2, 3 bytes in second word
        b1, b2, b3, b4 = get_utf8_4(ch)
        if (cnt & 7) == 1
            cnt + 4, mergebytes(k1, b1, b2, b3), u32(b4)
        elseif (cnt & 7) == 2
            cnt + 4, mergebytes(k1, b1, b2), b3 | up8(b4)
        else
            cnt + 4, k1 | up24(b1), mergebytes(b2, b3, b4)
        end
    end
end

const mmhash128_c = @static sizeof(Int) == 8 ? mmhash128_8_c : mmhash128_4
const mmhash128_a = @static sizeof(Int) == 8 ? mmhash128_8_a : mmhash128_4
const mmhash128_u = @static sizeof(Int) == 8 ? mmhash128_8_c : mmhash128_4 # Todo: fix unaligned

end # module MurmurHash3


text:
using Images, ImageView, Statistics

function draw_seam(img, seam)
    img_w_seam = copy(img)
    width = size(img, 1)
    for i = 1: width
        img_w_seam[i, seam[i]] = RGB(1, 1, 1)
    end
    return img_w_seam
end

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
function nqubits(g::Tuple)
  s = g[2]
  n = (s isa Number ? s : maximum(s))
  return n
end

nqubits(gates::Vector{<:Any}) = maximum((nqubits(gate) for gate in gates))

nlayers(circuit::Vector{<:Any}) = 1
nlayers(circuit::Vector{<:Vector{<:Any}}) = length(circuit)

ngates(circuit::Vector{<:Any}) = length(circuit)
ngates(circuit::Vector{<:Vector{<:Any}}) = length(vcat(circuit...))

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
-                               STANDARD CIRCUITS                              -
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    qft(N::Int; inverse::Bool = false)

Generate a list of gates for the quantum fourier transform circuit on `N` sites.
"""
function qft(N::Int; inverse::Bool=false)
  circuit = []
  if inverse
    for j in N:-1:1
      for k in N:-1:(j + 1)
        angle = -π / 2^(k - j)
        push!(circuit, ("CRz", (k, j), (ϕ=angle,)))
      end
      push!(circuit, ("H", j))
    end
  else
    for j in 1:(N - 1)
      push!(circuit, ("H", j))
      for k in (j + 1):N
        angle = π / 2^(k - j)
        push!(circuit, ("CRz", (k, j), (ϕ=angle,)))
      end
    end
    push!(circuit, ("H", N))
  end
  return circuit
end

"""
    ghz(N::Int)

Generate a list of gates for the GHZ state

ψ = (|0,0,…,0⟩ + |1,1,…,1⟩)/√2
"""
function ghz(N::Int)
  #gates = Tuple[("H",1)]
  circuit = []
  push!(circuit, ("H", 1))
  #gates = [[("H",1)...],[("CX",(1,2))...]]
  for j in 1:(N - 1)
    push!(circuit, ("CX", (j, j + 1)))
  end
  return circuit
end

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
-                               LAYER FUNCTIONS                                -
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    gatelayer(gatename::AbstractString, N::Int; kwargs...)

Create a uniform layer containing `N` identical quantum gates, idenfitied by
`gatename`. If additional parameteres are provided, they are identically added to all gates.
"""
function gatelayer(
  gatename::AbstractString, support::Union{Vector{<:Int},AbstractRange}; kwargs...
)
  return [isempty(kwargs) ? (gatename, n) : (gatename, n, values(kwargs)) for n in support]
end

gatelayer(gatename::AbstractString, N::Int; kwargs...) = gatelayer(gatename, 1:N; kwargs...)

"""
    gatelayer(bonds::Vector{Vector{Int}}, gatename::AbstractString)

Create a uniform layer of multi-qubit gates over a set of `bonds`.
"""
function gatelayer(gatename::AbstractString, bonds::Vector{Tuple}; kwargs...)
  return [
    isempty(kwargs) ? (gatename, bonds[n]) : (gatename, bonds[n], values(kwargs)) for
    n in 1:length(bonds)
  ]
end

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
-                               RANDOM CIRCUITS                                -
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    randomlayer(gatename::AbstractString, support::Union{Int,Vector{<:Vector{Int}}}; rng = Random.GLOBAL_RNG)

Generate a random layer built out of a set one or two qubit gates If `support::Int = N`, generates
`N` single-qubit gates `gatename`. If `support::Vector=bonds`, generates a set of two-qubit
gates on the couplings contained in `support`.
"""

function randomlayer(
  gatename::AbstractString,
  support::Union{Vector{<:Int},Vector{<:Tuple},AbstractRange};
  rng=Random.GLOBAL_RNG,
  kwargs...,
)
  layer = []
  for n in support
    pars = randomparams(gatename, length(n); rng = rng) # the 2^n is for the Haar dimension
    gatepars = (isempty(pars) ? (isempty(kwargs) ? nothing : values(kwargs)) : merge(pars,values(kwargs)))
    g = (isnothing(gatepars) ? (gatename, n) : (gatename, n, gatepars))
    push!(layer, g)
  end
  return layer
end

function randomlayer(gatename::AbstractString, support::Int; kwargs...)
  return randomlayer(gatename, 1:support; kwargs...)
end

"""
    randomlayer(gatenames::Vector{<:AbstractString}, support::Union{Int,Vector{<:Vector{Int}}};
                rng = Random.GLOBAL_RNG,
                weights::Union{Nothing,Vector{Float64}} = ones(length(gatenames))/length(gatenames))

Generate a random layer built out of one or two qubit gates, where `gatenames` is a set of possible
gates to choose from. By default, each single gate is sampled uniformaly over this set. If `weights`
are provided, each gate is sampled accordingly.
"""

function randomlayer(
  gatenames::Vector{<:AbstractString},
  support::Union{Vector{<:Int}, AbstractRange, Vector{<:Tuple}};
  rng=Random.GLOBAL_RNG,
  weights::Union{Nothing,Vector{Float64}}=ones(length(gatenames)) / length(gatenames),
  kwargs...,
)
  gate_id = StatsBase.sample(gatenames, StatsBase.Weights(weights), length(support))
  layer = []
  for (i,n) in enumerate(support)
    pars = randomparams(gate_id[i], length(n); rng = rng)
    gatepars = (isempty(pars) ? (isempty(kwargs) ? nothing : values(kwargs)) : merge(pars,values(kwargs)))
    g = (isnothing(gatepars) ? (gate_id[i], n) : (gate_id[i], n, gatepars))
    push!(layer, g)
  end
  return layer
end

function randomlayer(gatenames::Vector{<:AbstractString}, support::Int; kwargs...)
  return randomlayer(gatenames, 1:support; kwargs...)
end

"""
   randomcircuit(N::Int, depth::Int, coupling_sequence::Vector{<:Vector{<:Any}};
                 twoqubitgates::Union{String,Vector{String}} = "Haar",
                 onequbitgates::Union{Nothing,String,Vector{String}} = nothing,
                 layered::Bool = true,
                 rng = Random.GLOBAL_RNG)

Build a random quantum circuit with `N` qubits and depth `depth`.
"""
function randomcircuit(
  coupling_sequence::Vector;
  depth::Int = 1,
  twoqubitgates::Union{String,Vector{String}}="RandomUnitary",
  onequbitgates::Union{Nothing,String,Vector{String}}=nothing,
  layered::Bool=true,
  rng=Random.GLOBAL_RNG,
)
  #N = (coupling_sequence isa Vector{AbstractVector} ? maximum(vcat([maximum.(c) for c in coupling_sequence]...)) :
  #                                            maximum([maximum(c) for c in coupling_sequence]))

  N = 0
  coupling_sequence = coupling_sequence isa Vector{<:Tuple} ? [coupling_sequence] : coupling_sequence
  for seq in coupling_sequence
    for b in seq
      N = max(N,b[1],b[2])
    end
  end
  circuit = Vector[]
  for d in 1:depth
    layer = []
    # two-qubit gates
    bonds = coupling_sequence[(d - 1) % length(coupling_sequence) + 1]
    append!(layer, randomlayer(twoqubitgates, bonds; rng=rng))
    # one-qubit gates
    if !isnothing(onequbitgates)
      append!(layer, randomlayer(onequbitgates, N; rng=rng))
    end
    push!(circuit, layer)
  end
  layered && return circuit
  return vcat(circuit...)
end

"""
    randomcircuit(N::Int, depth::Int; kwargs...)

Generate a 1D random quantum circuit
"""
function randomcircuit(N::Int; kwargs...)
  return randomcircuit(lineararray(N); kwargs...)
end

"""
    randomcircuit(Lx::Int, Ly::Int, depth::Int; rotated::Bool = false, kwargs...)

Generate a 2D random quantum circuit
"""
function randomcircuit(size::Tuple; rotated::Bool=false, kwargs...)
  Lx, Ly = size
  return randomcircuit(squarearray(Lx, Ly; rotated=rotated); kwargs...)
end

function randomcircuit(L::Int, depth::Int; rotated::Bool=false, kwargs...)
  error("randomcircuit(N::Int, depth::Int; kwargs...) is depracated\n
         - for a 1d random circuit: randomcircuit(N::Int; depth = depth, kwargs...)\n
         - for a 2d random circuit: randomcircuit((Lx, Ly); depth = depth, kwargs...)")
end

ITensors.dag(single_gate::Tuple{String,Union{Int,Tuple}}) =
  (single_gate[1], single_gate[2], (dag = true,))

function ITensors.dag(single_gate::Tuple{String,Union{Int,Tuple},NamedTuple})
  prev_dag = get(single_gate[3], :dag, false)
  nt = Base.setindex(single_gate[3], !prev_dag, :dag)
  return (single_gate[1], single_gate[2], nt)
end

ITensors.dag(layer::Vector{<:Any}) =
  [ITensors.dag(g) for g in reverse(layer)]

ITensors.dag(circuit::Vector{<:Vector{<:Any}}) =
  [dag(layer) for layer in reverse(circuit)]


text:
"""
    Computes the volume of a polytope.

# Arguments

`polytope_vertices::AbstractArray{Float64, 2}`    Contains all the intersecting points between all the boundaries of the simplices. Each row represents some intersecting point. Matrix of size D-by-N, where D = total number of intersecting points. Each row is an intersecting point.
`convexp::AbstractArray{Float64, 2}`  Convex expansions of the vertices in `polytope_vertices` in terms of the original simplices. The size of this array is D-by-(2N + 2), where D = total number of intersecting points. (2N + 2) corresponds to the number of vertices in each simplex * 2. The first n+1 columns correspond to the convex expansion coefficients of the intersecting points in terms of the vertices generating simplex 1. The remaining n+1 to (2N+2) columns correspond to the convex expansion coefficients of the intersecting points in terms of the vertices generating simplex 2. The faces of each simplex are numbered according to the column labels.
"""
function VolumeComputation(polytope_vertices::AbstractArray{Float64, 2},
                           convexp::AbstractArray{Float64, 2})

    num_generators = size(polytope_vertices, 1) # The number of polytope generators (vertices)
    dim = size(polytope_vertices, 2)
    intvol = 0.0

    if num_generators == dim + 1 # Intersection is a simplex
        return abs(det(vcat(ones(1, dim + 1), copy(transpose(polytope_vertices)))))
    elseif num_generators > dim + 1 # Intersection is a polytope.
        # The polytope faces and how many of them are nonsingular.
        triangulation_polytope_faces, n_nonsingular_faces = TriangulationPolytopeFaces(convexp, num_generators, dim)

        if n_nonsingular_faces >= dim + 1
            polytope_centroid = ones(1, num_generators) * polytope_vertices / num_generators

            for a = 1:size(triangulation_polytope_faces, 1)
                Ind = round.(Int64, triangulation_polytope_faces[a, :])
                vertices = [polytope_vertices[vec(Ind), :]; polytope_centroid]
                intvol += abs(det([ones(dim + 1, 1) vertices]))
            end
        end
    end
    return intvol
end


text:
"""
Utility functions

Wouter Kouw
28-03-2020
"""

function key_from_value(d::Dict{String, Union{Float64, Symbol}}, k::Union{String, Symbol})
    "Given a value, find its key in a paired string dictionary."
    return collect(keys(d))[findfirst(collect(values(d)) .== k)]
end

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
include("../src/Padding.jl")
include("testing.jl")

importall Padding

for padder = [BytePadding, AnsiX923, Iso10126, Pkcs7, IsoIec7816_4]
    print("Testing ", padder, " ")
    for block_size = 1:10
        for test_value = map(x -> convert(Array{Uint8}, x), ["foo", "foobar", "baz", "this is a test", ""])
            padded = pad!(padder, block_size, copy(test_value))
            @assert_not_equal(test_value, padded)
            @assert_equal(0, length(padded) % block_size)

            unpadded = unpad!(padder, padded)
            @assert_equal(test_value, unpadded)
        end
    end
    println()
end

println()


text:
abstract AbstractBlockMatrix{T} <: AbstractMatrix{T}
typealias AbMat{T} AbstractMatrix{T}

"""Returns the blocks of a matrix (non-recursive)"""
blocks(A::AbstractMatrix) = Any[A]

"""Returns the specified matrix block"""
block(M::AbstractMatrix,inds...)  = getindex(blocks(M),inds...)

"""Returns the sizes of each of the blocks in the matrix (non-recursive)"""
sizes(M::AbstractBlockMatrix) = map(size,blocks(M))

collens(M::AbstractBlockMatrix) = Int64[s[1] for s in sizes(M)[:,1]]
rowlens(M::AbstractBlockMatrix) = Int64[s[2] for s in sizes(M)[1,:]]

"""Returns the size of the block matrix. Assumes row/column independent partitiong."""
function Base.size(M::AbstractBlockMatrix)
  return (sum(collens(M)),sum(rowlens(M)))
end

"""Recursive version of sizes, returning collection of sizes of non-block form matrices in block matrix"""
deepsizes(M::AbstractBlockMatrix) = map(deepsizes,blocks(M))
deepsizes(A::AbstractMatrix) = size(A)

function Base.getindex(M::AbstractBlockMatrix,i,j)
    cinds,rinds = cumsum([0;rowlens(M)]),cumsum([0;collens(M)])
    cblock,rblock = searchsortedfirst(rinds,i)-1,searchsortedfirst(cinds,j)-1
    return block(M,cblock,rblock)[i-cinds[cblock],j-rinds[rblock]]
end

function Base.full(M::AbstractBlockMatrix)
    fullM = similar(M)
    for i = 1:length(M)
        fullM[i] = M[i]
    end
    return fullM
end


text:
function icrgen_qc(winrate, n1, n2;uncommon_icr_probs = [0, 0.85, 0.10, 0.05], rare_icr_probs = [0,0,2/3,1/3])
	first_icr = Dict( (w,uncommon_icr_probs) for w in 0:7)
	second_icr = deepcopy(first_icr)
	second_icr[6] = rare_icr_probs
	second_icr[7] = rare_icr_probs
	third_icr = deepcopy(second_icr)
	for i in 4:5
		third_icr[i] = rare_icr_probs
	end

	avg_first_icr = sum( pdfnegbin_truncated(k; p = winrate, r = 3, maxwins = 7) .* first_icr[k] for k in 0:7 )
	avg_second_icr = sum( pdfnegbin_truncated(k; p = winrate, r = 3, maxwins = 7) .* second_icr[k] for k in 0:7 )
	avg_third_icr = sum( pdfnegbin_truncated(k; p = winrate, r = 3, maxwins = 7) .* third_icr[k] for k in 0:7 )

	if n1 > 0
		icrvec = [[ Categorical(avg_first_icr), Categorical(avg_second_icr), Categorical(avg_third_icr) ] for i in 1:n1]
		for i in 1:n2
			push!(icrvec, Categorical{Float64}[])
		end
		icrcycle = Cycle(  icrvec )
	else
		icrcycle = Cycle( [Categorical{Float64}[]] )
	end

	icrcycle
end

function icrgen_singleton(winrate, n1, n2;uncommon_icr_probs = [0, 0.85, 0.10, 0.05], rare_icr_probs = [0,0,2/3,1/3])
	first_icr = Dict( (w,rare_icr_probs) for w in 0:5)
	second_icr = Dict( (w, uncommon_icr_probs) for w in 0:5)

	avg_first_icr = sum( pdfnegbin_truncated(k; p = winrate, r = 2, maxwins = 5) .* first_icr[k] for k in 0:5 )
	avg_second_icr = sum( pdfnegbin_truncated(k; p = winrate, r = 2, maxwins = 5) .* second_icr[k] for k in 0:5 )
	if n1 > 0
		icrvec = [[ Categorical(avg_first_icr), Categorical(avg_second_icr)] for i in 1:n1]
		for i in 1:n2
			push!(icrvec, Categorical{Float64}[])
		end
		icrcycle = Cycle(  icrvec )
	else
		icrcycle = Cycle( [Categorical{Float64}[]] )
	end

	icrcycle
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
using Base.Test

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

@testset "mirrortypes/predicates.jl" begin
    include("mirrortypes/predicates.jl")
end

@testset "mirrortypes/types.jl" begin
    include("mirrortypes/types.jl")
end

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

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

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

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

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

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


text:
push!(LOAD_PATH,"../src/")
using Documenter, PointCloudRasterizers

makedocs(modules = [PointCloudRasterizers], sitename = "PointCloudRasterizers.jl")

deploydocs(
    repo = "github.com/Deltares/PointCloudRasterizers.jl.git",
)


text:
module Brikhead

using CodecZlib

export load_brikhead, load_brikhead_dir
export write_brikhead

struct AFNIAttribute
    name
    type
    count
    contents
end

# This will let us use hdr.attribute_name instead of hdr[:attribute_name].contents
function Base.getproperty(hdr::Dict{Symbol,AFNIAttribute}, sym::Symbol)
    if sym ∈ propertynames(hdr)
        # Dict has innate properties, this case handles them
        return getfield(hdr, sym)
    else
        return getindex(hdr, sym).contents
    end
end

function parse_afni_head(file)
    chunks = get_chunks(file)
    attributes = parse_chunk.(chunks)
    return pack_attributes_as_dict(attributes)
end

function get_chunks(file)
    f = open(file)
    lines = readlines(f)
    close(f)
    typelines = contains.(lines, r"type\s*=")
    chunks = Array{String,1}[]
    idx = findall(typelines)
    push!(idx, length(lines)+1)
    for i in 2:length(idx)
        push!(chunks, lines[idx[i-1]:idx[i]-1])
    end
    return chunks
end

function parse_chunk(chunk)
    type = parse_type(chunk)
    name = parse_name(chunk)
    count = parse_count(chunk)
    contents = parse_contents(chunk, type)
    return AFNIAttribute(name, type, count, contents)
end

function parse_type(chunk)
    line = first(chunk)
    if contains(line, "string")
        return String
    elseif contains(line, "integer")
        return Int
    elseif contains(line, "float")
        return Float64
    end
    error("No valid attribute found while parsing header in chunk:\n$(chunk)")
end

type2attribute(::Type{String}) = "string-attribute"
type2attribute(::Type{Int}) = "integer-attribute"
type2attribute(::Type{Float64}) = "float-attribute"

function parse_name(chunk)
    nameline = chunk[2]
    name = strip(split(nameline, "=")[2])
    return name
end

function parse_count(chunk)
    countline = chunk[3]
    count = parse(Int, split(countline, "=")[2])
    return count
end

function parse_contents(chunk, ::Type{String})
    # First character is ' and last is ~, so skip em
    return chunk[4][2:end-2]
end

function parse_contents(chunk, T::Type)
    contents_str = clean_chunk_contents(chunk)
    # Julia's builtin parse() takes care of everything
    contents = parse.(T, contents_str)
    return contents
end

function clean_chunk_contents(chunk)
    # Join all the lines to make parsing easier
    str_joined = strip(join(chunk[4:end]))
    # Numerical array are separated by spaces (tabs?)
    str_array = split(str_joined, " ")
    # Previous step converts multiple spaces into empty elements
    filter!(element -> !isempty(element), str_array)
    return str_array
end

function pack_attributes_as_dict(attributes::Vector{AFNIAttribute})
    dict = Dict{Symbol, AFNIAttribute}()
    for attribute in attributes
        key = Symbol(lowercase(attribute.name))
        dict[key] = attribute
    end
    return dict
end

function get_shape(hdr)
    (rank, num_timepoints) = hdr.dataset_rank[1:2]
    volume_dimensions = hdr.dataset_dimensions[1:rank]
    shape = (volume_dimensions..., num_timepoints)
    return shape
end

function get_resolution(hdr)
    voxel_sizes = abs.(hdr.delta)
    if haskey(hdr, :taxis_floats)
        temporal_resolution = hdr.taxis_floats[2]
    else
        temporal_resolution = 0
    end
    resolution = (voxel_sizes..., temporal_resolution)
    return resolution
end

function load_brikhead_dir(path)
		files = readdir(dirname(path); join = true, sort = true)
		brikfiles = filter(is_brik, files)
		return load_brikhead.(brikfiles)
end

function load_brikhead(file)
    brikhead = find_brikhead(file)
    hdr = parse_afni_head(brikhead[:head])
    img = parse_afni_brik(brikhead[:brik], hdr)
    return (; hdr, img)
end

function write_brikhead(file, brik)
    write_afni_head(file, brik)
    write_afni_brik(file, brik)
    return file
end

function write_afni_head(file, brik)
    file *= ".HEAD"
    rm(file; force = true)
    open(file, "w") do io
        for attribute in sort(collect(keys(brik.hdr)))
            write_chunk(io, brik.hdr[attribute])
        end
    end
    return
end

function write_chunk(io, chunk::AFNIAttribute)
    out = [
        "type = $(type2attribute(chunk.type))\n",
        "name = $(chunk.name)\n",
        "count = $(chunk.count)\n",
    ]
    if chunk.type == String
        contents = replace(chunk.contents, "\\n" => "\n")
        push!(out, "'" * contents * "~")
    else
        for (idx, str) in enumerate(string.(chunk.contents))
            push!(out, " $str")
            if rem(idx, 5) == 0
                push!(out, "\n")
            end
        end
    end
    write(io, join(out))
    write(io, "\n\n")
    return
end

function filename_without_ext(filepath)
    filename = basename(filepath)
    filename_body = first(split(filename, "."))
    return filename_body
end

function find_brikhead(file)
    if is_brik(file)
        brik = file
        head = find_head(file)
    elseif is_head(file)
        brik = find_brik(file)
        head = file
    else
        brik = find_brik(file)
        head = find_head(file)
    end
    return Dict(:brik => brik, :head => head)
end

find_head(file) = find_otherfile(file, is_head)
find_brik(file) = find_otherfile(file, is_brik)

function find_otherfile(file, filterfunc)
    dir = dirname(file)
    otherfiles = readdir(dir; join=true)
    filter!(filterfunc, otherfiles)
    for otherfile in otherfiles
        if filename_without_ext(otherfile) == filename_without_ext(file)
            return otherfile
        end
    end
    error("Couldn't find file paired with $(file)")
end

is_head(file) = filename_contains(file, "head")
is_brik(file) = filename_contains(file, "brik")
is_gz(file) = filename_contains(file, "gz")

function filename_contains(file, needle)
    return needle ∈ split(lowercase(basename(file)), ".")
end

const briktypes = Dict(
    0 => UInt8,
    1 => Int16,
    3 => Float32,
    5 => Float64
)

function parse_afni_brik(file, hdr)
    @assert length(unique(hdr.brick_types)) == 1
    vol_shape = get_shape(hdr)
    num_elements = prod(vol_shape)
    briktype = briktypes[hdr.brick_types[1]]
    iobuffer = open(file)
    if is_gz(file)
        iobuffer = GzipDecompressorStream(iobuffer)
    end
    filebytes = read(iobuffer)
    close(iobuffer)
    # ToDo? Account for: d.byteorder_string
    image_volume = reshape(reinterpret(briktype, filebytes), vol_shape)
    return permutedims(image_volume, (2,1,3,4))
end

function write_afni_brik(file, brik)
    file = file * ".BRIK.gz"
    rm(file; force = true)
    hdr = brik.hdr
    @assert length(unique(hdr.brick_types)) == 1
    vol_shape = get_shape(hdr)
    num_elements = prod(vol_shape)
    briktype = briktypes[hdr.brick_types[1]]
    img = permutedims(brik.img, (2,1,3,4))
    bytes = reinterpret(briktype, vec(img))
    iobuffer = open(file, "w")
    if is_gz(file)
        stream = GzipCompressorStream(iobuffer)
        write(stream, bytes)
        close(stream)
    else
        write(iobuffer, bytes)
    end
    close(iobuffer)
    return file
end

end


text:
# for test purposes

include("./domain.jl")

# createDatabaseItem - internal function
const agent = "AB9F"
createDatabaseItem(item::Any; agent=agent, action="CREATE") = DatabaseItem(time(), agent, action, item.id, item)


text:
struct SimulationResults
    global_index::MAPPING_DICT
    bus_lookup::Dict{Int, Int}
    system::PSY.System
    time_log::Dict{Symbol, Any}
    solution::SciMLBase.AbstractODESolution
    function SimulationResults(
        inputs::SimulationInputs,
        system::PSY.System,
        time_log,
        solution,
    )
        new(make_global_state_map(inputs), get_lookup(inputs), system, time_log, solution)
    end
end

get_global_index(res::SimulationResults) = res.global_index
get_bus_count(res::SimulationResults) = get_n_buses(res.system)
get_bus_lookup(res::SimulationResults) = res.bus_lookup
get_system(res::SimulationResults) = res.system
get_solution(res::SimulationResults) = res.solution


text:
"""
    flow_val(b::TapTransformer)

Calculates the From - To comp[lex power flow (Flow injected at the bus) of branch of type
TapTransformer

"""
function flow_val(b::TapTransformer)
    Y_t = 1 / (PowerSystems.get_r(b) + PowerSystems.get_x(b) * 1im)
    c = 1 / PowerSystems.get_tap(b)
    arc = PowerSystems.get_arc(b)
    V_from = arc.from.voltage * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.voltage * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = (V_from * Y_t * c^2) - (V_to * Y_t * c)
    flow = V_from * conj(I)
    return flow
end

"""
    flow_val(b::TapTransformer)

Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
Line

"""
function flow_val(b::Line)
    Y_t = (1 / (PowerSystems.get_r(b) + PowerSystems.get_x(b) * 1im))
    arc = PowerSystems.get_arc(b)
    V_from = arc.from.voltage * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.voltage * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = V_from * (Y_t + (1im * PowerSystems.get_b(b).from)) - V_to * Y_t
    flow = V_from * conj(I)
    return flow
end

"""
    flow_val(b::TapTransformer)

Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
Transformer2W

"""
function flow_val(b::Transformer2W)
    Y_t = 1 / (PowerSystems.get_r(b) + PowerSystems.get_x(b) * 1im)
    arc = PowerSystems.get_arc(b)
    V_from = arc.from.voltage * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.voltage * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = V_from * (Y_t + (1im * PowerSystems.get_primaryshunt(b))) - V_to * Y_t
    flow = V_from * conj(I)
    return flow
end

function flow_val(b::PhaseShiftingTransformer)
    error("Systems with PhaseShiftingTransformer not supported yet")
    return
end

function _update_branch_flow!(sys::System)
    for b in get_components(ACBranch, sys)
        S_flow = flow_val(b)
        b.activepower_flow = real(S_flow)
        b.reactivepower_flow = imag(S_flow)
    end
end

function _write_pf_sol!(sys::System, nl_result)
    result = round.(nl_result.zero; digits = 7)
    buses = enumerate(sort(collect(get_components(Bus, sys)), by = x -> get_number(x)))

    for (ix, bus) in buses
        if bus.bustype == BusTypes.REF
            P_gen = result[2 * ix - 1]
            Q_gen = result[2 * ix]
            injection = get_components(StaticInjection, sys)
            devices = [d for d in injection if d.bus == bus && !isa(d, ElectricLoad)]
            generator = devices[1]
            generator.activepower = P_gen
            generator.reactivepower = Q_gen
        elseif bus.bustype == BusTypes.PV
            Q_gen = result[2 * ix - 1]
            θ = result[2 * ix]
            injection_components = get_components(Generator, sys)
            devices = [d for d in injection_components if d.bus == bus]
            if length(devices) == 1
                generator = devices[1]
                generator.reactivepower = Q_gen
            end
            bus.angle = θ
        elseif bus.bustype == BusTypes.PQ
            Vm = result[2 * ix - 1]
            θ = result[2 * ix]
            bus.voltage = Vm
            bus.angle = θ
        end
    end

    _update_branch_flow!(sys)

    return
end

# TODO: Apply actions according to load type
function _get_load_data(sys::System, b::Bus)
    activepower = 0.0
    reactivepower = 0.0
    for l in get_components(ElectricLoad, sys)
        if (l.bus == b)
            if isa(l, FixedAdmittance)
                # Assume v rated = 1.0
                activepower += real(get_Y(l))
                reactivepower += imag(get_Y(l))
            else
                activepower += get_activepower(l)
                reactivepower += get_reactivepower(l)
            end
        end
    end
    return activepower, reactivepower
end

"""
    solve_powerflow!(system, solve_function, args...)

Solves a the power flow into the system and writes the solution into the relevant structs.
Updates generators active and reactive power setpoints and branches active and reactive
power flows (calculated in the From - To direction) (see
[`flow_val`](@ref))

Requires loading NLsolve.jl to run.

Supports passing NLsolve kwargs in the args. By default shows the solver trace.

Arguments available for `nlsolve`:

* `method` : See NLSolve.jl documentation for available solvers
* `xtol`: norm difference in `x` between two successive iterates under which
  convergence is declared. Default: `0.0`.
* `ftol`: infinite norm of residuals under which convergence is declared.
  Default: `1e-8`.
* `iterations`: maximum number of iterations. Default: `1_000`.
* `store_trace`: should a trace of the optimization algorithm's state be
  stored? Default: `false`.
* `show_trace`: should a trace of the optimization algorithm's state be shown
  on `STDOUT`? Default: `false`.
* `extended_trace`: should additifonal algorithm internals be added to the state
  trace? Default: `false`.

## Examples
```julia
using NLsolve
solve_powerflow!(sys, nlsolve)
# Passing NLsolve arguments
solve_powerflow!(sys, nlsolve, method = :newton)

```

"""
function solve_powerflow!(system::System, nlsolve; args...)
    buses = sort(collect(get_components(Bus, system)), by = x -> get_number(x))
    N_BUS = length(buses)

    # assumes the ordering in YBus is the same as in the buses.
    tempYb = Ybus(get_components(ACBranch, system), buses)
    a = collect(1:N_BUS)
    Yb = Ybus(tempYb.data, (a, a), (_make_ax_ref(a), _make_ax_ref(a))).data
    I, J, V = SparseArrays.findnz(Yb)
    neighbors = [Set{Int}([i]) for i in 1:N_BUS]
    for nz in eachindex(V)
        push!(neighbors[I[nz]], J[nz])
        push!(neighbors[J[nz]], I[nz])
    end
    x0 = zeros(N_BUS * 2)

    # Use vectors to cache data for closure
    # These should be read only
    P_GEN_BUS = fill(0.0, N_BUS)
    Q_GEN_BUS = fill(0.0, N_BUS)
    P_LOAD_BUS = fill(0.0, N_BUS)
    Q_LOAD_BUS = fill(0.0, N_BUS)
    P_net = fill(0.0, N_BUS)
    Q_net = fill(0.0, N_BUS)
    Vm = fill(0.0, N_BUS)
    θ = fill(0.0, N_BUS)

    state_variable_count = 1
    sources = get_components(StaticInjection, system, d -> !isa(d, ElectricLoad))

    for (ix, b) in enumerate(buses)
        bus_number = get_number(b)::Int
        bus_angle = get_angle(b)::Float64
        bus_voltage = get_voltage(b)::Float64
        P_GEN_BUS[ix] = 0.0
        Q_GEN_BUS[ix] = 0.0
        get_ext(b)["neighbors"] = neighbors[ix]
        for gen in sources
            if gen.bus == b
                P_GEN_BUS[ix] += get_activepower(gen)
                Q_GEN_BUS[ix] += get_reactivepower(gen)
            end
        end

        P_LOAD_BUS[ix], Q_LOAD_BUS[ix] = _get_load_data(system, b)

        if b.bustype == BusTypes.REF
            injection_components = get_components(StaticInjection, system, d -> d.bus == b)
            isempty(injection_components) &&
                throw(IS.ConflictingInputsError("The slack bus does not have any injection component. Power Flow can not proceed"))
            Vm[ix] = get_voltage(b)::Float64
            θ[ix] = get_angle(b)::Float64
            x0[state_variable_count] = P_GEN_BUS[ix]
            x0[state_variable_count + 1] = Q_GEN_BUS[ix]
            state_variable_count += 2
        elseif b.bustype == BusTypes.PV
            Vm[ix] = get_voltage(b)::Float64
            x0[state_variable_count] = Q_GEN_BUS[ix]
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        elseif b.bustype == BusTypes.PQ
            x0[state_variable_count] = bus_voltage
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        else
            throw(ArgumentError("Bustype not recognized"))
        end
    end

    @assert state_variable_count - 1 == N_BUS * 2
    bus_types = get_bustype.(buses)

    function pf!(F::Vector{Float64}, X::Vector{Float64})
        for (ix, b) in enumerate(bus_types)
            if b == BusTypes.REF
                # When bustype == REFERENCE Bus, state variables are Active and Reactive Power Generated
                P_net[ix] = X[2 * ix - 1] - P_LOAD_BUS[ix]
                Q_net[ix] = X[2 * ix] - Q_LOAD_BUS[ix]
            elseif b == BusTypes.PV
                # When bustype == PV Bus, state variables are Reactive Power Generated and Voltage Angle
                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
                Q_net[ix] = X[2 * ix - 1] - Q_LOAD_BUS[ix]
                θ[ix] = X[2 * ix]
            elseif b == BusTypes.PQ
                # When bustype == PQ Bus, state variables are Voltage Magnitude and Voltage Angle
                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
                Q_net[ix] = Q_GEN_BUS[ix] - Q_LOAD_BUS[ix]
                Vm[ix] = X[2 * ix - 1]
                θ[ix] = X[2 * ix]
            end
        end

        # F is active and reactive power balance equations at all buses
        for ix_f in a
            S_re = -P_net[ix_f]
            S_im = -Q_net[ix_f]
            for ix_t in neighbors[ix_f]
                gb = real(Yb[ix_f, ix_t])
                bb = imag(Yb[ix_f, ix_t])
                S_re +=
                    Vm[ix_f] *
                    Vm[ix_t] *
                    (gb * cos(θ[ix_f] - θ[ix_t]) + bb * sin(θ[ix_f] - θ[ix_t]))
                S_im +=
                    Vm[ix_f] *
                    Vm[ix_t] *
                    (gb * sin(θ[ix_f] - θ[ix_t]) - bb * cos(θ[ix_f] - θ[ix_t]))
            end
            F[2 * ix_f - 1] = S_re
            F[2 * ix_f] = S_im
        end
    end

    res = nlsolve(pf!, x0; args...)
    @info(res)
    if res.f_converged
        PowerSystems._write_pf_sol!(system, res)
        @info("PowerFlow solve converged, the results have been stored in the system")
        return res.f_converged
    end
    @error("The powerflow solver returned convergence = $(res.f_converged)")
    return res.f_converged
end


text:
# This example cross-predicts a field U from a field V.
# Both fields have to be represented as vectors of matrices.
# Where the fields come from does not matter, but to make the example
# runnable we load one of the test systems of TimeseriesPrediction.
#
# This example uses cubic shell embedding and a linear Barkley model.
#
# Importantly, the results are compared with the "real" evolution of the
# system.

# ### Simulate a test system
using PyPlot
using TimeseriesPrediction

testdir = dirname(dirname(pathof(TimeseriesPrediction)))*"/test"
@assert isdir(testdir)
include(testdir*"/system_defs.jl")

Ttrain = 500
Ttest = 10
T = Ttrain + Ttest

U, V = barkley(T;tskip=100, ssize=(50,50))
summary(U)

# ### Cross predict field U from field V
γ = 5
τ = 1
B = 1
k = 1
bc = PeriodicBoundary()

source_train = V[1: Ttrain]
target_train = U[1: Ttrain]
source_pred  = V[Ttrain  - γ*τ + 1:  T]
target_test  = U[Ttrain        + 1:  T]

em = cubic_shell_embedding(source_train, γ,τ,B,k,bc)
pcaem = PCAEmbedding(source_train, em; maxoutdim=5) # PCA speeds things up!

@time target_pred = crossprediction(source_train, target_train, source_pred, em;
progress = false)

err = [abs.(target_test[i]-target_pred[i]) for i=1:Ttest]

println("Maximum error: ", maximum(maximum(e) for e in err))

# ### Plot prediction

# Deduce field extremal values
source_max = maximum(maximum(s) for s in source_pred)
target_max = max(maximum(maximum(s) for s in target_test),
                 maximum(maximum(s) for s in target_pred))
source_min = minimum(minimum(s) for s in source_pred)
target_min = min(minimum(minimum(s) for s in target_test),
                 minimum(minimum(s) for s in target_pred))

# Plot various predicted frames (only the last one shown here)
for i in [1, length(err)÷2, length(err)]

    fig = figure(figsize=(10,10))
    ax1 = subplot2grid((2,2), (0,0))
    ax2 = subplot2grid((2,2), (0,1))
    ax3 = subplot2grid((2,2), (1,0))
    ax4 = subplot2grid((2,2), (1,1))
    im1 = ax1[:imshow](source_pred[i], cmap="viridis", vmin = source_min, vmax = source_max)
    im2 = ax2[:imshow](target_test[i], cmap="cividis", vmin = target_min, vmax = target_max)
    im3 = ax3[:imshow](target_pred[i], cmap="cividis", vmin = target_min, vmax = target_max)
    im4 = ax4[:imshow](err[i], cmap="inferno", vmin = 0, vmax = target_max - target_min)
    for (im, ax) in zip([im1,im2,im3,im4], [ax1,ax2,ax3,ax4])
        ax[:get_xaxis]()[:set_ticks]([])
        ax[:get_yaxis]()[:set_ticks]([])
        colorbar(im, ax = ax, fraction=0.046, pad=0.04)#, format="%.1f")
    end
    ax1[:set_title]("Source")
    ax2[:set_title]("Target Test")
    ax3[:set_title]("Target Cross-Pred.")
    ax4[:set_title]("absolute error")
    tight_layout(w_pad=0.6, h_pad=0.00001)
    suptitle("frame $i")
end
#md savefig("barkley_cross.png"); nothing # hide
#md # ![barkley_cross](barkley_cross.png)


text:
## ODEs of three species in a simple food chain

#### Define 3 species network
# i(rows) consumes resource j(columns), with i,j =( 1,2,3) corresponding to (basal, intermediate, top) species

A = [0 0 0; 1 0 0; 0 1 0]

# ## Trophic levels
# tl = [1, 2, 3]

# # Z is consumer-resource body-mass ratio
# Z = 100

# # basal species mass
# m = 0.01

#### Parameters
# # write a function to generate mass and foodweb
# function Network(m = 0.01, tl = [1, 2, 3]; Z)
#     # Body mass calculated from Z, trophic levels and basal mass
#     M = m * (Z .^ (tl .- 1))
#     FW = (Z = Z, M = M)
#     return FW
# end

# Z, M = Network(Z = 10)

# write a function to generate just body masses
function BodyMasses(m = 0.01, tl = [1, 2, 3]; Z)
    # Body mass calculated from Z, trophic levels and basal mass
    M = m .* (Z .^ (tl .- 1))
    return M
end

max_stars_repo_name: hanamayall/simple-feeding-DR

max_issues_repo_name: hanamayall/simple-feeding-DR

max_forks_repo_name: hanamayall/simple-feeding-DR


text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(DualComplex(0, 0.0))
end

@test begin
    !iszerodivisor(DualComplex(1,2,3,4))
end

@test begin
    iszerodivisor(DualComplex(0,0,3,4))
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(DualComplex(a))
end

@test begin
    a = 1
    !isreal(DualComplex(a, a))
end

@test begin
    z = random(DualComplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(DualComplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, DualComplex(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, H: 2, A: 3, HA: 4]"
    l == r
end

@test begin
    x = random(DualComplex{BigFloat})
    y = random(DualComplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(DualComplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(DualComplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Cmplex{BigInt})
    y = random(DualComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(Cmplex{BigInt})
    y = random(DualComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Cmplex{BigInt})
    y = random(DualComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    barA(barA(x)) == x
end

@test begin
    x = random(DualComplex{BigInt})
    barH(barH(x)) == x
end

@test begin
    x = random(DualComplex{BigInt})
    barH(barA(x)) == barA(barH(x))
end

@test begin
    x = random(DualComplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = barA(x * y)
    r = barA(x) * barA(y)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = barA(x + y)
    r = barA(x) + barA(y)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = barA(x - y)
    r = barA(x) - barA(y)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = barH(x * y)
    r = barH(x) * barH(y)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = barH(x + y)
    r = barH(x) + barH(y)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = barH(x - y)
    r = barH(x) - barH(y)
    l == r
end

@test begin
    x = random(DualComplex{Rational{BigInt}})
    y = random(DualComplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(DualComplex{Rational{BigInt}})
    y = random(DualComplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(DualComplex{Rational{BigInt}})
    y = random(DualComplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(DualComplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    n = norm(x)
    n > big(0) || n == big(0)
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = norm(x * y)
    r = norm(x) * norm(y)
    l == r
end

@test_throws ErrorException begin
    inv(DualComplex(0, 0.0))
end

@test_throws ErrorException begin
    random(DualComplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(DualComplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(DualComplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(DualComplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(DualComplex{BigInt})
    z == z + zero(DualComplex{BigInt})
end

@test begin
    z = random(DualComplex{BigInt})
    z == z * one(DualComplex{BigInt})
end

@test begin
    z = random(DualComplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(DualComplex{BigInt})
    z == z * one(z)
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = commutator(x, y)
    r = zero(DualComplex{BigInt})
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = alternatorL(x, y)
    r = zero(DualComplex{BigInt})
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = alternatorR(x, y)
    r = zero(DualComplex{BigInt})
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    l = flexator(x, y)
    r = zero(DualComplex{BigInt})
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = associator(x, y, z)
    r = zero(DualComplex{BigInt})
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(DualComplex{BigInt})
    z = random(DualComplex{BigInt})
    l = jacobiator(x, y, z)
    r = 3 * x * y * z
    l == r
end

@test begin
    l = DualComplex(Cmplex(1,2))
    r = DualComplex(1,2,0)
    l == r
end

@test begin
    l = DualComplex(1, 2, 3.0)
    r = DualComplex(1.0, 2.0, 3.0)
    l == r
end

@test begin
    l = DualComplex(1, 2, 3, 4.0)
    r = DualComplex(1.0, 2.0, 3.0, 4.0)
    l == r
end