text:
"""
Created on Wed 03 June 2020

@author: Yoann Pradat

Test rank by stability as implement in SigProfiler matlab package. Compare VNMF results to that of the package on the
same matrices.
"""

#### function for plotting factors
py"""
import matplotlib.cm as cm
import sys
if "." not in sys.path:
    sys.path.append(".")

from demo.src.util.plot_VWH import *

def make_plot_W(W_avg, W_std, stab_avg, stab_std, row_types, row_subtypes, filepath):
    row_groups_colors = {
        "C>A": "cornflowerblue",
        "C>G": "lightseagreen",
        "C>T": "coral",
        "T>A": "sandybrown",
        "T>C": "forestgreen",
        "T>G": "deeppink"
    }
    hmap_boundaries = [0, 0.001, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5]
    hmap_cmap = "Blues"

    fig = plot_W_stab(
        W_avg             = W_avg,
        W_std             = W_std,
        stab_avg          = stab_avg,
        stab_std          = stab_std,
        row_groups        = row_types,
        row_groups_colors = row_groups_colors,
        row_labels        = row_subtypes,
        hmap_boundaries   = hmap_boundaries,
        hmap_cmap         = hmap_cmap,
        figsize           = (16,10),
    )
    fig.savefig(filepath, bbox_inches="tight")
"""

#### nmf algorithm with identical parameters
nmf_global_params = NMFParams(
    max_iter      = 10_000,
    dist          = Uniform(0,1),
    stopping_crit = :conn,
    stopping_iter = 100,
    verbose       = false,
    seed          = 123
)

nmf_local_params = FIParams(
    β            = 1,
    scale_W_iter = false,
    scale_W_last = true,
    alg          = :mm
)

nmf_alg = VNMF._nmf_FI

#### read row types and row subtypes
row_types    = vec(readdlm(pwd() * "/tests/data/breast_21/types.txt", ',', String, '\n'))
row_subtypes = vec(readdlm(pwd() * "/tests/data/breast_21/subtypes.txt", ',', String, '\n'))
row_types_order = sortperm(row_types)
V_breast_21 = readdlm(pwd() * "/tests/data/breast_21/originalGenomes.txt", ',', Float64, '\n')

K_min  = 1
K_max  = 5
n_iter = 20

#### # 1. USE PERTURBED AND INIT MATRICES PRODUCED BY SIGPROFILER
#### # ################################################

folder = pwd() * "/tests/plot/breast_21/"
mkpath(folder)

rng = MersenneTwister(123)

@testset "breast_21" begin
    #### loop over candidate ranks
    for K = K_min:K_max
        nmf_global_params.rank =  K

        W_init = readdlm(pwd() * "/tests/data/breast_21/W_init_rank_$(K).txt", ',', Float64, '\n')
        H_init = readdlm(pwd() * "/tests/data/breast_21/H_init_rank_$(K).txt", ',', Float64, '\n')

        W_fit = readdlm(pwd() * "/tests/data/breast_21/W_fit_rank_$(K).txt", ',', Float64, '\n')
        H_fit = readdlm(pwd() * "/tests/data/breast_21/H_fit_rank_$(K).txt", ',', Float64, '\n')

        W_fit_julia = zeros(Float64, size(W_fit))
        H_fit_julia = zeros(Float64, size(H_fit))

        list_nmf_results = Array{NMFResults{Float64},1}()

        #### loop over bootstrap iterations
        for j = 1:n_iter
            v_pert = readdlm(pwd() * "/tests/data/breast_21/V_iter_$(j)_rank_$(K).txt", ',', Float64, '\n')

            w_init = reshape(W_init[:,((j-1)*K+1):j*K], size(W_init,1), K)
            h_init = reshape(H_init[((j-1)*K+1):j*K,:], K, size(H_init,2))

            nmf_results = nmf_alg(v_pert, w_init, h_init, nmf_global_params, nmf_local_params)
            push!(list_nmf_results, nmf_results)

            W_fit_julia[:, ((j-1)*K+1):j*K] .= nmf_results.W
            H_fit_julia[((j-1)*K+1):j*K, :] .= nmf_results.H
        end

        @test maximum(abs.(W_fit_julia - W_fit)) ≈ 0 atol = 1e-8
        @test maximum(abs.(H_fit_julia - H_fit)) ≈ 0 atol = 1e-8

        #### cluster factors across iterations
        clu_results = VNMF.evaluate_stability(list_nmf_results, dist_func = CosineDist(), rng=rng)

        filepath = joinpath(folder, "plot_W_rank_$(K)_sigpro.pdf")
        py"make_plot_W"(
            clu_results.W_avg[row_types_order,:],
            clu_results.W_std[row_types_order,:],
            clu_results.stab_avg,
            clu_results.stab_std,
            row_types[row_types_order],
            row_subtypes[row_types_order],
            filepath
           )
    end
end

#### # 2. USE PERTURBED MATRICES PRODUCED BY VNMF
#### # ################################################

#### params for rank selection
rs_params = RSParams(
    K_min     = K_min,
    K_max     = K_max,
    n_iter    = 200,
    pert_meth = :multinomial,
    seed      = 123
)

#### nmf algorithm
nmf_global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0, 1),
    max_iter      = 10_000,
    stopping_crit = :conn,
    stopping_iter = 100,
    verbose       = false,
)

nmf_local_params = FIParams(
    β            = 1,
    scale_W_iter = false,
    scale_W_last = true,
    alg          = :mm
)

nmf = NMF(
    solver        = nmf_FI,
    global_params = nmf_global_params,
    local_params  = nmf_local_params
)

#### run rank selection by stability procedure
@time rs_results = rank_by_stability(V_breast_21, rs_params, nmf)

folder = pwd() * "/tests/plot/breast_21/"

for clu_results in rs_results.list_clu_results
    K = length(clu_results.stab_avg)

    filepath = joinpath(folder, "plot_W_rank_$(K)_VNMF.pdf")
    py"make_plot_W"(
        clu_results.W_avg[row_types_order,:],
        clu_results.W_std[row_types_order,:],
        clu_results.stab_avg,
        clu_results.stab_std,
        row_types[row_types_order],
        row_subtypes[row_types_order],
        filepath
    )
end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ
using PrettyPrinting

X, y = @load_iris
@load DecisionTreeClassifier

dtc = DecisionTreeClassifier()
r   = range(dtc, :max_depth, lower=1, upper=5)

tm = TunedModel(model=dtc, ranges=[r, ], measure=cross_entropy)

m = machine(tm, X, y)
fit!(m)

fitted_params(m).best_model.max_depth

tm = TunedModel(model=dtc, ranges=r, operation=predict_mode,
                measure=misclassification_rate)
m = machine(tm, X, y)
fit!(m)
fitted_params(m).best_model.max_depth

r = report(m)
r.best_result

using PyPlot

figure(figsize=(8,6))
res = r.plotting # contains all you need for plotting
plot(res.parameter_values, res.measurements, ls="none", marker="o")

xticks(1:5, fontsize=12)
yticks(fontsize=12)
xlabel("Maximum depth", fontsize=14)
ylabel("Misclassification rate", fontsize=14)
ylim([0, 1])

X = (x1=rand(100), x2=rand(100), x3=rand(100))
y = 2X.x1 - X.x2 + 0.05 * randn(100);

dtr = @load DecisionTreeRegressor
forest = EnsembleModel(atom=dtr)

params(forest) |> pprint

r1 = range(forest, :(atom.n_subfeatures), lower=1, upper=3)
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0)
tm = TunedModel(model=forest, tuning=Grid(resolution=12),
                resampling=CV(nfolds=6), ranges=[r1, r2],
                measure=rms)
m = machine(tm, X, y)
fit!(m);

r = report(m)
r.best_result

figure(figsize=(8,6))

res = r.plotting

vals_sf = res.parameter_values[:, 1]
vals_bf = res.parameter_values[:, 2]

tricontourf(vals_sf, vals_bf, res.measurements)
xlabel("Number of sub-features", fontsize=14)
ylabel("Bagging fraction", fontsize=14)
xticks([1, 2, 3], fontsize=12)
yticks(fontsize=12)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

ext: jl


text:
using LinearAlgebraicRepresentation
using Plasm, SparseArrays
Lar = LinearAlgebraicRepresentation
using Debugger

V = [0.0 0.0 0.0; 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0; 0.5 0.25 0.25; 0.25 0.5 0.25; 0.25 0.25 0.5; 0.25 0.25 0.25; 1.25 0.25 0.25; 0.25 1.25 0.25; 0.25 0.25 1.25]
ld_bounds = ([1, 2, 3, 1, 4, 5, 2, 4, 6, 3, 5, 6, 7, 8, 9, 7, 8, 9, 7, 10, 12, 7, 11, 13, 14, 8, 10, 15, 8, 11, 16, 17, 9, 12, 15, 9, 13, 16, 18, 14, 17, 18], [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12], Int8[1, -1, 1, 1, -1, 1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1])
angles = Array{Int64,1}[[2, 1], [1, 3], [4, 1], [3, 2], [2, 4], [4, 3], [7, 4, 5, 6], [9, 8, 4, 5], [11, 4, 5, 10], [8, 6], [9, 7], [10, 6], [7, 11], [7, 12], [8, 10], [11, 9], [9, 12], [12, 11]]

sigma =1

function get_seed_cell(ld_cellsnum, count_marks)
	s = -1
	for i in 1:ld_cellsnum
		if count_marks[i] == 0
			return i
		elseif count_marks[i] == 1 && s < 0
			s = i
		end
	end
	return s
end

function nextprev(lld::Int64, ld::Int64, norp)
	as = angles[lld]
	#ne = findfirst(as, ld)  (findfirst(isequal(v), A), 0)[1]
	ne = (findfirst(isequal(ld), as), 0)[1]
	while true
		ne += norp
		if ne > length(as)
			ne = 1
		elseif ne < 1
			ne = length(as)
		end

		if count_marks[as[ne]] < 2
			break
		end
	end
	as[ne]
end

function _minimal_cycles(V::Lar.Points,ld_bounds::Lar.ChainOp)
#@show V
#@show ld_bounds
    lld_cellsnum, ld_cellsnum = size(ld_bounds)
    count_marks = zeros(Int8, ld_cellsnum)
    dir_marks = zeros(Int8, ld_cellsnum)
    d_bounds = spzeros(Int8, ld_cellsnum, 0)

    angles = Array{Array{Int64, 1}, 1}(undef, lld_cellsnum)

    for lld in 1:lld_cellsnum
        as = []
        for ld in ld_bounds[lld, :].nzind
            push!(as, (ld, angles_fn(lld, ld)))
        end
        sort!(as, lt=(a,b)->a[2]<b[2])
        as = map(a->a[1], as)
        angles[lld] = as
    end

    while (sigma = get_seed_cell(ld_cellsnum, count_marks)) > 0
        if verbose
            print(Int(floor(50 * sum(count_marks) / ld_cellsnum)), "%\r")
        end

        c_ld = spzeros(Int8, ld_cellsnum)
        if count_marks[sigma] == 0
            c_ld[sigma] = 1
        else
            c_ld[sigma] = -dir_marks[sigma]
        end
        c_lld = ld_bounds*c_ld
        while c_lld.nzind != []
			#corolla = spzeros(Int8, ld_cellsnum)
			corolla = zeros(Int64, ld_cellsnum)
            for tau in c_lld.nzind
                b_ld = ld_bounds[tau, :]
                pivot = intersect(c_ld.nzind, b_ld.nzind)[1]
                adj = nextprev(tau, pivot, sign(-c_lld[tau]))
                corolla[adj] = c_ld[pivot]
                if b_ld[adj] == b_ld[pivot]
                    corolla[adj] *= -1
                end
@show corolla
            end
            c_ld += corolla
            c_lld = ld_bounds*c_ld
@show c_lld
        end
        map(s->count_marks[s] += 1, c_ld.nzind)
        map(s->dir_marks[s] = c_ld[s], c_ld.nzind)
        d_bounds = [d_bounds c_ld]
@show d_bounds
    end
    return d_bounds
end

function minimal_cycles(angles_fn::Function, verbose=false)
    return _minimal_cycles
end

EF = SparseArrays.sparse(ld_bounds)


text:
@testset "QueryStrategies" begin

    dummy_data, _ = load_data("$(@__DIR__)/../../example/dummy.csv")

    @testset "multi kde" begin
        data = hcat(fill([1, 2, 3, 4,], 3)...)
        @test_throws KDEException multi_kde(data)
        data2 = hcat(data, hcat(fill([1, 2, 3, 4,], 2)...))
        @test_throws KDEException multi_kde(data2)
        data3 = hcat(data2, rand(4,20) * 4)
        @test isa(OneClassActiveLearning.multi_kde(data3), PyObject)
    end

    params = Dict{Symbol, Any}()
    @testset "initialize_qs" begin
        @testset "Standard" begin
            qs_types = [RandomQs]
            qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data), qs_types, params)
            for qs in qs_objs
                for labels in [labelmap(fill(:U, 10)), labelmap(fill(:Lin, 10)), labelmap(fill(:Lout, 10))]
                    scores = qs_score(qs, dummy_data, labels)
                    @test length(scores) == size(dummy_data, 2)
                end
            end
        end

        qs_types = [MinimumMarginQs, MinimumLossQs]
        for x in qs_types
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict{Symbol, Any}())
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => nothing))
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => 1.5))
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => -1))
        end
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => 0.05)), qs_types)
        for qs in qs_objs
            @testset "DataBasedQs $(typeof(qs))" begin
                @testset "only :U" begin
                        @test isdefined(qs, :bw_method)
                        @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:U, 10)))
                end

                @testset "only :Lin" begin
                    # number of :Lin must be sufficient that leave-one-out does not result in singular covariance matrix
                    scores = qs_score(qs, dummy_data, labelmap(fill(:Lin, 15)))
                    @test length(scores) == size(dummy_data, 2)
                end
                @testset "only :Lout" begin
                    @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:Lout, 10)))
                end
            end
        end

        qs_types = [ExpectedMinimumMarginQs, ExpectedMaximumEntropyQs]
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, params), qs_types)
        for qs in qs_objs
            @testset "DataBasedQs $(typeof(qs))" begin
                @testset "only :U" begin
                        @test isdefined(qs, :bw_method)
                        @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:U, 10)))
                end

                @testset "only :Lin" begin
                    # number of :Lin must be sufficient that leave-one-out does not result in singular covariance matrix
                    scores = qs_score(qs, dummy_data, labelmap(fill(:Lin, 15)))
                    @test length(scores) == size(dummy_data, 2)
                end
                @testset "only :Lout" begin
                    @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:Lout, 10)))
                end
            end
        end

        qs_types = [RandomOutlierQs, HighConfidenceQs, DecisionBoundaryQs]
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, params), qs_types)
        for qs in qs_objs
            @testset "ModelBasedQs $(typeof(qs))" begin
                scores = qs_score(qs, dummy_data, labelmap(fill(:Lin, 10)))
                @test length(scores) == size(dummy_data, 2)
            end
        end

        qs_types = [NeighborhoodBasedQs, BoundaryNeighborCombinationQs]
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, params), qs_types)
        for qs in qs_objs
            @testset "HybridQs $(typeof(qs))" begin
                scores = qs_score(qs, dummy_data, labelmap(fill(:U, 10)))
                @test length(scores) == size(dummy_data, 2)
            end
        end

        @testset "multiple classifiers" begin

            pools = fill(:U, size(dummy_data, 2))
            init_strategy = SVDD.FixedParameterInitialization(MLKernels.GaussianKernel(0.5), 0.5)

            models = []
            model = SVDD.SVDDneg(dummy_data, pools)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            model = SVDD.VanillaSVDD(dummy_data)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            model = SVDD.SSAD(dummy_data, pools)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            model = SVDD.RandomOCClassifier(dummy_data)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            qs_types = [MinimumMarginQs, MinimumLossQs]
            for qst in qs_types
                for m in models
                    @testset "initialize_qs $qst, $m" begin
                        # check if qs is instantiated
                        qs = OneClassActiveLearning.initialize_qs(qst, model, dummy_data, Dict(:p_inlier => 0.05))
                end
                end
            end

            qs_types = [ExpectedMinimumMarginQs, ExpectedMaximumEntropyQs] ∪ subtypes(ModelBasedQs) ∪ [NeighborhoodBasedQs, BoundaryNeighborCombinationQs]
            for qst in qs_types
                for m in models
                    @testset "initialize_qs $qst, $m" begin
                        # check if qs is instantiated
                        qs = OneClassActiveLearning.initialize_qs(qst, model, dummy_data, params)
                end
                end
            end
        end
    end
end


text:
module ParticleMDI
import Clustering.hclust
import Distributions
import NonUniformRandomVariateGeneration
import Plots
import Plots.PlotMeasures
import StatsBase

# Core code
include("pmdi.jl")
include("__pmdi.jl")
include("misc.jl")
include("update_hypers.jl")

# Datatype specific
include("datatypes/gaussian_cluster.jl")
include("datatypes/categorical_cluster.jl")
include("datatypes/binom_cluster.jl")
include("datatypes/negbinom_cluster.jl")

# Output analysis
include("output_analysis/acf_plots.jl")
include("output_analysis/phi_plots.jl")
include("output_analysis/nclust_plots.jl")
include("output_analysis/consensus_map.jl")
include("output_analysis/feature_select_plots.jl")

export pmdi,
        gaussian_normalise!, coerce_categorical,
        plot_phi_chain, plot_phi_matrix,
        plot_nclust_chain, plot_nclust_hist,
        plot_pmdi_data, get_feature_select_probs,
        generate_psm, consensus_map

end # module


text:
module Lu_elim_I
#
	using Lu
#
		ArgName = "elim_I"
		Descr = "Elimination without pivoting."
		Source = "[Trefethen] Numerical Linear Algebra (p151)."
		Detail = "O(n³)"
#
	type Dat{T}
		prob::Lu.Problem{T}
		U::Lu.MatrixTypesT(T)
		L::Matrix{T}
	#
		Dat() = new()
	end
#
	function construct_dat(T::DataType)
		return Dat{T}()
	end
	function fill_dat{T}(prob::Lu.Problem, dat::Dat{T})
		dat.prob = prob
		dat.U = deepcopy(prob.A); dat.L = eye(T, prob.r);
	end
	function solve_dat{T}(dat::Dat{T}, sol::Lu.Solution{T})
		U = dat.U; L = dat.L; m = dat.prob.c;
		for k in 1:m-1
			for j in k+1:m
				L[j,k] = U[j,k]/U[k,k]
				U[j,k:m] = U[j,k:m] - L[j,k] * U[k,k:m]
			end
		end
		sol.L = tril(L); sol.U = triu(U); sol.solved = all(isfinite(L)) && all(isfinite(U));
	end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
# Functions determining local data ranges in the different pencil
# configurations.

function local_data_range(p, P, N)
    @assert 1 <= p <= P
    a = (N * (p - 1)) ÷ P + 1
    b = (N * p) ÷ P
    a:b
end

# "Complete" dimensions not specified in `dims` with ones.
# Example: if N = 5, dims = (2, 3) and vals = (42, 12), this returns
# (1, 42, 12, 1, 1).
function complete_dims(::Val{N}, dims::Dims{M}, vals::Dims{M}) where {N, M}
    @assert N >= M
    i = 0
    vals_all = ntuple(Val(N)) do n
        if findfirst(dims .== n) === nothing
            1  # this dimension is not included in `dims`, so we put a 1
        else
            i += 1
            vals[i]
        end
    end
    @assert i == M
    vals_all :: Dims{N}
end

# Get axes (array regions) owned by all processes in a given pencil
# configuration.
function get_axes_matrix(decomp_dims::Dims{M},
                         proc_dims::Dims{M},
                         size_global::Dims{N}) where {N, M}
    axes = Array{ArrayRegion{N}, M}(undef, proc_dims)

    # Number of processes in every direction, including those where
    # decomposition is not applied.
    procs = complete_dims(Val(N), decomp_dims, proc_dims)

    for I in CartesianIndices(proc_dims)
        coords = complete_dims(Val(N), decomp_dims, Tuple(I))
        axes[I] = local_data_range.(coords, procs, size_global)
    end

    axes
end


text:
# The following routines multiply
#
#  C[:,j] = α*A*B[:,j] + β* C[:,j]
#
#  But are divided into cases to minimize the calculated
#  rows.  We use A is n x ν, B is ν x m

# Equivalent to
#
#  C[1:min(Cl+j,n),j] = α*A[1:min(Cl+j,n),1:min(Bl+j,ν)]*B[1:min(Bl+j,ν),j]
#                       + β*C[1:min(Cl+j,n),j]
#
@inline function A11_Btop_Ctop_gbmv!(α,β,
                               n,ν,m,j,
                               sz,
                               a,Al,Au,sta,
                               b,Bl,Bu,stb,
                               c,Cl,Cu,stc)
   # A=BandedMatrix(pointer_to_array(a,(Al+Au+1,ν)),n,Al,Au)
   # B=BandedMatrix(pointer_to_array(b,(Bl+Bu+1,m)),ν,Bl,Bu)
   # C=BandedMatrix(pointer_to_array(c,(Cl+Cu+1,m)),n,Cl,Cu)
   #
   # nr=1:min(Cl+j,n)
   # νr=1:min(Bl+j,ν)
   #
   # cj = α*A[nr,νr]*B[νr,j] + β*C[nr,j]
   #
   # for k in nr
   #     C[k,j]=cj[k-first(nr)+1]
   # end
   #
   # c

   gbmv!('N',min(Cl+j,n), min(Bl+j,ν),
          Al, Au,
          α,
          a, sta,
          b+sz*((j-1)*stb+Bu-j+1), 1, β,
          c+sz*((j-1)*stc+Cu-j+1), 1)
end

# Equivalent to
#
#  C[1:min(Cl+j,n),j] = α*A[1:min(Cl+j,n),p:min(p+Bl+Bu+1,ν)]*
#                               B[p:min(p+Bl+Bu+1,ν),j]
#                       + β*C[1:min(Cl+j,n),j]
# for p = j-B.u
#

@inline function Atop_Bmid_Ctop_gbmv!(α,β,
                               n,ν,m,j,
                               sz,
                               a,Al,Au,sta,
                               b,Bl,Bu,stb,
                               c,Cl,Cu,stc)
   # p=j-Bu

   # A=BandedMatrix(pointer_to_array(a,(Al+Au+1,ν)),n,Al,Au)
   # B=BandedMatrix(pointer_to_array(b,(Bl+Bu+1,m)),ν,Bl,Bu)
   # C=BandedMatrix(pointer_to_array(c,(Cl+Cu+1,m)),n,Cl,Cu)
   #
   # nr=1:min(Cl+j,n)
   # νr=p:min(p+Bl+Bu+1,ν)
   #
   # cj = α*A[nr,νr]*B[νr,j] + β*C[nr,j]
   #
   # for k in nr
   #     C[k,j]=cj[k-first(nr)+1]
   # end
   #
   # c

   gbmv!('N', min(Cl+j,n), min(Bl+Bu+1,ν-j+Bu+1),
           Al+j-Bu-1, Au-j+Bu+1,
           α,
           a+sz*(j-Bu-1)*sta, sta,
           b+sz*(j-1)*stb, 1, β,
           c+sz*((j-1)*stc+Cu-j+1), 1)
end

# Equivalent to
#
#  C[nr,j] = α*A[nr,νr]*B[νr,j] + β*C[nr,j]
#
# for p  = j-B.u
#     nr = p-Au : min(p+Al,n)
#     νr = p    : min(p+Bl+Bu+1,ν)
#
#

@inline function Amid_Bmid_Cmid_gbmv!(α,β,
                               n,ν,m,j,
                               sz,
                               a,Al,Au,sta,
                               b,Bl,Bu,stb,
                               c,Cl,Cu,stc)
   p = j-Bu
   # A = BandedMatrix(pointer_to_array(a,(Al+Au+1,ν)),n,Al,Au)
   # B = BandedMatrix(pointer_to_array(b,(Bl+Bu+1,m)),ν,Bl,Bu)
   # C = BandedMatrix(pointer_to_array(c,(Cl+Cu+1,m)),n,Cl,Cu)
   #
   # nr= j-Cu : min(j+Cl,n)
   # νr= p    : min(p+Bl+Bu+1,ν)
   # if !isempty(nr) && !isempty(νr)
   #     cj = α*A[nr,νr]*B[νr,j] + β*C[nr,j]
   #
   #     for k in nr
   #         C[k,j]=cj[k-first(nr)+1]
   #     end
   # end
   #
   # c
   gbmv!('N', min(Cl+Cu+1,n-j+Cu+1), min(Bl+Bu+1,ν-p+1),
           Al+Au, 0,
           α,
           a+sz*(j-Bu-1)*sta, sta,
           b+sz*(j-1)*stb, 1, β,
           c+sz*(j-1)*stc, 1)
end

# Equivalent to
#
#  C[nr,j] =  β*C[nr,j]
#
# for nr= max(1,j-Cu) : min(j+Cl,n)
#

@inline function Anon_Bnon_C_gbmv!(α,β,
                               n,ν,m,j,
                               sz,
                               a,Al,Au,sta,
                               b,Bl,Bu,stb,
                               c,Cl,Cu,stc)
   # C = BandedMatrix(pointer_to_array(c,(Cl+Cu+1,m)),n,Cl,Cu)
   #
   # nr= max(1,j-Cu) : min(j+Cl,n)
   #
   # if !isempty(nr)
   #     cj =  β*C[nr,j]
   #
   #     for k in nr
   #         C[k,j]=cj[k-first(nr)+1]
   #     end
   # end
   #
   # c

   BLAS.scal!(Cu+Cl+1,β,c+sz*(j-1)*stc,1)
end

# function Amid_Bbot_Cmid_gbmv!(α,β,
#                                n,ν,m,j,
#                                sz,
#                                a,Al,Au,sta,
#                                b,Bl,Bu,stb,
#                                c,Cl,Cu,stc)
#    gbmv!('N', Cl+Cu+1,
#            Al+Au, 0,
#            α,
#            a+sz*(j-Bu-1)*sta, ν-j+1, sta,
#            b+sz*(j-1)*stb, β,
#            c+sz*(j-1)*stc)
# end
#
#
# function Abot_Bbot_Cbot_gbmv!(α,β,
#                                n,ν,m,j,
#                                sz,
#                                a,Al,Au,sta,
#                                b,Bl,Bu,stb,
#                                c,Cl,Cu,stc)
#    gbmv!('N', n-j+C.u+1,
#            A.l+A.u, 0,
#            α,
#            a+sz*(j-B.u-1)*sta, B.l+B.u+1-(j-ν+B.l), sta,
#            b+sz*(j-1)*stb, β,
#            c+sz*(j-1)*stc)
# end

_fill_lmul!(β, A::AbstractArray{T}) where T = β == zero(T) ? fill!(A, β) : lmul!(β, A)

function _num_zeroband_u(A)
    Al, Au = bandwidths(A)
    for b = 0:Al+Au
        any(!iszero, view(A,band(Au-b))) && return b
    end
    return Al+Au+1
end

function _num_zeroband_l(A)
    Al, Au = bandwidths(A)
    for b = 0:Al+Au
        any(!iszero, view(A,band(-Al+b))) && return b
    end
    return Al+Au+1
end

function gbmm!(tA::Char, tB::Char, α::T, A::AbstractMatrix{T}, B::AbstractMatrix{T}, β::T, C::AbstractMatrix{T}) where {T<:BlasFloat}
    if tA ≠ 'N' || tB ≠ 'N'
        error("Only 'N' flag is supported.")
    end
    n,ν = size(A)
    m = size(B,2)

    @assert n == size(C,1)
    @assert ν == size(B,1)
    @assert m == size(C,2)

    Al, Au = bandwidths(A)
    Bl, Bu = bandwidths(B)
    C̃l, C̃u = bandwidths(C)
    Cl,Cu = Al+Bl,Au+Bu

    # prune zero bands
    if C̃u < Cu
        Au_r, Bu_r = _num_zeroband_u(A), _num_zeroband_u(B)
        C̃u ≥  Cu - Au_r - Bu_r || throw(BandError(C, Cu - Au_r - Bu_r))
        A_data = bandeddata(A)
        B_data = bandeddata(B)

        if Au-Au_r < -Al || Bu - Bu_r < -Bl
            _fill_lmul!(β, C)
            return C
        end

        Ã = _BandedMatrix(@views(A_data[Au_r+1:end,:]), n, Al, Au-Au_r)
        B̃ = _BandedMatrix(@views(B_data[Bu_r+1:end,:]), ν, Bl, Bu-Bu_r)
        return materialize!(MulAdd(α, Ã, B̃, β, C))
    end

    if C̃l < Cl
        Al_r, Bl_r = _num_zeroband_l(A), _num_zeroband_l(B)
        C̃l ≥  Cl - Al_r - Bl_r || throw(BandError(C, Cl - Al_r - Bl_r))
        A_data = bandeddata(A)
        B_data = bandeddata(B)

        if Al-Al_r < -Au || Bl - Bl_r < -Bu
            _fill_lmul!(β, C)
            return C
        end

        Ã = _BandedMatrix(@views(A_data[1:end-Al_r,:]), n, Al-Al_r, Au)
        B̃ = _BandedMatrix(@views(B_data[1:end-Bl_r,:]), ν, Bl-Bl_r, Bu)
        return materialize!(MulAdd(α, Ã, B̃, β, C))
    end

    A_data = bandeddata(A)
    B_data = bandeddata(B)
    C̃_data = bandeddata(C)

    # scale extra bands
    _fill_lmul!(β, view(C̃_data, 1:C̃u-Cu,:))
    _fill_lmul!(β, view(C̃_data, (C̃u+Cl+1)+1:size(C̃_data,1),:))
    C_data = view(C̃_data, (C̃u-Cu+1):(C̃u+Cl+1), :) # shift to bands we will write to
    _gbmm!(α, A_data, B_data, β, C_data, (n,ν,m), (Al, Au), (Bl, Bu), (Cl, Cu))

    C
end

function _gbmm!(α::T, A_data, B_data, β, C_data, (n,ν,m), (Al, Au), (Bl, Bu), (Cl, Cu)) where T
    a = pointer(A_data)
    b = pointer(B_data)
    c = pointer(C_data)
    sta = stride(A_data,2)
    stb = stride(B_data,2)
    stc = stride(C_data,2)
    sz = sizeof(T)

    # Multiply columns j where B[1,j]≠0: A is at 1,1 and C[1,j]≠0
    for j = 1:min(m,1+Bu)
        A11_Btop_Ctop_gbmv!(α,β,
                            n,ν,m,j,
                            sz,
                            a,Al,Au,sta,
                            b,Bl,Bu,stb,
                            c,Cl,Cu,stc)
    end

    # Multiply columns j where B[k,j]=0 for k<p=(j-B.u-1), A is at 1,1+p and C[1,j]≠0
    # j ≤ ν + B.u since then 1+p ≤ ν, so inside the columns of A

    for j = 2+Bu:min(1+Cu,ν+Bu,m)
        Atop_Bmid_Ctop_gbmv!(α,β,
                            n,ν,m,j,
                            sz,
                            a,Al,Au,sta,
                            b,Bl,Bu,stb,
                            c,Cl,Cu,stc)
    end

    # multiply columns where A and B are mid and C is bottom
    for j = 2+Cu:min(m,ν+Bu,n+Cu)
        Amid_Bmid_Cmid_gbmv!(α,β,
                            n,ν,m,j,
                            sz,
                            a,Al,Au,sta,
                            b,Bl,Bu,stb,
                            c,Cl,Cu,stc)
    end

    # scale columns of C by β that aren't impacted by α*A*B
    _fill_lmul!(β, view(C_data, :, ν+Bu+1:min(m,n+Cu)))
end


text:
# =============== String Conversion ===============

# Convert vector of strings to vector of integers to enable histogram for string vector
function stringvec2intvec{Ts<:AbstractString}(s::Array{Ts,1})
    u = unique(s)
    x = zeros(Int, length(s))
    for i = 1:length(u)
        x[find(s .== u[i])] = i
    end
    return x
end

# =============== Bin Calcs ===============

# Calculate number of bins in histogram for vector of floating point numbers

# Freedman-Diaconis rule
function nBins_fdr{Tf<:AbstractFloat}(x::Array{Tf,1})
    n = length(x)
    iqr = quantile(x, 0.75) - quantile(x, 0.25)
    maxmin = maximum(x) - minimum(x)
    fdr = 2 * iqr * n^(-1/3)
    nbins = ceil(Int, maxmin/fdr)
end

# Scott's Rule
function nBins_sr{Tf<:AbstractFloat}(x::Array{Tf,1})
    n = length(x)
    iqr = quantile(x, 0.75) - quantile(x, 0.25)
    maxmin = maximum(x) - minimum(x)
    sr = 3.5*std(x)*(length(x))^(-1/3)
    nbins = ceil(Int, maxmin/sr)
end

# Square Root method
function nBins_sqrt{Tf<:AbstractFloat}(x::Array{Tf,1})
    n = length(x)
    nbins = ceil(Int, sqrt(n))
end

# =============== Entropy ===============

# Entropy fuctions for vectors of integers, real numbers and strings

function entropy{Ti<:Integer}(x::Array{Ti,1})
    n = length(x)                            # n = total number of values
    nBins = length(unique(x))                # nBins = number of unique values
    c = hist(x, nBins)[2]                    # c = vector showing count of values of each outcome (i.e. in each bucket)
    p = c ./ n                               # p = probability of each outcome (i.e. bucket)
    #p = P[P .!= 0]                          # remove zeros [not necessary]
    H = -sum(p .* log2(p))                   # H = entropy of x
end

function entropy{Tf<:AbstractFloat}(x::Array{Tf,1})
    n = length(x)                            # n = total number of values
    nBins = nBins_fdr(x)                     # nBins = number of bins using Freedman-Diaconis rule
    c = hist(x, nBins)[2]                    # c = vector showing count of values in each bucket
    P = c ./ n                               # P = probability of each bucket
    p = P[P .!= 0]                           # remove zeros
    H = -sum(p .* log2(p))                   # H = entropy of x
end

function entropy{Ts<:AbstractString}(s::Array{Ts,1})
    x = stringvec2intvec(s)                  # Change to vector of integers
    entropy(x)                               # Calculate entropy using Integer method
end

# =============== Joint Entropy ===============

# Joint Entropy functions for vectors of integers, real numbers and strings

function jointEntropy{Ti<:Integer}(x::Array{Ti,1}, y::Array{Ti,1})
    n = length(x)                             # n = total number of values (same length for both x and y)
    nBinsx = length(unique(x))                # nBinsx = number of unique values of vector x
    nBinsy = length(unique(y))                # nBinsy = number of unique values of vector y
    C = hist2d([x y], nBinsx, nBinsy)[3]      # C = matrix showing count of values in each joint outcome (i.e. in each bucket)
    P = C ./ n                                # P = probability of each joint outcome (i.e. bucket)
    p = P[P .!= 0]                            # remove zeros
    H = -sum(p .* log2(p))                    # H = joint entropy of x and y
end

function jointEntropy{Tf<:AbstractFloat, Ti<:Integer}(x::Array{Tf,1}, y::Array{Ti,1})
    n = length(x)                             # n = total number of values (same length for both x and y)
    nBinsx = nBins_fdr(x)                     # nBinsx = number of bins of vector x using Freedman-Diaconis rule
    nBinsy = length(unique(y))                # nBinsy = number of unique values of vector y
    C = hist2d([x y], nBinsx, nBinsy)[3]      # C = matrix showing count of values in each bucket
    P = C ./ n                                # P = probability of each bucket
    p = P[P .!= 0]                            # remove zeros
    H = -sum(p .* log2(p))                    # H = joint entropy of x and y
end

function jointEntropy{Ti<:Integer, Tf<:AbstractFloat}(x::Array{Ti,1}, y::Array{Tf,1})
    n = length(x)                             # n = total number of values (same length for both x and y)
    nBinsx = length(unique(x))                # nBinsx = number of unique values of vector x
    nBinsy = nBins_fdr(y)                     # nBinsy = number of bins of vector y using Freedman-Diaconis rule
    C = hist2d([x y], nBinsx, nBinsy)[3]      # C = matrix showing count of values in each bucket
    P = C ./ n                                # P = probability of each bucket
    p = P[P .!= 0]                            # remove zeros
    H = -sum(p .* log2(p))                    # H = joint entropy of x and y
end

function jointEntropy{Tf<:AbstractFloat}(x::Array{Tf,1}, y::Array{Tf,1})
    n = length(x)                             # n = total number of values (same length for both x and y)
    nBinsx = nBins_fdr(x)                     # nBinsx = number of bins of vector x using Freedman-Diaconis rule
    nBinsy = nBins_fdr(y)                     # nBinsy = number of bins of vector y using Freedman-Diaconis rule
    C = hist2d([x y], nBinsx, nBinsy)[3]      # C = matrix showing count of values in each bucket
    P = C ./ n                                # P = probability of each bucket
    p = P[P .!= 0]                            # remove zeros
    H = -sum(p .* log2(p))                    # H = joint entropy of x and y
end

function jointEntropy{Ts<:AbstractString}(s1::Array{Ts,1}, s2::Array{Ts,1})
    x = stringvec2intvec(s1)
    y = stringvec2intvec(s2)
    jointEntropy(x, y)
end

function jointEntropy{Ts<:AbstractString, Tr<:Real}(s::Array{Ts,1}, y::Array{Tr,1})
    x = stringvec2intvec(s)
    jointEntropy(x, y)
end

function jointEntropy{Tr<:Real, Ts<:AbstractString}(x::Array{Tr,1}, s::Array{Ts,1})
    y = stringvec2intvec(s)
    jointEntropy(x, y)
end

# =============== Mutual Information ===============

# Mutual Information function for vectors of integers, real numbers and strings
function mutualInformation(x, y)
    entropy(x) + entropy(y) - jointEntropy(x, y)
end

# Mutual Information function for vectors within dataframes
function mutualInformation(MIdf::DataFrame, dx::Symbol, dy::Symbol)
    df = MIdf[:, [dx, dy]]
    complete_cases!(df)
    x = convert(Array, df[dx])
    y = convert(Array, df[dy])
    mutualInformation(x, y)
    #entropy(x) + entropy(y) - jointEntropy(x, y)
end

# Mutual information table from dataframe
function mutualInformation(MIdf::DataFrame, target::Symbol)
    variables = names(MIdf)
    targetindex = findfirst(variables, target)
    splice!(variables, targetindex)
    MIvec = Real[]
    for var in variables
        df = MIdf[:, [var, target]]
        complete_cases!(df)
        x = convert(Array, df[var])
        y = convert(Array, df[target])
        push!(MIvec, mutualInformation(x, y))
    end
    df = DataFrame(Name=variables)
    insert!(df, 2, MIvec, :MI)
    sort(df, cols=:MI, rev=true)
end


text:
FIXTURES = joinpath(@__DIR__, "fixtures")
SBML_FILE = joinpath(FIXTURES, "_model_v1.xml")

# test getmodel
@test_nowarn getmodel(SBML_FILE)
model = getmodel(SBML_FILE)

# test getparameters
parameters = getparameters(model)
true_parameters = Pair{ModelingToolkit.Num,Float64}[a0 => 1.0,
                                                   b0 => 0.0,
                                                   k1 => 0.8,
                                                   k2 => 0.6,
                                                   compartment => 1.0]
@test repr(parameters) == repr(true_parameters)

# test getinitialconditions
@parameters t
u0 = getinitialconditions(model)
true_u0 = Pair{ModelingToolkit.Num,Float64}[A => 1.0,
                                            B => 0.0]
@test repr(u0) == repr(true_u0)

# test getodes
rxs = getreactions(model)
truerxs = Reaction[Reaction(compartment*k1*A, [A], [B], [1.0], [1.0], only_use_rate=true), Reaction(compartment*k2*B, [B], [A], [1.0], [1.0], only_use_rate=true)]
@test repr(rxs) == repr(truerxs)
rs  = ModelingToolkit.ReactionSystem(rxs, t, [item.first for item in true_u0], [item.first for item in true_parameters])
odesys = convert(ModelingToolkit.ODESystem, rs)
@test repr(ModelingToolkit.get_iv(odesys)) == "t"
@test repr(ModelingToolkit.get_states(odesys)) == "Term{Real}[A(t), B(t)]"

# test sbml2odesystem
sys = sbml2odesystem(SBML_FILE)
true_eqs = ModelingToolkit.Equation[
            _Differential(A) ~ -1.0 * (compartment * k1 * A) + 1.0 * (compartment * k2 * B),
            _Differential(B) ~ 1.0 * (compartment * k1 * A) - 1.0 * (compartment * k2 * B)]
true_sys = ModelingToolkit.ODESystem(true_eqs)
@test ModelingToolkit.get_iv(sys) == ModelingToolkit.get_iv(true_sys)
@test isequal(ModelingToolkit.get_states(sys), ModelingToolkit.get_states(true_sys))
@test ModelingToolkit.get_default_u0(sys) == Dict(true_u0)
@test collect(values(ModelingToolkit.get_default_p(sys))) == collect(values(Dict(true_parameters)))

# test sbml2odeproblem
prob = sbml2odeproblem(SBML_FILE,(0.0,10.0))
@test prob.u0 == [1.0, 0.0]  # "Pair{ModelingToolkit.Num,Float64}[A => 1.0,B => 0.0]"
@test prob.tspan == (0.0, 10.0)
@test prob.p == [1.0, 0.8, 0.6]  # Todo: try removing the Set() when new ODEproblem version is released.
@test_nowarn OrdinaryDiffEq.solve(prob,OrdinaryDiffEq.Tsit5())
@test_nowarn sbml2odeproblem(SBML_FILE,(0.0,1.0),jac=false)

# test simulatesbml
sol1 = simulatesbml(SBML_FILE,(0.0,10.0))
@test sol1[end] ≈ [0.42857284277223784, 0.5714271572277619]
@test_nowarn simulatesbml(SBML_FILE,(0.0,1.0),saveat=0.5,jac=false,solver=OrdinaryDiffEq.Vern9())
@test_nowarn simulatesbml(joinpath(FIXTURES, "model_Vinod_FEBS2015.xml"),(0.0,1.0))


text:
# tests of basic utilities for working with HTML

import Gumbo: HTMLNode, NullNode

# convenience constructor works
@test HTMLElement(:body) == HTMLElement{:body}(HTMLNode[],
                                               NullNode(),
                                               Dict{AbstractString,AbstractString}())

# accessing tags works
@test HTMLElement(:body) |> tag == :body

let
    elem = HTMLElement{:body}(HTMLNode[], NullNode(), Dict("foo" => "bar"))
    @test getattr(elem, "foo") == "bar"
    @test getattr(elem, "foo", "baz") == "bar"
    @test getattr(elem, "bar", "qux") == "qux"
    @test getattr(() -> "qux", elem, "bar") == "qux"
end


text:
# n = mult(3, 4) #> ErrorException("mult not defined"))
function mult(x, y)
    println("x is $x and y is $y")
    x * y
end

n = mult(3, 4) #> 12

function mult(x, y)
    println("x is $x and y is $y")
    if x == 1
        return y
    end
    x * y
end

function multi(n, m)
    n * m, div(n, m), n % m
end
x, y, z = multi(8, 2)  #> (16,4,0)

# variable number of arguments:
function varargs(n, m, args...)
    println("arguments: $n $m $args")
end
varargs(1, 2, 3, 4) #> arguments: 1 2 (3, 4)

function varargs2(args...)
    println("arguments2: $args")
end

x = (3, 4)
varargs2(1, 2, x...) # args is (1,2,3,4)
x = [10, 11, 12]
varargs2(1, 2, x...) # args is (1,2,10,11,12)

function insert_elem(arr)
    push!(arr, -10)
end

arr = [ 2, 3, 4]
insert_elem(arr)
println("arr is now $arr") #> arr is now [ 2, 3, 4, -10 ]

# typing the arguments:
function mult(x::Float64, y::Float64)
    println("x is $x and y is $y")
    x * y
end
# following gives error:
# mult(5, 6)
# ERROR: `mult` has no method matching mult(::Int64, ::Int64)

# one-line assignment syntax:
mult(x, y) = x * y
f(x, y) = x^3 - y + x*y
println(f(3,2) ) # 31

∑(x,y) = x + y
∑(3, 4) #> 7


text:
#-----------------------------------------------------------------------------------------------------
# Functions for Operations.

# Operation is the collection of "Operator"s, together with a "Basis".
# Functions on "Operation"s are design to mimic the behavior of matrices.
#-----------------------------------------------------------------------------------------------------
"""
    Operator{T1<:Number, T2<:Integer}

Type that contains 2 fields:

1. `mat`: Matrices representation of local operator.
2. `inds`: Indeces of sites the operator acts on.
"""
struct Operator{T1<:Number, T2<:Integer}
    mat::Matrix{T1}
    inds::Vector{T2}
end
#-----------------------------------------------------------------------------------------------------
# Basic functions on "Operator":

# 1. *: Multiplication by a number.
# 2. /: Division by a number.
#-----------------------------------------------------------------------------------------------------
*(c::Number, o::Operator) = Operator(c * o.mat, o.inds)
/(o::Operator, c::Number) = Operator(o.mat / c, o.inds)
#-----------------------------------------------------------------------------------------------------
# Operator to fill vector/matrix.

# To finally get the multiplication function, here we focus on a single row manipulation:

# 1. Start with a given product state(digits), we cut a segment from it;
# 2. Iterate all possible digits in the segment, and read the matrix element of operator;
# 3. Add to the vector/matrix element.
#-----------------------------------------------------------------------------------------------------
"""
    addtovec!(vec::AbstractVector, opt::Operator, basis::Basis, coeff::Number=1)

Elementary multiplication step. Calculate single vector element of an operator * vector.

- `vec`: Vector to fill
- `opt`: Operator
- `basis`: Indicate the specific row
- `coeff`: The vector element of input state
"""
function addtovec!(vec::AbstractVector, opt::Operator, basis::Basis, coeff::Number=1)
    basis_view = view(basis, opt.inds)              # Create a view on the segment of the digits
    index_view = index(basis_view)                  # Get the initial index of the viewed segment
    opt_column = opt.mat[:, index_view] * coeff     # Get the elements in row of index_view
    row_number = length(opt_column)
    for k = 1:row_number
        change!(basis_view, k)
        i = index(basis)
        vec[i] += opt_column[k]
    end                                             # Fill the vector
    change!(basis_view, index_view)                 # Reset the segment
end
#-----------------------------------------------------------------------------------------------------
"""
    addtovecs!(vecs::AbstractMatrix, opt::Operator, basis::Basis, coeff::AbstractVector)

Elementary multiplication step. Calculate single row of matrix element of an operator * matrix.

- vecs : Vectors to fill
- opt  : Operator
- basis: Indicate the specific row
- coeff: The row of matrix elements of input states
"""
function addtovecs!(vecs::AbstractMatrix, opt::Operator, basis::Basis, coeff::AbstractVector{<:Number})
    basis_view = view(basis, opt.inds)              # Create a view on the segment of the digits
    index_view = index(basis_view)                  # Get the initial index of the viewed segment
    opt_column = opt.mat[:, index_view]             # Get the elements in row of index_view
    row_number = length(opt_column)                 # Fill the matrix
    for k = 1:row_number
        change!(basis_view, k)
        i = index(basis)
        vecs[i, :] += opt_column[k] * coeff
    end
    change!(basis_view, index_view)                 # Reset the segment
end
#-----------------------------------------------------------------------------------------------------
# Type: Operation
#-----------------------------------------------------------------------------------------------------
export Operation
"""
    Operation{OptType<:AbstractVector{<:Operator}, BasType<:Basis}

Type that contains 2 fields:

1. `opts`: List of operators
2. `basis`: Basis for the system
"""
struct Operation{T1 <: Number, T2 <: Integer, T3}
    opts::Vector{Operator{T1, T2}}
    basis::Basis{T3}
end
eltype(::Operation{T1, T2, T3}) where T1 where T2 where T3 = T1
function size(opt::Operation)
    basis = opt.basis
    dim = basis.base ^ basis.len
    (dim, dim)
end
#-----------------------------------------------------------------------------------------------------
# Basis Operation initiation
#-----------------------------------------------------------------------------------------------------
export operation
"""
    operation(mats, inds, len=0; base=0)

Canonical construction method for `Operation` object.

- `mats`: List of matrix representation of operators.
- `inds`: List of sites the operastors act on.
- `len`: Specify the size of the system. Default `len=0` will induce auto deduction.
- `base`: Quantum number of each sites. Default `base=0` will induce auto deduction.
"""
function operation(
    mats::AbstractVector{<:AbstractMatrix},
    inds::AbstractVector{<:AbstractVector},
    len::Integer=0;
    base::Integer=0
)
    B = begin
        b = base==0 ? round(Int64, size(mats[1], 1)^(1/length(inds[1]))) : Int64(base)
        l = len==0 ? maximum(maximum.(inds)) : Int64(len)
        basis(b, l)
    end
    mat_type = promote_type(eltype.(mats)...)
    ind_type = promote_type(eltype.(inds)...)
    O = [Operator(Array{mat_type}(mats[i]), Array{ind_type}(inds[i])) for i=1:length(mats)]
    Operation(O, B)
end
#-----------------------------------------------------------------------------------------------------
# Basic Functions for type Operation
#-----------------------------------------------------------------------------------------------------
*(c::Number, o::Operation) = Operation(c .* o.opts, o.basis)
/(o::Operation, c::Number) = Operation(o.opts ./ c, o.basis)
+(opt1::Operation, opt2::Operation) = Operation(vcat(opt1.opts, opt2.opts), opt1.basis)
-(o1::Operation, o2::Operation) = o1 + ((-1) * o2)
function sum(ol::AbstractVector{<:Operation})
    basis = ol[1].basis
    opts = vcat([oi.opts for oi in ol]...)
    Operation(opts, basis)
end
#-----------------------------------------------------------------------------------------------------
# Operation to fill vector/matrix:

# The method is basically iterate each operator to the given digits-represented basis
#-----------------------------------------------------------------------------------------------------
"""
    addtovec!(vec::AbstractVector, opt::Operation, coeff::Number=1)

Elementary multiplication step. Calculate single vector element of an operation * vector.
The row information is stored in the basis of operation.

- `vec`: Vector to fill.
- `opt`: Operation.
- `coeff`: The vector element of input state.
"""
function addtovec!(vec::AbstractVector, opt::Operation, coeff::Number=1)
    basis = opt.basis
    opts = opt.opts
    num_of_opts = length(opts)
    for i = 1:num_of_opts
        addtovec!(vec, opts[i], basis, coeff)
    end # Multiply by each operator and add them all to vector
end
#-----------------------------------------------------------------------------------------------------
"""
    addtovecs!(vecs::AbstractMatrix, opt::Operation, coeff::AbstractVector)

Elementary multiplication step. Calculate single row of matrix element of an operation * matrix.
The row information is stored in the basis of operation.

- `vecs`: Vectors to fill
- `opt`: Operation
- `coeff`: The row of matrix elements of input states
"""
function addtovecs!(vecs::AbstractMatrix, opt::Operation, coeff::AbstractVector{<:Number})
    basis = opt.basis
    opts = opt.opts
    num_of_opts = length(opts)
    for i = 1:num_of_opts
        addtovecs!(vecs, opts[i], basis, coeff)
    end # Multiply by each operator and add them all to vectors
end
#-----------------------------------------------------------------------------------------------------
# Multiplication

# The idea is to iterate all product state basis, and get all the vector/matrix elements
#-----------------------------------------------------------------------------------------------------
"""
    mul!(vec::AbstractVector, opt::Operation, state::AbstractVector)

Full multiplication for operation and vector.

- `vec`: Vector to fill.
- `opt`: Operation.
- `state`: Input vector.
"""
function mul!(vec::AbstractVector, opt::Operation, state::AbstractVector)
    basis = opt.basis
    for j = 1:length(state)
        change!(basis, j)
        addtovec!(vec, opt, state[j])
    end
end
#-----------------------------------------------------------------------------------------------------
"""
    mul!(mat::AbstractVector, opt::Operation, state::AbstractVector)

Full multiplication for operation and matrix.

- `mat`: Matrix to fill.
- `opt`: Operation.
- `states`: Input states(matrix).
"""
function mul!(mat::AbstractMatrix, opt::Operation, states::AbstractMatrix)
    basis = opt.basis
    for j = 1:size(states, 1)
        change!(basis, j)
        addtovecs!(mat, opt, states[j, :])
    end
end
#-----------------------------------------------------------------------------------------------------
"""
    *(opt::Operation, vec_or_mat::AbstractVecOrMat)

General multiplication for operation and vector/matrix.

- `opt`: Operation.
- `vec_or_mat`: Input vector/matrix.
"""
function *(opt::Operation, vec_or_mat::AbstractVecOrMat)
    ctype = promote_type(eltype(opt), eltype(vec_or_mat))
    out = zeros(ctype, size(vec_or_mat))
    mul!(out, opt, vec_or_mat)
    out
end
#-----------------------------------------------------------------------------------------------------
# Fill matrix
#-----------------------------------------------------------------------------------------------------
export fillmat!
"""
    fillmat!(mat::AbstractMatrix, opt::Operation)

Fill the zero matrix with matrix element from operation.

- `mat`: Matrix to fill.
- `opt`: Operation.
"""
function fillmat!(mat::AbstractMatrix, opt::Operation)
    basis = opt.basis
    row_number = size(mat, 1)
    for j = 1:row_number
        change!(basis, j)
        addtovec!(view(mat, :, j), opt)
    end
end
#-----------------------------------------------------------------------------------------------------
"""
    Array(opt::Operation)

Return the full matrix form of `opt`.
"""
function Array(opt::Operation)
    mat = zeros(eltype(opt), size(opt))
    fillmat!(mat, opt)
    mat
end


text:
struct FurutaParams{T}
    α::T
    β::T
    γ::T
    δ::T
    Ku::T
    fϕ::T
    fθ::T
    h::T
    max_speed::T
    max_torque::T
    noise::T
    friction::Symbol
end

mutable struct SimulatedFurutaPendulum{T, R<:AbstractRNG} <: SimulatedProcess
    params::FurutaParams{T}
    x::Vector{T} # θ, θdot, ϕ, ϕdot = base angle & speed, arm angle & speed
    u::T
    rng::R
end

"""
    SimulatedFurutaPendulum(;
        J = 1.54e-4,
        M = 0,
        ma = 0,
        mp = 5.44e-3,
        la = 4.3e-2,
        lp = 6.46e-2,
        friction = :none,
        fϕ = 0.0238,
        fθ = 0.000008,
        R = 0.13,
        Ka = 0.03,
        h = 0.002,
        x0 = zeros(4),
        max_speed = 100,
        max_torque = 1,
        noise = 0.01,
        friction = :viscous,
        rng = Random.GLOBAL_RNG,
    )

Creates an instance of a simulator for the furutapendulum.
It is made up of two inertial bodies, a center pillar with
intertia `J` connected to a horizontal arm with length `la`
and homogenuously distributed mass `ma`. The other end of the
arm is connected to the pendulum which has length `lp`, homogenuously
distributed mass `mp` and point distributed mass `M` at the end.
The `ϕ`-joint is in the center pillar while the `θ`-joint is for the
pendulum arm. There is significant `friction` in the `ϕ`-joint which
is modeled using either coulomb and viscous friction `:viscous`,
or coulomb friction and stiction `:stiction` (for more info see Gäfvert
https://portal.research.lu.se/portal/files/4453844/8727127.pdf).
The simulation takes steps of `h` seconds, and the initial value `x0`
is for `[ϕ, ϕdot, θ, θdot]`.
There is `max_speed` and `max_torque` limiting the speed and torque, and
`noise` that adds white noise with this scaling factor on top of the
measurements.
"""
function SimulatedFurutaPendulum(;
        J = 1.54e-4,
        M = 0,
        ma = 0,
        mp = 5.44e-3,
        la = 4.3e-2,
        lp = 6.46e-2,
        fϕ = 0.0238,
        fθ = 0.000008,
        R = 0.13,
        Ka = 0.03,
        h = 0.002, # Simulation time, TODO find what this actually is
        max_speed = 100,
        max_torque = 0.04,
        x0 = [0.0, 0.0, π, 0.0], # Start down as default
        noise = 0.01,
        friction::Symbol = :viscous,
        rng = Random.GLOBAL_RNG,
    )
    g = 9.81
    α = J+(M+ma/3+mp)*la^2
    β = (M+mp/3)*lp^2
    γ = (M+mp/2)*la*lp
    δ = (M+mp/2)*g*lp
    Ku = Ka / R

    params = FurutaParams{Float64}(α, β, γ, δ, Ku, fϕ, fθ, h, max_speed, max_torque, noise, friction)
    SimulatedFurutaPendulum{Float64, typeof(rng)}(params, x0, 0.0, rng)
end

function step!(p::SimulatedFurutaPendulum, dt)
    x = p.x
    u = p.u
    h = p.params.h

    steps = round(Int, dt / h)
    steps ≈ dt / h || throw(ArgumentError("`dt`=$(dt) need to be a multiple of the simulation step `h`=$(h)."))
    for _ in 1:steps
        #RK3/8
        k1 = f(p, x, u)
        k2 = f(p, x + h * k1 / 3, u)
        k3 = f(p, x + h * (-k1 / 3 + k2), u)
        k4 = f(p, x + h * (k1 - k2 + k3), u)
        x .+= h .* (k1 .+ 3 .* k2 .+ 3 .* k3 .+ k4) ./ 8
    end

    x[1] = mod2pi(x[1])
    x[3] = mod2pi(x[3])

    # Should this be here?
    x[2] = clamp(x[2] , -p.params.max_speed, p.params.max_speed)
    x[4] = clamp(x[4] , -p.params.max_speed, p.params.max_speed)

    nothing
end

# See https://portal.research.lu.se/portal/files/4453844/8727127.pdf
# for source of dynamic equations
function f(p::SimulatedFurutaPendulum, x, u)
    # Base angle/speed, arm angle/speed
    ϕ, ϕdot, θ, θdot = x

    @unpack α, β, γ, δ, Ku, fϕ, fθ = p.params

    # clamp(Ku * p.u, -p.params.max_torque, p.params.max_torque)
    τϕm = Ku * u # Torque from motor (approximating that transients between voltage and current are fast)
    τϕf = -fϕ * ϕdot
    τϕ = τϕm + τϕf
    # theta joint
    τθ = -fθ * θdot

    ψ = 1 / (α*β - γ^2 + (β^2 + γ^2)*sin(θ)^2)
    dϕ    = ϕdot
    dϕdot = ψ*(β*γ*(sin(θ)^2-1)*sin(θ)*ϕdot^2 - 2*β^2*cos(θ)*sin(θ)*ϕdot*θdot
        + β*γ*sin(θ)*θdot^2 - γ*δ*cos(θ)*sin(θ) + β*τϕ - γ*cos(θ)*τθ)
    dθ    = θdot
    dθdot = ψ * (β*(α+β*sin(θ)^2)*cos(θ)*sin(θ)*ϕdot^2
        + 2*β*γ*(1-sin(θ)^2)*sin(θ)*ϕdot*θdot - γ^2*cos(θ)*sin(θ)*θdot^2
        + δ*(α+β*sin(θ)^2)*sin(θ) - γ*cos(θ)*τϕ + (α+β*sin(θ)^2)*τθ)

    return [dϕ, dϕdot, dθ, dθdot]
end

# To allow reading values on the same low level interface as for the hardware
# using read(furuta.port)
# Currently only implemented for energy readings
struct ValueContainer
    v::Float64
end

function Base.getproperty(p::SimulatedFurutaPendulum, s::Symbol)
    K0 = 0.0000652 # kinetic energy constant base: 0.5*moment of inertia, kgm^2
    K1 = 0.00000387 # kinetic energy constant pendulum

    if s === :arm_energy
        ϕ, ϕdot, θ, θdot = p.x
        arm_energy = (K1 * θdot^2 + 0.0054 * 9.8 * 0.07/2 * (1 - cos(θ)))
        return ValueContainer(arm_energy)
    elseif s === :total_energy
        ϕ, ϕdot, θ, θdot = p.x
        arm_energy = (K1 * θdot^2 + 0.0054 * 9.8 * 0.07/2 * (1 - cos(θ)))
        base_energy = K0 * ϕdot^2
        return ValueContainer(base_energy + arm_energy)
    else
        return getfield(p, s)
    end
end

Base.read(c::ValueContainer) = c.v


text:
using NeRCA
using ROy
using Sockets

calib = Calibration("scripts/NeRCA_00000043_20022019.detx")
calib = Calibration("/home/tgal/data/detx/NeRCA_-00000001_20171212.detx")

client = NeRCA.CHClient(ip"127.0.0.1", 5553, ["IO_EVT"])

try
    for message in client
        event = read(IOBuffer(message.data), NeRCA.DAQEvent)
        hits = event.snapshot_hits
        chits = NeRCA.calibrate(calib, hits)
        thits = event.triggered_hits
        cthits = NeRCA.calibrate(calib, thits)
        println(length(cthits))
    end
catch InterruptException
    println("ok")
end


text: using Pkg if isfile("Project.toml") && isfile("Manifest.toml") Pkg.activate(".") end


text:
using InteractiveUtils
using PanelCutting
using Unitful
using UnitfulUS
# using UnitfulCurrency
using Test
using Logging
using NahaGraphs

#=
include("test_graph_rewrite.jl")
=#

@testset "AllOf" begin
    ao = AllOf([1, 2, 3], (:a, :b))
    @test length(ao) == 5
    @test ao[2] == 2
    @test ao[5] == :b
    @test collect(ao) == [1, 2, 3, :a, :b]
end

function allsubtypes(t, result=Set())
    push!(result, t)
    for t1 in subtypes(t)
	allsubtypes(t1, result)
    end
    return result
end

include("example_makers.jl")

@testset "cut math" begin
    at = 8u"inch"
    plw = [2u"ft" 1u"ft"]
    pxy = [1u"inch" 2u"inch"]
    vmult(a, b) = map(*, a, b)
    lw(axis) = (at * unit_vector(axis) +
        vmult(plw, unit_vector(other(axis))))
    @test lw(LengthAxis()) == [at plw[2]]
    @test lw(WidthAxis()) == [plw[1] at]
end

@testset "PanelIdentity" begin
    for paneltype in filter(isconcretetype, allsubtypes(AbstractPanel))
        @test example(paneltype) != example(paneltype)
    end
end

@testset "PanelIntegrity" begin
    # Make sure we can't associate the wrong kinds or sizes of panels
    # with one another
end

@testset "cut: stock too short" begin
    axis = LengthAxis()
    panel1 = BoughtPanel(AvailablePanel("30 by 60", 30u"inch", 60u"inch", 20))
    panels = cut(panel1, axis, panel1.length * 1.1)
    @test isempty(panels)
end

# A bunch of common test assertions relating a Panel to the
# CuttablePanel it was cut from.
function panel_test(panel::Panel, kerf)
    parent = panel.cut_from
    if panel.x == parent.x && panel.y == parent.y
        # The first panel returned by cut, the intended panel
        if panel.cut_axis isa LengthAxis
            @test panel.length == panel.cut_at
            @test panel.width == parent.width
        else
            @test panel.width == panel.cut_at
            @test panel.length == parent.length
        end
    else
        # the off-cut
        offset = panel.cut_at + kerf
        if panel.cut_axis isa LengthAxis
            @test panel.length == parent.length - offset
            @test panel.width == parent.width
            @test panel.x == parent.x + offset
            @test panel.y == parent.y
        else
            @test panel.width == parent.width - offset
            @test panel.length == parent.length
            @test panel.x == parent.x
            @test panel.y == parent.y + offset
        end
    end
end

@testset "trivial cut" begin
    axis = LengthAxis()
    kerf = (1/8)u"inch"
    cut_cost = money(1)
    panel1 = BoughtPanel(AvailablePanel("30 by 60", 30u"inch", 60u"inch", 20))
    panels = cut(panel1, axis, panel1.length; kerf=kerf, cost=cut_cost)
    @test length(panels) == 1
    panel2 = panels[1]
    @test panel2 isa Panel
    @test panel2.cut_from == panel1
    @test panel2.cut_at == panel1.length
    @test panel2.cut_axis == axis
    @test panel2.cost == panel1.cost
    @test panel2.x == panel1.x
    @test panel2.y == panel1.y
    panel_test(panel2, kerf)
end

@testset "close shave" begin
    kerf= (1/8)u"inch"
    cut_cost = 1
    panel1 = BoughtPanel(AvailablePanel("30 by 60", 30u"inch", 60u"inch", 20))
    panels = cut(panel1, LengthAxis(),
                 panel1.length - kerf;
                 kerf=kerf, cost=cut_cost)
    @test length(panels) == 1
    panel2 = panels[1]
    @test panel2.cut_from == panel1
    @test panel2 isa Panel
    @test panel2.cost == panel1.cost + cut_cost
    panel_test(panel2, kerf)
end

@testset "Cut" begin
    kerf= (1/8)u"inch"
    cut_cost = 1
    panel1 = BoughtPanel(AvailablePanel("30 by 60", 30u"inch", 60u"inch", 20))
    panel2, panel3 = cut(panel1, LengthAxis(), 25u"inch";
                         kerf=kerf, cost=cut_cost)
    total_cost = panel1.cost + cut_cost
    @test panel2.length == 25u"inch"
    @test panel2.width == 30u"inch"
    @test panel2.cut_from == panel1
    @test panel3.cut_from == panel1
    @test panel2.cut_at == panel3.cut_at
    @test panel2.cut_axis == panel3.cut_axis
    panel_test(panel2, kerf)
    panel_test(panel3, kerf)
    total_area = area(panel2) + area(panel3)
    @test panel1.length == panel2.length + kerf + panel3.length
    @test panel2.cost == total_cost * area(panel2) / total_area
    @test panel3.cost == total_cost * area(panel3) / total_area
end

# Each element of searcher.cheapest.bought should be the progenitor of
# at leat one searcher.cheapest.finished panel.
function unuesd_bought_panels(state::SearchState)
    finished_bought = Set(map(progenitor, state.finished))
    symdiff(Set(state.bought), finished_bought)
end

@testset "Search: fits in one" begin
    supplier = Supplier(
        name="test",
        cost_per_cut=1.50,
        kerf=(1/8)u"inch",
        # 4ft long, 2ft wide:
        available_stock=[example(AvailablePanel)])
    wanted = 3 * WantedPanel(label="panel",
                             length=1u"ft",
                             width=20u"inch")
    searcher = Searcher(supplier, wanted)
    search(searcher)
    @test length(wanted) == length(searcher.cheapest.finished)
    bought = Set(map(progenitor, searcher.cheapest.finished))
    let
        fname = joinpath(@__DIR__, "fits-in-one.dot")
        dotgraph(fname, makePanelGraph(searcher.cheapest))
        rundot(fname)
        @info(fname)
    end
    @info("bought = $bought")
    @test length(bought) == 1
    @test isempty(unuesd_bought_panels(searcher.cheapest))
end

@testset "Search: buy more" begin
    # Make sure we buy another AvailablePanel when none of our scraps
    # are big enough.
    supplier = Supplier(
        name="test",
        cost_per_cut=1.50,
        kerf=(1/8)u"inch",
        available_stock=[example(AvailablePanel)])
    wanted = 3 * WantedPanel(label="panel",
                             length=1.5u"ft",
                             width=20u"inch")
    searcher = Searcher(supplier, wanted)
    search(searcher)
    @test length(wanted) == length(searcher.cheapest.finished)
    bought = Set(map(progenitor, searcher.cheapest.finished))
    @test length(bought) == 2
    @test isempty(unuesd_bought_panels(searcher.cheapest))
end

# Each element of searcher.cheapest.bought should be the progenitor of
# at leat one searcher.cheapest.finished panel.
function unused_bought_panels(state::SearchState)
    symdiff(Set(state.bought),
            Set(map(progenitor, state.finished)))
end

@testset "Search: too many buys" begin
    box_length = 12.5u"inch"
    box_width = 2u"inch"
    box_height = 1.875u"inch"
    stock_thickness = 0.125u"inch"
    rabbet_depth = 0.0625u"inch"
    compartment_count = 4
    wanted_panels = [
        WantedPanel(
	    label="bottom",
	    length=box_length,
	    width=box_width),
        (2 * WantedPanel(
	    label="side",
	    length=box_length,
	    width=box_height - stock_thickness + rabbet_depth))...,
        ((compartment_count + 1) *
	    WantedPanel(label="divider",
		        width=(box_height - stock_thickness + rabbet_depth),
		        length=(box_width - 2 * stock_thickness + 2 * rabbet_depth))...)
    ]
    supplier = Supplier(
	name = "www.midwestproducts.com",
	kerf = 0.036u"inch",    # the measured set of my Japanese pull saw
	cost_per_cut=0.10,
	available_stock = [
	    AvailablePanel(
		label="basswood 24 × 4 × 1/8",
		length=24u"inch",
		width=4u"inch",
		cost=4)
	])
    searcher = Searcher(supplier, wanted_panels)
    search(searcher)
    @test isempty(unused_bought_panels(searcher.cheapest))
    let
        fname = joinpath(@__DIR__, "too-many-buys.dot")
        dotgraph(fname, makePanelGraph(searcher.cheapest))
        rundot(fname)
        @info(fname)
    end
    bought = Set(map(progenitor, searcher.cheapest.finished))
    @test length(bought) == 2
end

@testset "OrFlipped" begin
    wanted = example(WantedPanel)
    both = orFlipped(wanted)
    @test both[2].was === wanted
    @test both[1].width == both[2].length
    @test both[1].length == both[2].width
    @test both[2].label == "flipped $(wanted.label)"
end

@testset "multiply WantedPanel" begin
    wanted = 2 * example(WantedPanel)
    @test length(wanted) == 2
    @test wanted[1].length == wanted[2].length
    @test wanted[1].width == wanted[2].width
    labels = map(wanted) do p; p.label; end
    index = 1
    while labels[1][index] == labels[2][index]
        index += 1
    end
    @test map(labels) do l; parse(Int, l[index:end]); end == 1:2
end

@testset "multiply orFlipped" begin
    wanted = reduce(vcat,
                    map(orFlipped,
                        2 * example(WantedPanel)))
    @test length(wanted) == 4
end

@testset "Graph" begin
    g = PanelGraph()
    injest(g, :a => :b)
    injest(g, :b => :c)
    injest(g, :a => 1)
    injest(g, 1 => 1.1)
    injest(g, 1 => 1.2)
    @test query(g, :a, Symbol) == Set([:a => :b])
    @test query(g, Symbol, :c) == Set([:b => :c])
    @test length(query(g, Any, Number)) == 3
    @test length(query(g, :b, :c)) == 1
end

include("test_readme_examples.jl")


text:
#  notebook(dir=pwd(),detached=true)

import Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using Revise
# include("jac-main.jl")
using JAC


text:
@testset "Circle/Sphere construction" begin

    for (type, dim_allowed) in zip([Circle, Sphere], [2, 3])
        for dim in 1:4
            point = zeros(Int64, dim)

            if dim == dim_allowed
                obj = type(point, 1)
                @test isa(obj.point, SVector{dim, Int64})
            else
                message = "The dimension of the point must be $dim_allowed"
                @test_throws ArgumentError(message) type(point, 1)
            end
        end
    end
end


text:
# Various utility functions used to make nice plots

# Custom colors
mycolors = [
    PlotThemes.RGB255(49, 114, 174),  # blue
    PlotThemes.RGB255(68, 156, 118),  # green
    PlotThemes.RGB255(117, 112, 173), # purple
    PlotThemes.RGB255(224, 107, 97),  # red
    PlotThemes.RGB255(221, 162, 66),  # yellow
    PlotThemes.RGB255(202, 98, 159),  # magenta
    PlotThemes.RGB255(114, 179, 224), # cyan
]

# Custom theme
mytheme = PlotThemes.PlotTheme(Dict([
            :grid => false,
            :fontfamily => "Computer Modern",
            :label => "",
            :guidefontsize => 12,
            :titlefontsize => 12,
            :legendfontsize => 12,
            :tickfontsize => 8,
            :framestyle => :box,
            :linewidth => 3,
            :palette => mycolors,
            :colorgradient => :viridis,
            :foregroundcolorlegend => nothing,
            :backgroundcolorlegend => nothing,
            ]))

function custom_theme()
    return mytheme
end

function custom_plots()
    add_theme(:mytheme, custom_theme())
    theme(:mytheme)
end

"""
    circle_shape(x, y, r)

Constucts circle object to be used in plots.

Example usage:
plot!(circle_shape(x, y, r), seriestype=[:shape], fillcolor=:red, fillalpha=0.2,
        linecolor=:black, lw=3, edgecolor=:black, label="")
"""
function circle_shape(x, y, r)
    θ = LinRange(0, 2 * π, 500)
    return x .+ r * cos.(θ), y .+ r * sin.(θ)
end

function circle_shape(O::CircularObstacle)
    θ = LinRange(0, 2 * π, 500)
    return O.c[1] .+ O.r * cos.(θ), O.c[2] .+ O.r * sin.(θ)
end


text:
function  problem(bha, traj, wobs, incls;
  penalty=1.0e9, computefunction=PtFEM.p44, pdir=ProjDir)

  prob = Dict{Symbol, Any}(
    :struc_el => Frame(bha[:noofelements], bha[:noofnodes], 2, 1, 1, Line(2, 3)),
    :properties => bha[:properties],
    :etype => bha[:elements].npind,
    :xcoords => Float64.(0:bha[:noofelements]),
    :ycoords => fill(0.0,bha[:noofnodes]),
    :gnum => [(1:bha[:noofelements])'; (2:bha[:noofnodes])'],
    :support => [bha[:noofnodes], (0, 0, 1)],
    :loaded_nodes => [(1, [10000.0 0.0 0.0])],
    :penalty => penalty,
    #:eq_nodal_forces_and_moments => [],
    :traj => traj,
    :wobrange => wobs,
    :inclinationrange => incls,
    :computefunction => computefunction,
    :pdir => pdir
  )

  prob
end

#=
Needs to generate something like:

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(7, 8, 2, 1, 1, Line(2, 3)),
  :properties => [
    1.0e10 1.0e6 20.0;
    1.0e10 1.0e6 50.0;
    1.0e10 1.0e6 80.0
     ],
  :etype => [1, 2, 2, 1, 3, 3, 1],
  :x_coords => [0.0,  0.0, 10.0, 20.0, 20.0, 35.0, 50.0, 50.0],
  :y_coords => [0.0, 15.0, 15.0, 15.0,  0.0, 15.0, 15.0,  0.0],
  :g_num => [
    1 2 3 5 4 6 8;
    2 3 4 4 6 7 7],
  :support => [
    (1, [0 0 0]),
    (5, [0 0 0]),
    (8, [0 0 0])
    ],
  :loaded_nodes => [
    (1, [10000.0 -600.0 -600.0]),
    (2, [0.0 -600.0 0.0]),
    (3, [0.0 -600.0 0.0]),
...
    (10, [0.0 -600.0 0.0]),
    (11, [0.0 -600.0 -600.0])
    ],
  :penalty => 1e19,
  :eq_nodal_forces_and_moments => [
    (1, [0.0 -60.0 -60.0 0.0 -60.0 60.0]),
    (2, [0.0 -120.0 -140.0 0.0 -120.0 140.0]),
    (3, [0.0 -20.0 -6.67 0.0 -20.0 6.67])
  ]
)
=#


text:
# This script trains and saves the 300-40 MLP for the DSGE example
using MXNet
# size of layers
layer1 = 300
layer2 = 40
outputs = 9
include("dataprep.jl") # load training and testing data
# set up data providers
batchsize = 2048 # can adjust this later, but must be defined now for next line
trainprovider = mx.ArrayDataProvider(:data => X, batch_size=batchsize, shuffle=false, :label => Y)
evalprovider = mx.ArrayDataProvider(:data => XT, batch_size=batchsize, shuffle=false, :label => YT)
# set up 2 layer MLP with 2 outputs
data = mx.Variable(:data)
label = mx.Variable(:label)
net  = @mx.chain    mx.FullyConnected(data = data, num_hidden=layer1) =>
                    mx.Activation(act_type=:relu) =>
                    mx.FullyConnected(num_hidden=layer2) =>
                    mx.Activation(act_type=:relu) =>
                    mx.FullyConnected(num_hidden=outputs)
# squared error loss is appropriate for regression, don't change
cost = mx.LinearRegressionOutput(data = net, label=label)
# final model definition, don't change, except if using gpu
model = mx.FeedForward(cost, context=mx.cpu())
# set up the optimizer: select one, explore parameters, if desired
# optimizer = mx.SGD(lr=0.01, momentum=0.9, weight_decay=0.00001)
optimizer = mx.ADAM()
# train, reporting loss for training and evaluation sets
# initial training with small batch size, to get to a good neighborhood
batchsize = 256
mx.fit(model, optimizer, initializer=mx.NormalInitializer(0.0,0.05), eval_metric=mx.MSE(), trainprovider, eval_data=evalprovider, n_epoch = 100)
# more training with larger sample, saving the final fitted model
batchsize = 2048
mx.fit(model, optimizer, eval_metric=mx.MSE(), trainprovider, eval_data=evalprovider, n_epoch = 20, callbacks=[mx.do_checkpoint("dsge30040", frequency=20, save_epoch_0=false)])


text:
using ForwardDiff: derivative, gradient

"""
    OneDimensionalQuantumSystem(β, m, potential)

Creates an 1D quantum system with temperature β and Hamiltonian ``H=p^2/2m+V(x)``.
"""
struct OneDimensionalQuantumSystem
    β::Real
    m::Real
    potential::Function
    force::Function
    hessian::Function
    OneDimensionalQuantumSystem(β::Real, m::Real, potential::Function) = begin
        force = x -> -derivative(potential, x)
        hessian = x -> -derivative(force, x)
        new(β, m, potential, force, hessian)
    end
end

const harmonicPotential = x::Real -> .5x^2
const anharmonicPotential = x::Real -> .5x^2 + .1x^3 + .01x^4
const quarticPotential = x::Real -> .25x^4

const stiffnessMatrix = (system::OneDimensionalQuantumSystem, n::Integer) -> begin
    tmp = zeros(n, n)
    for i in 1:n
        tmp[i,i] += 1
        tmp[mod1(i+1, n), mod1(i+1, n)] += 1
        tmp[i, mod1(i+1, n)] -= 1
        tmp[mod1(i+1, n), i] -= 1
    end
    Symmetric(tmp * system.m * n / system.β^2)
end


text:
import RowEchelon

#=
Homework 2
Section 1.4
Problem 11
=#

# Create the first Matrix
A = [1 3 -4 ;
	 1 5 2 ;
	 -3 -7 6]

# Create the second matrix
B = [-2 ; 4 ; 12]

# Create an augmented matrix
Augment = [1 3 -4 -2 ;
		   1 5 2 4 ;
		   -3 -7 6 12]

# Solve the matrix
Augment_rref = RowEchelon.rref(Augment)

# Print the results
println()
println("Original Matrix, A:")
println(A)
println("Original Matrix, B")
println(B)
println("Augmented Matrix:")
println(Augment)
println("Row Reduced Echelon Matrix:")
println(Augment_rref)
println("Solution as a vector:")
println(Augment_rref[:,4])
println()


text:
using DelimitedFiles
using PyPlot

#rebound = readdlm("rebound_times.txt",',',comments=true)
#nbodygrad = readdlm("nbodygradient_times.txt",',')
#nbodygrad_pair = readdlm("nbodygradient_times_pair.txt",',')
using DelimitedFiles
#rebound = readdlm("reb_times.txt",',')
rebound_ias15 = readdlm("rebound_times_ias15.txt",',',comments=true)
rebound_whfast = readdlm("rebound_times_whfast.txt",',',comments=true)
nbodygrad = readdlm("nbg_times.txt",',',comments=true)

using PyPlot

clf()
semilogy(rebound_ias15[:,1],rebound_ias15[:,3],label="IAS15, gradient",c="r",linewidth=2)
semilogy(rebound_ias15[:,1],rebound_ias15[:,3],".",c="r",markersize=15)
semilogy(nbodygrad[:,1],nbodygrad[:,6],label="NBG, tt w/ gradient",c="b",linewidth=2,linestyle=":")
semilogy(nbodygrad[:,1],nbodygrad[:,6],".",c="b",markersize=15)
semilogy(nbodygrad[:,1],nbodygrad[:,3],label="NBG, gradient",c="b",linewidth=2)
semilogy(nbodygrad[:,1],nbodygrad[:,3],".",c="b",markersize=15)
semilogy(nbodygrad[:,1],nbodygrad[:,5],label="NBG, gradient, adjacent Kepler",c="g",linewidth=2)
semilogy(nbodygrad[:,1],nbodygrad[:,5],".",c="g",markersize=15)
semilogy(rebound_whfast[:,1],rebound_whfast[:,3],label="WHFAST, gradient",c="m",linewidth=2)
semilogy(rebound_whfast[:,1],rebound_whfast[:,3],".",c="m",markersize=15)

semilogy(rebound_ias15[:,1],rebound_ias15[:,2],label="IAS15, no gradient",linestyle="dashed",c="r",linewidth=2)
semilogy(rebound_ias15[:,1],rebound_ias15[:,2],".",c="r",markersize=15)
semilogy(nbodygrad[:,1],nbodygrad[:,2],label="NBG, no gradient",linestyle="dashed",c="b",linewidth=2)
semilogy(nbodygrad[:,1],nbodygrad[:,2],".",c="b",markersize=15)
semilogy(nbodygrad[:,1],nbodygrad[:,4],label="NBG, no gradient, adjacent Kepler",linestyle="dashed",c="g",linewidth=2)
semilogy(nbodygrad[:,1],nbodygrad[:,4],".",c="g",markersize=15)
semilogy(rebound_whfast[:,1],rebound_whfast[:,2],label="WHFAST, no gradient",linestyle="dashed",c="m",linewidth=2)
semilogy(rebound_whfast[:,1],rebound_whfast[:,2],".",c="m",markersize=15)
grid(linestyle=":")
legend(loc="upper_left",fontsize=8)
axis([0,11,3e-3,1e2])
xlabel("Number of planets",fontsize=20)
ylabel("Integration time [sec]",fontsize=20)
#minorticks_on()
#tick_params(labeltop="true",labelright="true")
tick_params(axis="both",direction="inout",which="both",right="true",top="true",labeltop="true",labelright="true")
#secondary_yaxis("right")
#PyPlot.savefig("rebound_vs_nbodygradient.pdf",bbox_inches="tight")
PyPlot.savefig("reb_vs_nbg.pdf", bbox_inches="tight")

max_forks_repo_name: langfzac/nbg-papers


text:
using Winston

println("-----------------------------------")
println("  Welcome to G(t) Plot.jl  ")
println("-----------------------------------")

mt_int = ARGS[1]
mt_float = ARGS[2]
xi = ARGS[3]

Data = readcsv("Data/Gauge_$(mt_int)_$(mt_float)_$(xi).csv")

t = Data[:,1];
# λ = Data[:,2];
# yt = Data[:,3];
# g1 = Data[:,4];
# g2 = Data[:,5];
# g3 = Data[:,6];
G = Data[:,7];

# Gauge Plot
p = FramedPlot(
    title="G(t) Plots",
    xlabel="t",
    ylabel="Gauge");
C = Curve(t, G)
setattr(C, "label", "G(t)")
lgnd = Legend(.9, .9, [C])
add(p, C, lgnd)
savefig(p, "Fig/G_$(mt_int)_$(mt_float)_$(xi).svg", (1000, 600))
run(`inkscape -z Fig/G_$(mt_int)_$(mt_float)_$(xi).svg -e Fig/G_$(mt_int)_$(mt_float)_$(xi).png -d 300 --export-background=WHITE`)
run(`rm Fig/G_$(mt_int)_$(mt_float)_$(xi).svg`)

max_stars_repo_name: Axect/RGE

max_issues_repo_name: Axect/RGE

max_forks_repo_name: Axect/RGE


text:
using ScRNAseq
using DelimitedFiles
cd(dirname(@__FILE__))

X=readdlm("X.txt",',',Int16)
genelist=vec(readdlm("genelist.txt",String))

X,genelist=ScRNAseq.QualityControl.selectg(X,genelist)
X1=ScRNAseq.Transformation.pearsonresiduals(X)
X2=ScRNAseq.Normalization.norm_libsize(X)
Y=ScRNAseq.Embedding.umap(X1)

# Y2=ScRNAseq.Embedding.tsne(X)
using Plots
scatter3d(Y[:,1],Y[:,2],Y[:,3])


text:
#
#  Initial values from:
#         "Explicit symplectic approximation of nonseparable
#          Hamiltonians: Algorithm and long time performance"
#                 Molei Tao
#

function InitialNLS(T=Float64)
#
# N=5
#

    q=[parse(BigFloat,"3.0"), parse(BigFloat,"0.01"), parse(BigFloat,"0.01"),
       parse(BigFloat,"0.01"), parse(BigFloat,"0.01")]

    p=[parse(BigFloat,"1.0"), parse(BigFloat,"0.0"), parse(BigFloat,"0.0"),
       parse(BigFloat,"0.0"), parse(BigFloat,"0.0")]

    u0=convert.(T,[q; p])

    return u0

end


text:
module Day11
export day11

_string_iterator(str) = Iterators.take(str, length(str))

parseLine(line) = parse.(Int64, _string_iterator(line))

function readInput(inputfile="inputs/input11.txt")
    return hcat([parseLine(line) for line in eachline(inputfile)]...)
end

const ADJACENT = let temp = vec(CartesianIndex.(-1:1, (-1:1)'))
    deleteat!(temp, 5)
    temp
end # delete (0,0)

function step!(grid)
    grid .+= 1

    flashed = 0
    for I in CartesianIndices(grid)
        if grid[I] > 9
            flashed += _octopus_flash!(grid, I)
        end
    end
    return flashed
end

function _octopus_flash!(grid, loc)
    grid[loc] = 0
    flashed = 1
    for adj in ADJACENT
        neighbour = loc + adj
        checkbounds(Bool, grid, neighbour) || continue
        if grid[neighbour] > 0 # exclude already flashed
            grid[neighbour] += 1
        end
        if grid[neighbour] > 9
            flashed += _octopus_flash!(grid, neighbour)
        end
    end
    return flashed
end

function day11Part1!(grid, steps=100)
    #return sum(step!(grid) for _ in 1:steps) # less allocs but not really faster
    flashes = [step!(grid) for _ in 1:steps]
    @assert !any(flashes .== length(grid))
    return sum(flashes)
end

function day11Part2!(grid)
    flashed = 0
    round = 0 # from part 1
    while flashed != length(grid)
        flashed = step!(grid)
        round += 1
    end
    return round
end

function day11(inputdata = readInput("inputs/input11.txt"))
    return (day11Part1!(inputdata, 100), 100+day11Part2!(inputdata))
end

end #module


text:
"""
    nearestneighbor(x, y)

Nearest neighbor interpolation with the spline values in the center of intervals.

- aliases: `spline0`, `stepspline`
"""
function nearestneighbor(x::Vector, y::Vector)
    NearestNeighborSpline(y, x)
end

const spline0 = nearestneighbor
const stepspline = nearestneighbor
const nspline = nearestneighbor

"""
    linearspline(x, y)

Linear spline interpolation.

- aliases: `spline1`, `lspline`
"""
function linearspline(x::Vector, y::Vector)
    #n = length(x)
    b = Vector{Float64}(undef, length(x))

    @inbounds for i in 1:length(x)-1
        b[i] = (y[i+1] - y[i]) / (x[i+1] - x[i])
    end
    #yy = y[begin+1:end].-y[begin:end-1]
    #xx = x[begin+1:end].-x[begin:end-1]
    LinearSpline(b, y, x)
end

const spline1 = linearspline
const lspline = linearspline

"""
    quadraticspline(x::Vector, y::Vector, b1 = 0.0)

Quadratic spline interpolation.

- aliases: `spline2`, `qspline`

Optional argument `b1` sets the 1st derivative at the first knot. Default value = `0.0`.
"""
function quadraticspline(x::Vector, y::Vector, b1 = 0.0)
    n = length(x)
    b = Vector{Float64}(undef, n-1)
    c = Vector{Float64}(undef, n-1)
    b[1] = b1

    Δx = Vector{Float64}(undef, n-1)
    Δy = Vector{Float64}(undef, n-1)

    @inbounds for i in 1:n-1
        Δx[i] = x[i+1] - x[i]
        Δy[i] = y[i+1] - y[i]
    end

    @inbounds for i in 1:n-1
        b[i+1] = -b[i] + 2 * Δy[i] / Δx[i]
    end

    @inbounds for i in 1:n-1
         c[i] = (b[i+1] - b[i]) / (2 * Δx[i])
    end
    QuadraticSpline(c, b, y, x)
end

const spline2 = quadraticspline
const qspline = quadraticspline

"""
    cubicspline(x::Vector, y::Vector, style...)

- alias: `spline3`, `cspline`

Creates a cubic spline. Boundary conditions are defined by `style`.

`style...`: settings for boundary conditions.

- `:natural`: zero curvature at the start and end points of the spline.
    - aliases: `cubicspline(x, y)`, `naturalspline(x, y)`
- `:clamped, ds, de`: predefined 1st derivative of start (ds) and end points (de).
    - alias: `clampedspline(x, y, ds, de)`
- `:periodic`: periodic boundary conditions. `y` values at start and end points have to be equal,
conditions at start and end points are identical.
    - alias: `periodicspline(x, y)`
- `:parabolic`, `:notaknot`: parabolic runout spline, also known as not-a-knot spline
- `:cubic`: cubic runout spline
- `:constrained`: constrained cubic spline, see Kruger’s paper
- `:hermite`: hermite spline
    - `:kruger`: constrained cubic spline, see Kruger’s paper. No overshooting.
    - `:free`: hermite spline with free derivative conditions. No overshooting.
    - `:predef`: predefined derivatives at knots

# Examples:
```jldoctest
julia> x=[0,1,2,3,4]; y=[-1,1,5,2,-3]; spline3(x,y,:natural)
CubicSpline([1.0, -3.0, 2.0, 0.0], [0.0, 3.0, -6.0, 0.0], [1.0, 4.0, 1.0, -5.0], [-1.0, 1.0, 5.0, 2.0])
```
"""
function cubicspline(x::Vector, y::Vector, style...)
    if style[1] == :periodic
        y[1] != y[end] ? error("Periodic spline error: different y values at endpoints!") : nothing
    elseif style[1] == :constrained
        return hermitespline(x, y, :kruger)
    elseif style[1] == :hermite
        if length(style) == 1
            return hermitespline(x, y, :kruger)
        elseif length(style) == 2
            return hermitespline(x, y, style[2])
        else
            throw(DomainError(style[1], "Spline style does not exist."))
        end
    else
        n = length(x)
        #Δx= Vector{Float64}(undef, n-1)
        Δx = zeros(Float64, n-1)
        #Δy= Vector{Float64}(undef, n-1)
        Δy = zeros(Float64, n-1)
        #λ = Vector{Float64}(undef, n)    # upper diagonal
        λ = zeros(Float64, n)
        δ = ones(n) * 2 # main diagonal
        #μ = Vector{Float64}(undef, n)    # lower diagonal
        μ = zeros(Float64, n)
        #r = Vector{Float64}(undef, n)    # derivatives
        r = zeros(Float64, n)
        #M = Vector{Float64}(undef, n)
        #M = zeros(Float64, n)

        @inbounds for i in 1:n-1
            Δx[i] = x[i+1] - x[i]
            Δy[i] = y[i+1] - y[i]
        end
        if style == :notaknot
            A = zeros(n, n)
            #=@inbounds=# for i in 2:n-1
                A[i, i-1] = Δx[i-1]
                A[i, i] = 2.0 * (Δx[i-1] + Δx[i])
                A[i, i+1] = Δx[i+1]
                r[i] = 6.0 * (Δy[i] / Δx[i] - Δy[i-1] / Δx[i-1])
            end
            #A[begin, begin] = -Δx[begin + 2]
            #A[begin, begin + 1] = Δx[begin + 1] + Δx[begin + 2]
            #A[begin, begin + 2] = -Δx[begin + 1]
            A[1,1] = -(x[3] - x[2])
            A[1,2] = x[3] - x[1]
            A[1,3] = -(x[2] - x[1])
            #A[end, end - 2] = -Δx[end - 1]
            #A[end, end - 1] = Δx[end - 2] + Δx[end - 1]
            #A[end, end]   = Δx[end - 2]
            A[n-1, n-3] = -(x[n-1] - x[n-2])
            A[n-1, n-2] = x[n-1] - x[n-3]
            A[n-1, n-1] = -(x[n-2] - x[n-3])
            #A[n-1, n-1] = -(x[n-3] - x[n-2])
            #A[n-1, n-2] = x[n-3] - x[n-1]
            #A[n-1, n-3] = -(x[n-2] - x[n-1])
            #for i in 2:n-1
            #    r[i] = 6.0 / (Δx[i-1] + Δx[i]) * (Δy[i] / Δx[i] - Δy[i-1] / Δx[i-1])
            #    #r[i] = 6.0 * (Δy[i] / Δx[i] - Δy[i-1] / Δx[i-1])
            #end
        else
            #=@inbounds=# for i in 2:n-1
                λ[i] = Δx[i] / (Δx[i-1] + Δx[i])
                μ[i] = 1.0 - λ[i]
                r[i] = (6.0 / (Δx[i-1] + Δx[i])) * (Δy[i] / Δx[i] - Δy[i-1] / Δx[i-1])
            end
        end

        e = length(Δx)
        r = rbound(r, e, Δx, Δy, style...)
        if style ≠ :notaknot
            λ, μ = λµbound(λ, μ, e, Δx, style...)
            A = Tridiagonal(μ[2:n], δ, λ[1:n-1])
            M = A \ r
            M = Mbound(M,style...)
        else
            M = A \ r
        end
        b = zeros(Float64, n-1) # Vector{Float64}(undef, n-1)
        c = zeros(Float64, n-1) # Vector{Float64}(undef, n-1)
        d = zeros(Float64, n-1) # Vector{Float64}(undef, n-1)

        #=@inbounds=# for i in 1:n-1
            d[i] = (M[i+1] - M[i]) / (6 * Δx[i])
            c[i] = 0.5 * M[i] #
            b[i] = Δy[i] / Δx[i] - Δx[i] / 6.0 * (2.0 * M[i] + M[i+1])
        end
    end
    if style[1] in (:hermite, :constrained)
        nothing
    else
        return CubicSpline(d, c, b, y, x)
    end
end

const spline3 = cubicspline
const cspline = cubicspline

"""
    rbound(r,Δx,Δy,style...)

Returns vector `r` according to boundary conditions defined by `style...`
Helper function for `cubicspline`.

Returns: `r`
"""
function rbound(r, e, Δx, Δy, style...)
    if style[1] in (:natural, :notaknot)
        r[1]   = 0.0
        r[end] = 0.0
    elseif style[1] == :parabolic
        r[1] = r[2]
        r[end] = r[end-1]
    elseif style[1] == :clamped
        ds     = style[2]
        de     = style[3]
        r[1]   = 6.0 / Δx[1] * (Δy[1] / Δx[1] - ds)    #ds = start derivative
        r[end] = -6.0 / Δx[e] * (Δy[e] / Δx[e] - de)   #de = end derivative
    elseif style[1] == :periodic
        r[end] = 6.0 / (Δx[1] + Δx[e]) * (Δy[1] / Δx[1] - Δy[e] / Δx[e])
        r[1]   = r[end]
    end
    return r
end

"""
    Mbound(M,style...)

Returns vector `M` (moments) according to boundary conditions defined by `style...`

Helper function for `cubicspline`.
"""
function Mbound(M, style...)
    if style[1] == :natural
        M[1]   = 0.0
        M[end] = 0.0
    elseif style[1] == :periodic
        M[1] = M[end]
    elseif style[1] == :parabolic
        M[1] = M[2]
        M[end] = M[end-1]
    elseif style[1] == :notaknot
        #M[1] = M[2]
        #M[end] = M[end-1]
    end
    return M
end

function λµbound(λ,μ,e,Δx,style...)
    if style[1] == :natural
        λ[1]   = 0.0
        μ[end] = 0.0
    elseif style[1] == :clamped
        λ[1]   = 1.0
        μ[end] = 1.0
    elseif style[1] == :parabolic
        λ[1] = λ[2]
        μ[end] = λ[end-1]
    elseif style[1] == :periodic
        λ[end] = Δx[1] / (Δx[1] + Δx[e-1])
        μ[end] = 1 - Δx[1] / (Δx[1] + Δx[e])
    end
    return λ, µ
end

"""
    cubicspline(x, y)
Fall back function, setting the default boundary condition to `:natural`

- alias: `naturalspline`
"""
function cubicspline(x, y)
    cubicspline(x, y, :natural)
end

"""
    naturalspline(x::Vector, y::Vector)

Equivalent to `cubicspline(x::Vector, y::Vector, :natural)`.
"""
naturalspline(x, y) = cubicspline(x, y, :natural)

"""
    clampedspline(x::Vector, y::Vector, ds = 0.0, de = 0.0)

Equivalent to `cubicspline(x::Vector, y::Vector, ds, de)`.

Returns a clamped cubic spline.

Default arguments: 1st and last derivatives = `0.0`
"""
function clampedspline(x, y, ds = 0.0, de = 0.0)
    cubicspline(x, y, :clamped, ds, de)
end

"""
    periodicspline(x::Vector, y::Vector)

Equivalent to `cubicspline(x::Vector, y::Vector, :periodic)`.

`y` values of start and end points have to be equal.
"""
function periodicspline(x, y)
    cubicspline(x, y, :periodic)
end

notaknotspline(x, y) = cubicspline(x, y, :notaknot)

"""
    hermitespline(x::Vector, y::Vector, style::Symbol = :kruger, args...)

Creates cubic hermite interpolation spline with defined 1st derivatives at all knots.

Allowed `style`s:
- `:kruger`: Default value. Defines 1st derivatives using harmonic means to prevent overshoot at extremal knots.
- `:free`: Simple mean values for 1st derivatives.
- `:predef`: user defined 1st derivatives at all knots.
"""
function hermitespline(x::Vector, y::Vector, style::Symbol = :kruger, args...)
    n = length(x)
    # 1st derivatives
    b = Vector{Float64}(undef, n)
    # finite differences
    Δx = Vector{Float64}(undef, n-1)
    Δy = Vector{Float64}(undef, n-1)

    @inbounds for i in 1:n-1
        Δx[i] = x[i+1] - x[i]
        Δy[i] = y[i+1] - y[i]
    end

    #generate 1st derivatives at kntos

    # calculate derivatives based on C.J.C. Kruger,
    # [`Constrained Cubic Spline Interpolation for Chemical Engineering Applications`](https://pages.uoregon.edu/dgavin/software/spline.pdf), 2002
    # which utilizes the harmonic mean of the neighboring slopes of each knot.
    if style == :kruger
        @inbounds for i in 2:n-1
            if Δy[i-1] * Δy[i] > 0.0
                b[i] = 2.0 / (Δx[i] / Δy[i] + Δx[i-1] / Δy[i-1]) #harmonic mean
            end
        end
        # boundary knots
        b[1] = 1.5 * Δy[1]/Δx[1] - b[2] / 2
        b[n] = 1.5 * Δy[n-1]/Δx[n-1] - b[n-1] / 2

    # 1st derivatives at knots are arithmetic means of neighboring slopes
    elseif style == :free
        @inbounds for i in 2:n-1
            b[i] = Δy[i] / Δx[i]
        end
        # boundary knots
        b[1] = Δy[1] / Δx[1]
        b[n] = Δy[n-1] / Δx[n-1]
    elseif style == :predef
        if length(args) == 0
            throw(ArgumentError("Derivative vector missing."))
        end
        if length(args[1]) ≠ n-1
            throw(DimensionMismatch("x, y, and derivative vectors must have the same length"))
        end
        @inbounds for i in 1:length(args[1])
            b[i] = args[1][i]
        end
    else
        throw(DomainError(style, "Derivative constraint style does not exist!"))
    end

    HermiteSpline(b, y, x)
end

const hspline = hermitespline

#    for i in 1:n-1
#        Φ0 = 2 * Δx[i]^3 - 3 * Δx[i]^2 + 1
#        Φ1 = Δx[i]^3 - 2 * Δx[i]^2 + Δx[i]
#        Ψ0 = -2 * Δx[i]^3 + 3 * Δx[i]^3
#        Ψ1 = Δx[i]^3 - Δx[i]^2
#
#        Φ0_ = 6 * Δx[i]^2 - 6 * Δx[i]
#        Φ1_ = 3 * Δx[i]^2 - 4 * Δx[i] + 1
#        Ψ0_ = -6 * Δx[i]^2 + 6 * Δx[i]
#        Ψ1_ = 3 * Δx[i]^2 - 2 * Δx[i]
#
#        Φ0__ = 12 * Δx[i] - 6
#        Φ1__ = 6 * Δx[i] - 4
#        Ψ0__ = -12 * Δx[i] + 6
#        Ψ1__ = 6 * Δx[i] - 2
#
#        Φ0___ = 12
#        Φ1___ = 6
#        Ψ0___ = -12
#        Ψ1___ = 6
#        #xb = (xs[i] - spl.x[j]) / (spl.x[j+1] - spl.x[j])
#        c[i] = (y[i] * Φ0_ + Δx[i] * b[i] * Φ1_ + y[i+1] * Ψ0_ + Δx[i] * b[i+1] * Ψ1_) / Δx[i]
#        #d[i] = (y[i] * Φ0__ + Δx[i] * b[i] * Φ1__ + y[i+1] * Ψ0__ + Δx[i] * b[i+1] * Ψ1__) / Δx[i]
#    end

    #CubicSpline(d, c, b, y, x)

"""
    nakspline(x::Vector, y::Vector)

Creates not-a-knot spline.
"""
function nakspline(x::Vector, y::Vector)
    n = length(x)
    Δx = zeros(Float64, n)
    Δy = zeros(Float64, n)
    r = zeros(Float64, n-1)
    for i in 1:n-1
        Δx[i] = x[i+1] - x[i]
        Δy[i] = y[i+1] - y[i]
    end
    A = zeros(n-1, n-1)

    for i in 2:n-2
        A[i,i-1] = Δx[i]
        A[i,i]   = 2 * (Δx[i] + Δx[i+1])
        A[i,i+1] = Δx[i+1]
        r[i] = 3 * Δy[i+1] / Δx[i+1] - 3 * Δy[i] / Δx[i]
    end

    A[1,1] = Δx[2]
    A[1,2] = -(Δx[1] + Δx[2])
    A[1,3] = Δx[1]

    A[n-1, n-3] = Δx[n-1]
    A[n-1, n-2] = -(Δx[n-2] + Δx[n-1])
    A[n-1, n-1] = Δx[n-2]

    r[1] = r[2]
    r[n-1] = r[n-2]

    #return A,r

    M = A \ r

    b = zeros(Float64, n-1) # Vector{Float64}(undef, n-1)
    c = zeros(Float64, n-1) # Vector{Float64}(undef, n-1)
    d = zeros(Float64, n-1) # Vector{Float64}(undef, n-1)
    for i in 1:n-2
        #d[i] = (M[i+1] - M[i]) / (6 * Δx[i])
        d[i] = (M[i+1] - M[i]) / (3 * Δx[i])
        c[i] = 0.5 * M[i] #
        #b[i] = Δy[i] / Δx[i] - Δx[i] / 6.0 * (2.0 * M[i] + M[i+1])
        b[i] = Δy[i] / Δx[i] - (M[i+1] + 2 * M[i]) * Δx[i] / 3
    end
    return CubicSpline(d, c, b, y, x)
end


text:
using StanSample, MCMCChains, Distributions

# Define the Stan language model

m8_4s = "
data{
  int N;
  vector[N] y;
}
parameters{
  real sigma;
  real alpha;
  real a1;
  real a2;
}
model{
  real mu;
  mu = a1 + a2;
  y ~ normal( mu , sigma );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m_8_4s = SampleModel("m8.4s", m8_4s);

# Input data for cmdstan

m8_4_data = Dict("N" => 100, "y" => rand(Normal(0, 1), 100));

# Sample using cmdstan

rc = stan_sample(m_8_4s, data=m8_4_data);

rethinking = "
        mean   sd  5.5% 94.5% n_eff Rhat
alpha 0.06 1.90 -2.22  2.49  1321    1
sigma 2.15 2.32  0.70  5.21   461    1
";

# Describe the draws

if success(rc)
  chn = read_samples(m_8_4s; output_format=:mcmcchains)
  describe(chn)
end


text:
using Folds.Testing: parse_tests, test_with_sequential
using FoldsKernelAbstractions

# Taken from FoldsCUDA.jl
rawdata = """
count(isodd(x) for x in 1:10)
extrema((x - 5)^2 for x in 1:10)
mapreduce(x -> x^2, +, 1:10)
mapreduce(*, +, 1:10, 11:20)
mapreduce(*, +, 1:10, 11:20, 21:30)
maximum(0:9)
maximum(9:-1:0)
maximum([2, 3, 0, 3, 4, 0, 5, 7, 4, 2])
minimum(0:9)
minimum(9:-1:0)
minimum([2, 3, 0, 3, 4, 0, 5, 7, 4, 2])
prod(1:2:10)
sum(1:10)
sum(x^2 for x in 1:11)
sum(x^2 for x in 1:11 if isodd(x); init = 0)
"""

# Works in FoldsCUDA but not with CUDADevice:
#=
prod(y for x in 1:11 if isodd(x) for y in 1:x:x^2; init = 1)
=#

tests = parse_tests(rawdata, @__MODULE__)

get_executors(device) = [
    KAEx(device; groupsize = 1, basesize = 1),
    KAEx(device; groupsize = 1, basesize = 3),
    KAEx(device; groupsize = 3, basesize = 1),
    KAEx(device; groupsize = 3, basesize = 3),
]


text:
using Test, Flags

@testset "Flags.jl" begin

f = Flag("a", 1, "hello world")
@test formathelp(f) == "  -a=1: hello world"

end