text:
using Statistics
using Profile

@enum BinInit begin
    LightSpeed  # set limits of initial bin +- 3.1e8 m/s
    MinMaxVel  # set limits of initial bin to min/max velocities of particles in question
end

@enum BinSplit begin
    MiddleSplit  # split bin among middle
    MeanVelSplit  # split bin along mean velocity in bin
end

@enum BinRefinement begin
    Greedy  # refine the first bin that is suitable
    Heaviest  # refine the heaviest bin
end

struct Constants
    k::Float64
end

constants = Constants(1.38064852e-23)

mutable struct Particle
    w::Float64
    v::Array{Float64,1}
    octant::Int
end

mutable struct Bin
    w::Float64  # total weight in bin
    nparticles::Int  # total number of particles in bin
    v::Array{Float64,1}  # velocity of bin
    c::Array{Float64,1}  # second moment in bin

    v_lo::Array{Float64,1}  # min velocity of bin
    v_hi::Array{Float64,1}  # max velocity of bin
    v_split::Array{Float64,1}  # velocities along which octant splitting is performed

    index_start::Int
    index_end::Int # indices of particle array - start, end of particles belonging to this bin

    level::Int
end

function maxwellian(T, mass)
    # sample 3-D velocity from Maxwellian
    return sqrt(constants.k * T / mass) * randn((3,))
end

function compute_octant!(particle, v_split)
    # compute index of octree bin particle is in
    particle.octant = 1
    if particle.v[1] > v_split[1]
        particle.octant += 1
    end
    if particle.v[2] > v_split[2]
        particle.octant += 2
    end
    if particle.v[3] > v_split[3]
        particle.octant += 4
    end
end

function compute_bin_values!(particle_list, bin)
    # compute density, velocity and 2nd order moments in bin
    bin.w = 0.0
    bin.nparticles = bin.index_end - bin.index_start + 1
    bin.v = [0.0, 0.0, 0.0]
    bin.c = [0.0, 0.0, 0.0]

    for i=bin.index_start:bin.index_end
        bin.w += particle_list[i].w
        bin.v += particle_list[i].w .* particle_list[i].v
    end

    if (bin.w > 0.0)
        bin.v = bin.v / bin.w
    end

    for i=bin.index_start:bin.index_end
        bin.c += particle_list[i].w * ((particle_list[i].v .- bin.v).^2)
    end

    if (bin.w > 0.0)
        bin.c = sqrt.(bin.c / bin.w)
    end
end

function compute_bin_weight_and_split!(particle_list, bin, split_type)
    # when refining, we just need to keep track of the weight and number of particles in a bin
    # also computes the velocities to split the bin along

    bin.w = 0.0
    bin.nparticles = bin.index_end - bin.index_start + 1
    bin.v = [0.0, 0.0, 0.0]

    for i=bin.index_start:bin.index_end
        bin.w += particle_list[i].w
        bin.v += particle_list[i].w .* particle_list[i].v
    end

    if ((split_type == MiddleSplit) || (bin.w <= 0.0))
        bin.v_split = (bin.v_lo + bin.v_hi) / 2
    else
        bin.v_split = bin.v / bin.w
    end
end

function refine_bin!(particle_list, bins_list, bin_id, n_bins, split_type)
    # in order to refine a bin, we need to keep particle sorted
    # by octant
    # we use a counting sort
    bin_counts = [0, 0, 0, 0, 0, 0, 0, 0]

    for particle in particle_list[bins_list[bin_id].index_start:bins_list[bin_id].index_end]
        compute_octant!(particle, bins_list[bin_id].v_split)
    end

    # a better/faster version would operate on lists of pointers to particles, not particle structs themselves
    particles_in_bin = particle_list[bins_list[bin_id].index_start:bins_list[bin_id].index_end]

    for particle in particle_list[bins_list[bin_id].index_start:bins_list[bin_id].index_end]
        bin_counts[particle.octant] += 1
    end

    bin_counts_copy = bin_counts[:]

    for i=2:8
        bin_counts[i] += bin_counts[i-1]
    end

    bin_counts_acc_copy = bin_counts[:]

    for particle in particle_list[bins_list[bin_id].index_start:bins_list[bin_id].index_end]
        particles_in_bin[bin_counts[particle.octant]] = particle
        bin_counts[particle.octant] -= 1
    end

    particle_list[bins_list[bin_id].index_start:bins_list[bin_id].index_end] = particles_in_bin[:]

    new_bins = 0

    for i=1:8
        if bin_counts_copy[i] > 0
            new_bins += 1

            # set new bin limits, split along z axis
            if i<=4
                vz_min = bins_list[bin_id].v_lo[3]
                vz_max = bins_list[bin_id].v_split[3]
            else
                vz_min = bins_list[bin_id].v_split[3]
                vz_max = bins_list[bin_id].v_hi[3]
            end

            # set new bin limits, split along y axis
            if i%2 == 0
                vy_min = bins_list[bin_id].v_split[2]
                vy_max = bins_list[bin_id].v_hi[2]
            else
                # lower vy values
                vy_min = bins_list[bin_id].v_lo[2]
                vy_max = bins_list[bin_id].v_split[2]
            end

            if ((i==3) || (i==4) || (i==7) || (i==8))
                vx_min = bins_list[bin_id].v_split[1]
                vx_max = bins_list[bin_id].v_hi[1]
            else
                vx_min = bins_list[bin_id].v_lo[1]
                vx_max = bins_list[bin_id].v_split[1]
            end

            if i==1
                istart = bins_list[bin_id].index_start
            else
                istart = bin_counts_acc_copy[i-1] + bins_list[bin_id].index_start
            end
            iend = bin_counts_acc_copy[i] + bins_list[bin_id].index_start - 1

            # split are zeros for now, we will compute them
            push!(bins_list, Bin(0.0, 0, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0],
                                 [vx_min, vy_min, vz_min], [vx_max, vy_max, vz_max], [0.0, 0.0, 0.0],
                                 istart, iend,
                                 bins_list[bin_id].level+1))

            compute_bin_weight_and_split!(particle_list, bins_list[n_bins+new_bins], split_type)
        end
    end

    bins_list[bin_id:n_bins + new_bins-1] = bins_list[bin_id+1:n_bins + new_bins]

    resize!(bins_list, n_bins + new_bins - 1)
    # shift bins to left to overwrite bin
    return n_bins + new_bins - 1
end

function init_bin(particle_list, bin_init, split_type)
    # create initial bin and set initial velocity splitting values
    vx_min = 3.1e8
    vx_max = -3.1e8
    vy_min = 3.1e8
    vy_max = -3.1e8
    vz_min = 3.1e8
    vz_max = -3.1e8

    if bin_init == LightSpeed
        print("initial bounds set to +- 1.03c\n")
    else
        print("initial bounds set to min-max velocities\n")

        for particle in particle_list
            if particle.v[1] < vx_min
                vx_min = particle.v[1]
            elseif particle.v[1] > vx_max
                vx_max = particle.v[1]
            end

            if particle.v[2] < vy_min
                vy_min = particle.v[2]
            elseif particle.v[2] > vy_max
                vy_max = particle.v[2]
            end

            if particle.v[3] < vz_min
                vz_min = particle.v[3]
            elseif particle.v[3] > vz_max
                vz_max = particle.v[3]
            end
        end
    end

    # compute bin weight
    bin = Bin(0.0, 0, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0],
              [vx_min, vy_min, vz_min], [vx_max, vy_max, vz_max], [0.0, 0.0, 0.0],
              1, length(particle_list), 0)
    compute_bin_weight_and_split!(particle_list, bin, split_type)

    # should not merge particles with velocities of different signs
    bin.v_split = [0.0, 0.0, 0.0]
    return bin
end

function octree_merge!(particle_list, bin_list, target_particle_count, max_refinement, aggressiveness, split_type, bin_refinement)
    if (length(particle_list) > target_particle_count)
        do_refinement = true
        target_weight = bin_list[1].w / target_particle_count

        while (do_refinement)
            curr_num_particles = length(bin_list)*2
            if (curr_num_particles >= target_particle_count)
                do_refinement = false
            else
                maxweight = 0.0
                refine_bin_id = -1

                if (bin_refinement == Greedy)
                # based on octree paper, attempt to refine bins greedily
                    for (i, bin) in enumerate(bin_list)
                        # check that we haven't reached maximum recursion level
                        if (bin.level < max_refinement)
                            if ((bin.nparticles > 2) && (bin.w > aggressiveness * target_weight))
                                refine_bin_id = i
                                break
                            end
                        end
                    end
                else
                # refine heaviest bin first
                    for (i, bin) in enumerate(bin_list)
                        # check that we haven't reached maximum recursion level
                        if ((bin.level < max_refinement) && (bin.w > maxweight))
                            if ((bin.nparticles > 2) && (bin.w > aggressiveness * maxweight))
                                refine_bin_id = i
                                maxweight = bin.w
                            end
                        end
                    end
                end
                if (refine_bin_id == -1)
                    print("No bins found to refine")
                    do_refinement = false
                else
                    refine_bin!(particle_list, bin_list, refine_bin_id, length(bin_list), split_type)
                end
            end
        end
    end

    for bin in bin_list
        compute_bin_values!(particle_list, bin)
    end
end

function compute_cell!(particle, N_cell, v_lo, v_hi, v_step)
    # for a cell-based merge (assuming N_cell is the same in the x/y/z directions)
    # compute index of cell particle is in
    # if particle is outside the cell region, the index is based on the octant it is in
    if ((particle.v[1] > v_lo[1]) && (particle.v[1] < v_hi[1]) &&
        (particle.v[2] > v_lo[2]) && (particle.v[2] < v_hi[2]) &&
        (particle.v[3] > v_lo[3]) && (particle.v[3] < v_hi[3]))

        i_x = (particle.v[1] - v_lo[1]) / v_step[1]
        i_y = (particle.v[2] - v_lo[2]) / v_step[2]
        i_z = (particle.v[3] - v_lo[3]) / v_step[3]

        particle.octant = trunc(Int, i_x) * (N_cell^2) + trunc(Int, i_y) * N_cell + trunc(Int, i_z)

        particle.octant += 1
    else
        particle.octant = N_cell^3 + 1

        if (particle.v[1] > 0.0)
            particle.octant += 1
        end
        if (particle.v[2] > 0.0)
            particle.octant += 2
        end
        if (particle.v[3] > 0.0)
            particle.octant += 4
        end
    end
end

function cell_merge!(particle_list, N_cell, v_lo, v_hi)
    # merge all particles in a cell down to 2 particles

    bin_list = Array{Bin}(undef, N_cell^3+8)
    v_step = (v_hi - v_lo) / N_cell

    for i in 1: N_cell^3+8
        bin_list[i] = Bin(0.0, 0, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0],
                          [0.0, 0.0, 0.0], 0, 0, 0)
    end

    for particle in particle_list
        compute_cell!(particle, N_cell, v_lo, v_hi, v_step)

        bin_list[particle.octant].w += particle.w

        bin_list[particle.octant].v += particle.w * particle.v
    end

    for bin in bin_list
        if (bin.w > 0.0)
            bin.v = bin.v / bin.w
        end
    end

    for particle in particle_list
        bin_list[particle.octant].c += particle.w * ((particle.v .- bin_list[particle.octant].v).^2)
    end

    for bin in bin_list
        if (bin.w > 0.0)
             bin.c = sqrt.(bin.c / bin.w)
        end
    end

    return bin_list
end

function create_new_particles(bin_list)
    # create new particles from the merging quantities
    new_particles = []

    for bin in bin_list
        if (bin.w > 0)
            v_sign = sign.(0.5 .- rand(3))

            push!(new_particles, Particle(bin.w / 2, (bin.v .+ (v_sign .* bin.c)), 1))
            push!(new_particles, Particle(bin.w / 2, (bin.v .- (v_sign .* bin.c)), 1))
        end
    end

    return new_particles
end

function compute_stats(particle_list, moment_nos)
    # compute stats so that we can check mass, momentum and energy conservation
    ndens = 0.0
    vel = [0.0, 0.0, 0.0]

    for particle in particle_list
        ndens += particle.w
        vel += particle.w .* particle.v
    end

    vel = vel ./ ndens

    print("number density: ", ndens, "\n")
    print("velocity: ", vel, "\n")

    moments = fill(0.0, (length(moment_nos), ))

    energy = 0.0

    for particle in particle_list

        energy += particle.w * sum((particle.v .- vel) .^ 2)
    end
    println("Energy: ", energy / ndens)

    for particle in particle_list

        for (i, moment_no) in enumerate(moment_nos)
            moments[i] += particle.w * sum((particle.v .- vel) .^ moment_no)
        end
    end

    for (i, moment_no) in enumerate(moment_nos)
        print("moment no.", moment_no, ": ", moments[i] / ndens, "\n")
    end
end

function create_particles(mass, T, n_particles)
    particle_list = Array{Particle}(undef, n_particles)

    for i=1:n_particles
        particle_list[i] = Particle(1.0, maxwellian(T, mass), 0)
    end
    return particle_list
end

function main()
    ar_mass = 6.633521356992e-26;
    T = 1000.0  # temperature of gas
    n_particles = 10000 # number of particles to be sampled

    v_ref = sqrt(2 * constants.k * T / ar_mass)

    moment_nos = [4,6,8,10]

    # octree merging settings
    target_np = 128 # target number of particles
    aggr = 1.0
    max_refinement = 5  # maximum number of times a bin can be split

    inittype = MinMaxVel
    # inittype = LightSpeed

    splittype = MiddleSplit
    # splittype = MeanVelSplit

    refinement = Greedy
    # refinement = Heaviest

    # cell-based merging settings
    # number of cells in each x/y/z velocity direction
    N_cell = 4

    # extent of merging region ([-max_vel_to_v_ref * v_ref, max_vel_to_v_ref * v_ref])
    # scaled w.r.t. v_ref
    max_vel_to_v_ref = 4.5

    particle_list = create_particles(ar_mass, T, n_particles);
    println("Original particles")
    compute_stats(particle_list, moment_nos)
    println()

    println("Octree merge")
    binlist_test = [init_bin(particle_list, inittype, splittype)];
    @time octree_merge!(particle_list, binlist_test, target_np, 5, aggr, splittype, refinement);
    newparticles = create_new_particles(binlist_test);

    compute_stats(newparticles, moment_nos)
    println()

    println("Cell-based merge")
    @time binlist_test = cell_merge!(particle_list, N_cell,
                                     [-v_ref * max_vel_to_v_ref, -v_ref * max_vel_to_v_ref, -v_ref * max_vel_to_v_ref],
                                     [v_ref * max_vel_to_v_ref, v_ref * max_vel_to_v_ref, v_ref * max_vel_to_v_ref]);
    newparticles = create_new_particles(binlist_test);

    compute_stats(newparticles, moment_nos)
end

main()

max_stars_repo_name: knstmrd/particlemerging

max_issues_repo_name: knstmrd/particlemerging

max_forks_repo_name: knstmrd/particlemerging


text:
using DistributionFits
using Documenter
using StatsPlots, Statistics # used in help

# use plot on headleass server:
# https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block
ENV["GKSwstype"] = "100"

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

makedocs(;
    modules=[DistributionFits],
    authors="Thomas Wutzler <twutz@bgc-jena.mpg.de> and contributors",
    repo="https://github.com/bgctw/DistributionFits.jl/blob/{commit}{path}#{line}",
    sitename="DistributionFits.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bgctw.github.io/DistributionFits.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Dependencies" => "set_optimize.md",
        "Gamma" => "gamma.md",
        "LogNormal" => "lognormal.md",
        "LogitNormal" => "logitnormal.md",
        "Weibull" => "weibull.md",
        #"Details" => "z_autodocs.md",
    ],
)

deploydocs(;
    repo="github.com/bgctw/DistributionFits.jl",
    devbranch="main",
)

ext: jl


text:
import GeoInterface
using Test, LibGEOS

version = LibGEOS.GEOSversion()
@info "GEOS version $version"

if version != LibGEOS.GEOS_CAPI_VERSION
    @warn string(
        "The GEOS version in use is different from the version used to ",
        "generate libgeos_api.jl. Consider regenerating if the difference is more than a ",
        "patch release.",
    ) run_version = version gen_version = LibGEOS.GEOS_CAPI_VERSION
end

@testset "LibGEOS" begin
    include("test_geos_functions.jl")
    include("test_geos_operations.jl")
    include("test_geo_interface.jl")
    include("test_regressions.jl")
    include("test_invalid_geometry.jl")
    include("test_strtree.jl")
end


text:
# Hermit polynomial related functions for mode coupling analysis
# 2014 - murozturk@ku.edu.tr

# Returns values of hermite polynomials for x
# 0th hermite polynomial is always 1 so it is not passed around
# degree number of values are returned as an array

function hermit (x::Float64, degree::Int=32)
	Hx = Array(Float64, degree+1)
	Hx[1:2] = [1. x]
	for i=3:degree+1
		Hx[i] = x*Hx[i-1] - (i-2)*Hx[i-2]
	end
	return Hx[2:end]
end
# vectorized version
function hermit (x::Vector{Float64} )
	degree = length(hermit(0.))
	Hx = Array(Float64, length(x), degree)
	for i=1:length(x)
		Hx[i,:]=hermit(x[i])
	end
	return Hx
end

# Returns probabilites based on hermit polynomials

function probHermit (x::Float64, coef::Vector{Float64})
	#exp(-(x*x)/2)/sqrt(2*pi)*(1+sum(hermit(x).*coef))
	# sum(a.*b)= a'*b
	# 1/sqrt(2pi)= 0.3989422804014327
	# hermit(x)[1:2] always zero
	exp(-(x*x)/2)*(1+ (hermit(x)[3:end]'*coef[3:end])[1] )*0.3989422804014327

end
# vectorized version
function probHermit (x::Vector{Float64}, coef::Vector{Float64})
	Float64[probHermit(r,coef) for r in x]
end

# Calculates hermit coefficients for a given distribution
# C_i = E[ H_i ] / i!

function coefHermit( r::Matrix{Float64} )

	println("Calculating Hermite coefficients..")

	nSamples, nModes = size(r)
	nBins=sqrt(nSamples)
	degree = length(hermit(0.))

	coef = Array(Float64, degree, nModes)

	for mode=1:nModes
		rHat, rDist = dist(r[:,mode])
		rHermit = hermit(rHat)
		for d=1:degree
			h = rHermit[:,d]
			# expectation of dth Hermite for this mode
			Ehr = sum(h.*rDist)
			# divided by (d)!
			coef[d, mode] = divideByFactorial(Ehr, d)
		end
	end

	return  coef
end

# Returns samples from a Hermite distribution

function sampleHermit(nSamples::Int, coef::Vector{Float64})
	sample::Float64 = randn()
	#samples = randn(nSamples)

	samples = Array(Float64, nSamples)
	dice::Float64   = rand()
	for i=1:nSamples
		println(string(i,"/",nSamples))
		sample = randn()
		dice   = rand()
		while dice > probHermit(sample, coef)
			sample  = randn()
			dice    = rand()
		end
		samples[i]=sample
	end
	return samples
end

function sampleHermit(nSamples::Int, coef::Matrix{Float64})
	nModes = size(coef,2)

	sample::Float64 = randn()
	samples = Array(Float64, nSamples, nModes)

	println ("Sampling f1 distribution...")

	dice::Float64 = 1.

	for mode=1:nModes
		println(string(mode,"/",nModes));
		for i=1:nSamples
			dice   = 1.0
			while dice > probHermit(sample, coef[:,mode])
				#sample  = randn() #gaussian kernel
				sample = rand()*6-3 #uniform kernel
				dice   = rand()
			end
			samples[i,mode] = sample
		end
	end

	return samples
end

# Returns x/f!

function divideByFactorial (x::Float64, factorial::Int)
	factorial==0 && return x
	factorial==1 && return x
	divideByFactorial(x/factorial, factorial-1)
end

# returns prob. distributions of a given sample

function dist (x::Vector{Float64})
	nSamples = length(x)
	nBins=int(sqrt(nSamples))

	bins= linspace(minimum(x), maximum(x), nBins+1)

	px= hist(x,bins)[2]/nSamples

	index = (bins[1:nBins] + bins[2:nBins+1] ) /2

	return index, px
end

# this results in a common bin for all modes
# but actual n nBins is not well defined
# only good for quick comparisons

function dist (X::Matrix{Float64})
	nSamples, nModes = size(X)
	nBins=int(sqrt(nSamples))

	bins=linspace(minimum(X),maximum(X),nBins+1)

	pX = hist(X,bins)[2] /nSamples

	rHat  = ( bins[1:nBins] + bins[2:nBins+1] ) / 2

	return rHat , pX
end

# Calculates the expectation values for a given pair of modes

function f2_exp ( m1::Vector{Float64}, m2::Vector{Float64}, ind1::Int64 , ind2::Int64 , degree::Int = 32)

        bin_number=int(ceil(sqrt(length(m1))));

        combined = hcat(m1,m2);

        b1, b2, bins = hist2d(combined, linspace(minimum(m1),maximum(m1),bin_number), linspace(minimum(m2), maximum(m2),bin_number))

        # Normalization for a fixed bin base area

        base_area = (b1[2]-b1[1]) * (b2[2]-b2[1]);

        area = 0;

        for i=1:size(bins,1);
                for j=1:size(bins,2)
                        area += bins[i,j];
                end
        end

        bins /= (area);

        HxHy=zeros(495,5);

	# To keep track of the mode indices

        HxHy[:,3] = max(ind1,ind2);

        HxHy[:,4] = min(ind1,ind2);

        # Main loop

        for s1=1:bin_number-1

                for s2=1:bin_number-1

                        vec1 = modeCouplingAnalysis.hermit(b1[s1]);
                        vec2 = modeCouplingAnalysis.hermit(b2[s2]);

                        n = 1;

                        for v=3:degree

                                for p=1:(v-1)

                                        HxHy[n,5] += bins[s1,s2] * vec1[v-p] * vec2[p]
                                        HxHy[n,1] = p;

                                        HxHy[n,2] = v-p;

                                        n += 1;
                                end

                        end

                end

        end

        return HxHy

end

function combination (x1::Int64, x2::Int64)

	# Calculates and returns (x1,x2)
	if x1 >= x2
	return factorial(x1)/factorial(x2)/factorial(x1-x2)
	end
end

function coef_f2( f2_exp :: Array )

        # Calculates the f2 coefficients as in Eq. 8

        coef_f2 = f2_exp;

        for i=1:length(f2_exp)

                coef_f2[n,5] = combination(f2_exp[n,2] + f2_exp[n,1], f2_exp[n.1]) * f2_exp[n,5] ;

        end

        return coef_f2;

end


text:
using Pkg.TOML
using Test

cd(joinpath(@__DIR__, ".."))

contents = readdir()
for c in contents
    if isdir(c)
        print("Entering $c\n")
        cd(c)
        files = readdir(".")
        for f in files
            if splitext(f)[2] == ".toml"
                print("parsing $f\n")
                @test_nowarn TOML.parsefile(f)
            end
        end
        cd("..")
    end
end

max_stars_repo_name: JuliaPOMDP/Registry

max_issues_repo_name: JuliaPOMDP/Registry

max_forks_repo_name: JuliaPOMDP/Registry


text:
#ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля без внутренних перегородок и маркеров.
#РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.
function crest!(r::Robot)
    for side in (HorizonSide(i) for i=0:3)
        putmarkers!(r,side)
        go_by_markers(r,reverse(side))
    end
    putmarker!(r)
end
function putmarkers!(r::Robot,side::HorizonSide)
while isborder(r,side)==false
    move!(r,side)
    putmarker!(r)
end
end
function go_by_markers(r::Robot,side::HorizonSide)
while ismarker(r)==true
    move!(r,side)
end
end
function reverse(side::HorizonSide)
    HorizonSide(mod(Int(side)+2, 4))
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
"""
    RLS(y, x, λ)

Recursive Least Squares algorithm.

# Arguments
- `y::REALVEC`
- `x::REALMAT`
- `λ::Real = 1`

# Examples
```julia
w = [1, -1]
σ = 2
x = rand(100, 2)
y = x * w + randn(100) * σ
model = RLS(y, x, 0.99)
```
"""
function RLS(y::REALVEC,
             x::REALMAT,
             λ::Real = 0.99)

    @assert λ <= 1 && λ > 0 "Forgetting factor, λ, must be ]0;1]."

    T, D = size(x)

    w = zeros(T + 1, D)
    P = zeros(T + 1, D, D)
    μ = zeros(T)
    σ² = zeros(T)

    P[1, :, :] = Matrix(1000.0 * I, D, D)
    for t in 1:T
        # Predict
        μ[t]  = predict(x[t, :], w[t, :])

        # Update
        w[t + 1, :], P[t + 1, :, :] = update(y[t], x[t, :], w[t, :], P[t, :, :], μ[t], λ)
    end
    return (predictions = μ, coefficients = w, covariance = P)
end

function RLS(y::REALVEC,
             x::REALVEC,
             λ::Real = 0.99)
    RLS(y, reshape(x, :, 1), λ)
end

function predict(x, w)
    return dot(x, w)
end

function update(y, x, w, P, μ, λ)
    Q = λ + x' * P * x
    K = P * x / Q
    w = w + (y - μ) * K
    P = (I - K * x') * P / λ
    return (w, P)
end

predict(x::REALMAT, w::Vector{<:AbstractFloat}) = x * w
predict(x::REALVEC, w::Vector{<:AbstractFloat}) = dot(x, w)


text:
"""
    shoot(ϕ,xspan,g₁,g₂,init)

Shooting method to solve a two-point boundary value problem with
ODE u'' = `ϕ`(x,u,u') for x in `xspan`, left boundary condition
`g₁`(u,u')=0, and right boundary condition `g₂`(u,u')=0. The
value `init` is an initial estimate for vector [u,u'] at x=a.

Returns vectors for the nodes, the solution u, and derivative u'.
"""
function shoot(ϕ,xspan,g₁,g₂,init,tol=1e-5)
    # ODE posed as a first-order equation in 2 variables.
    shootivp = (v,p,x) -> [ v[2]; ϕ(x,v[1],v[2]) ]

    # Evaluate the difference between computed and target values at x=b.
    function objective(s)
        IVP = ODEProblem(shootivp,s,float.(xspan))
        sol = solve(IVP,Tsit5(),abstol=tol/10,reltol=tol/10)
        x = sol.t;  y = sol;
        return [g₁(s...),g₂(y[end]...)]
    end

    # Find the unknown quantity at x=a by rootfinding.
    x = [];  y = [];   # these values will be overwritten
    s = levenberg(objective,init,xtol=tol)[:,end]

    # Use the stored last solution of the IVP.
    u,du_dx = eachrow(y)
    return x,u,du_dx
end

"""
    diffmat2(n,xspan)

Compute 2nd-order-accurate differentiation matrices on `n`+1 points
in the interval `xspan`. Returns a vector of nodes and the matrices
for the first and second derivatives.
"""
function diffmat2(n,xspan)
    a,b = xspan
    h = (b-a)/n
    x = [ a + i*h for i in 0:n ]   # nodes

    # Define most of Dₓ by its diagonals.
    dp = fill(0.5/h,n)        # superdiagonal
    dm = fill(-0.5/h,n)       # subdiagonal
    Dₓ = diagm(-1=>dm,1=>dp)

    # Fix first and last rows.
    Dₓ[1,1:3] = [-1.5,2,-0.5]/h
    Dₓ[n+1,n-1:n+1] = [0.5,-2,1.5]/h

    # Define most of Dₓₓ by its diagonals.
    d0 =  fill(-2/h^2,n+1)    # main diagonal
    dp =  ones(n)/h^2         # super- and subdiagonal
    Dₓₓ = diagm(-1=>dp,0=>d0,1=>dp)

    # Fix first and last rows.
    Dₓₓ[1,1:4] = [2,-5,4,-1]/h^2
    Dₓₓ[n+1,n-2:n+1] = [-1,4,-5,2]/h^2

    return x,Dₓ,Dₓₓ
end

"""
    diffcheb(n,xspan)

Compute Chebyshev differentiation matrices on `n`+1 points in the
interval `xspan`. Returns a vector of nodes and the matrices for the
first and second derivatives.
"""
function diffcheb(n,xspan)
    x = [ -cos( k*π/n ) for k in 0:n ]    # nodes in [-1,1]

    # Off-diagonal entries.
    c = [2; ones(n-1); 2];    # endpoint factors
    dij = (i,j) -> (-1)^(i+j)*c[i+1]/(c[j+1]*(x[i+1]-x[j+1]))
    Dₓ = [ dij(i,j) for i in 0:n, j in 0:n ]

    # Diagonal entries.
    Dₓ[isinf.(Dₓ)] .= 0         # fix divisions by zero on diagonal
    s = sum(Dₓ,dims=2)
    Dₓ -= diagm(s[:,1])         # "negative sum trick"

    # Transplant to [a,b].
    a,b = xspan
    x = @. a + (b-a)*(x+1)/2
    Dₓ = 2*Dₓ/(b-a)             # chain rule

    # Second derivative.
    Dₓₓ = Dₓ^2
    return x,Dₓ,Dₓₓ
end

"""
    bvplin(p,q,r,xspan,lval,rval,n)

Use finite differences to solve a linear bopundary value problem.
The ODE is u''+`p`(x)u'+`q`(x)u = `r`(x) on the interval `xspan`,
with endpoint function values given as `lval` and `rval`. There will
be `n`+1 equally spaced nodes, including the endpoints.

Returns vectors of the nodes and the solution values.
"""
function bvplin(p,q,r,xspan,lval,rval,n)
    x,Dₓ,Dₓₓ = diffmat2(n,xspan)

    P = diagm(p.(x))
    Q = diagm(q.(x))
    L = Dₓₓ + P*Dₓ + Q     # ODE expressed at the nodes

    # Replace first and last rows using boundary conditions.
    z = zeros(1,n)
    A = [ [1 z]; L[2:n,:]; [z 1] ]
    b = [ lval; r.(x[2:n]); rval ]

    # Solve the system.
    u = A\b
    return x,u
end

"""
    bvp(ϕ,xspan,lval,lder,rval,rder,init)

Finite differences to solve a two-point boundary value problem with
ODE u'' = `ϕ`(x,u,u') for x in `xspan`, left boundary condition
`g₁`(u,u')=0, and right boundary condition `g₂`(u,u')=0. The value
`init` is an initial estimate for the values of the solution u at
equally spaced values of x, which also sets the number of nodes.

Returns vectors for the nodes and the values of u.
"""
function bvp(ϕ,xspan,g₁,g₂,init)
    n = length(init) - 1
    x,Dₓ,Dₓₓ = diffmat2(n,xspan)
    h = x[2]-x[1]

    function residual(u)
        # Residual of the ODE at the nodes.
        du_dx = Dₓ*u                   # discrete u'
        d2u_dx2 = Dₓₓ*u                # discrete u''
        f = d2u_dx2 - ϕ.(x,u,du_dx)

        # Replace first and last values by boundary conditions.
        f[1] = g₁(u[1],du_dx[1])/h
        f[n+1] = g₂(u[n+1],du_dx[n+1])/h
        return f
    end

    u = levenberg(residual,init)
    return x,u[end]
end

"""
    fem(c,s,f,a,b,n)

Use a piecewise linear finite element method to solve a two-point
boundary value problem. The ODE is (`c`(x)u')' + `s`(x)u = `f`(x) on
the interval [`a`,`b`], and the boundary values are zero. The
discretization uses `n` equal subintervals.

Return vectors for the nodes and the values of u.
"""
function fem(c,s,f,a,b,n)
    # Define the grid.
    h = (b-a)/n
    x = @. a + h*(0:n)

    # Templates for the subinterval matrix and vector contributions.
    Ke = [1 -1; -1 1]
    Me = (1/6)*[2 1; 1 2]
    fe = (1/2)*[1; 1]

    # Evaluate coefficent functions and find average values.
    cval = c.(x);   cbar = (cval[1:n]+cval[2:n+1]) / 2;
    sval = s.(x);   sbar = (sval[1:n]+sval[2:n+1]) / 2;
    fval = f.(x);   fbar = (fval[1:n]+fval[2:n+1]) / 2;

    # Assemble global system, one interval at a time.
    K = zeros(n-1,n-1);  M = zeros(n-1,n-1);  f = zeros(n-1);
    K[1,1] = cbar[1]/h;  M[1,1] = sbar[1]*h/3;  f[1] = fbar[1]*h/2;
    K[n-1,n-1] = cbar[n]/h;  M[n-1,n-1] = sbar[n]*h/3;  f[n-1] = fbar[n]*h/2;
    for k in 2:n-1
        K[k-1:k,k-1:k] += (cbar[k]/h) * Ke
        M[k-1:k,k-1:k] += (sbar[k]*h) * Me
        f[k-1:k] += (fbar[k]*h) * fe
    end

    # Solve system for the interior values.
    u = (K+M) \ f
    u = [0; u; 0]      # put the boundary values into the result
    return x,u
end


text:
struct ScaledPSDCone <: MOI.AbstractVectorSet
    side_dimension::Int
end

Base.copy(x::ScaledPSDCone) = ScaledPSDCone(x.side_dimension)

MOI.side_dimension(x::ScaledPSDCone) = x.side_dimension

function MOI.dimension(x::ScaledPSDCone)
    return x.side_dimension^2
end

function MOI.Utilities.set_with_dimension(::Type{ScaledPSDCone}, dim)
    return ScaledPSDCone(isqrt(dim))
end

struct ScaledPSDConeBridge{T,G} <: MOI.Bridges.Constraint.SetMapBridge{
    T,
    ScaledPSDCone,
    MOI.PositiveSemidefiniteConeTriangle,
    MOI.VectorAffineFunction{T},
    G,
}
    constraint::MOI.ConstraintIndex{MOI.VectorAffineFunction{T},ScaledPSDCone}
end

function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{ScaledPSDConeBridge{T}},
    ::Type{G},
    ::Type{MOI.PositiveSemidefiniteConeTriangle},
) where {T,G<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction{T}}}
    return ScaledPSDConeBridge{T,G}
end

function MOI.Bridges.map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::MOI.PositiveSemidefiniteConeTriangle,
)
    return ScaledPSDCone(set.side_dimension)
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::ScaledPSDCone,
)
    return MOI.PositiveSemidefiniteConeTriangle(set.side_dimension)
end

# Map ConstraintFunction from MOI -> SCS
function MOI.Bridges.map_function(BT::Type{<:ScaledPSDConeBridge{T}}, func::MOI.VectorOfVariables) where {T}
    new_f = MOI.Utilities.operate(*, Float64, 1.0, func)
    return MOI.Bridges.map_function(BT, new_f)
end
function MOI.Bridges.map_function(::Type{<:ScaledPSDConeBridge}, f::MOI.VectorAffineFunction)
    n = MOI.output_dimension(f)
    d = MOI.Utilities.side_dimension_for_vectorized_dimension(n)
    constants = copy(f.constants)
    scale_coefficients!(constants)
    constants = triangle_to_square(constants, d)
    terms = copy(f.terms)
    scale_coefficients!(terms)
    triangle_to_square_indices!(terms, d)
    return MOI.VectorAffineFunction(terms, constants)
end

# Used to map the ConstraintPrimal from SCS -> MOI
function MOI.Bridges.inverse_map_function(::Type{<:ScaledPSDConeBridge}, square)
    triangle = square_to_triangle(square)
    unscale_coefficients!(triangle)
    return triangle
end

# Used to map the ConstraintDual from SCS -> MOI
function MOI.Bridges.adjoint_map_function(::Type{<:ScaledPSDConeBridge}, square)
    triangle = square_to_triangle(square)
    n = isqrt(length(square))
    for i in 1:n, j in 1:(i-1)
        # Add lower diagonal dual. It should be equal to upper diagonal dual
        # but `unscale_coefficients` will divide by 2 so it will do the mean
        triangle[triangle_map(i, j)] += square[square_map(i, j, n)]
    end
    unscale_coefficients!(triangle)
    return triangle
end

# `dimension` -> `side_dimension`, see
# http://jump.dev/MathOptInterface.jl/v0.8.1/apireference/#MathOptInterface.PositiveSemidefiniteConeTriangle
triangle_side_dimension(n) = div(isqrt(1 + 8n) - 1, 2)
square_side_dimension(n) = isqrt(n)

# Matrix indices -> Index in vectorized form
function triangle_map(i::Integer, j::Integer)
    if i < j
        return triangle_map(j, i)
    else
        # See http://jump.dev/MathOptInterface.jl/v0.8.1/apireference/#MathOptInterface.PositiveSemidefiniteConeTriangle
        return div((i - 1) * i, 2) + j
    end
end
function square_map(i::Integer, j::Integer, n::Integer)
    if i < j
        return square_map(j, i, n)
    else
        return i + (j - 1) * n
    end
end

function copy_upper_triangle(x, n, map_from, map_to)
    y = zeros(eltype(x), map_to(n, n))
    for i in 1:n, j in 1:i
        y[map_to(i, j)] = x[map_from(i, j)]
    end
    return y
end
function square_to_triangle(x, n=square_side_dimension(length(x)))
    return copy_upper_triangle(x, n, (i, j) -> square_map(i, j, n),
                               triangle_map)
end
function triangle_to_square(x, n=triangle_side_dimension(length(x)))
    return copy_upper_triangle(x, n, triangle_map, (i, j) -> square_map(i, j, n))
end

function triangle_to_square_indices!(
    x::Vector{<:MOI.VectorAffineTerm},
    n,
)
    map = square_to_triangle(1:n^2, n)
    for i in eachindex(x)
        x[i] = MOI.VectorAffineTerm(
            map[x[i].output_index],
            x[i].scalar_term,
        )
    end
end

_row(i, t::MOI.VectorAffineTerm) = t.output_index
_row(i, β) = i

_prod(α, t::MOI.VectorAffineTerm) = MOI.Utilities.operate_term(*, α, t)
_prod(α, β) = α * β

# Scale coefficients depending on rows index on symmetric packed upper triangular form
# coef: List of coefficients
# rev: if true, we unscale instead (e.g. divide by √2 instead of multiply for PSD cone)
# rows: List of row indices
# d: dimension of set
function _scale_coefficients!(coef::AbstractVector, rev::Bool)
    scaling = rev ? 0.5 : 2.0
    for i in 1:length(coef)
        if !MOI.Utilities.is_diagonal_vectorized_index(_row(i, coef[i]))
            coef[i] = _prod(scaling, coef[i])
        end
    end
end
# Scale the coefficients in `coef` with respective rows in `rows` for a set `s`
function scale_coefficients!(coef)
    return _scale_coefficients!(coef, false)
end
# Unscale the coefficients of `coef` in symmetric packed upper triangular form
function unscale_coefficients!(coef)
    return _scale_coefficients!(coef, true)
end


text:
#=
"""
--------------------------------------------------------------------------------
Function: XP_Bisection
--------------------------------------------------------------------------------
Description: Bisects box in either X or P dimension.
--------------------------------------------------------------------------------
Inputs:
h           function - h(z,p)
hj          function - Jacobian of h(z,p) w.r.t. z
CurrNode    Vector{Interval{Float64}} - Current X node
P           Vector{Interval{Float64}} - Current P node
NextNode    Vector{Interval{Float64}} - New X node
NextNodeP   Vector{Interval{Float64}} - New P node 1
Ptemp       Vector{Interval{Float64}} - New P node 2
PIflag      Bool - Partial Inclusion Flag
Iflag       Bool - Inclusion flag
Eflag       Bool - Exclusion flag
ptol        Float64 - tolerance for P dimension size
pbisect     Bool - bisect based on relative width?
Pstart      Vector{Interval{Float64}} - initial P size
--------------------------------------------------------------------------------
Returns:
bp          Int64 - Cut type flag
CurrNode    Vector{Interval{Float64}} - Current X node
P           Vector{Interval{Float64}} - Current P node
NextNode    Vector{Interval{Float64}} - New X node
NextNodeP   Vector{Interval{Float64}} - New P node 1
Ptemp       Vector{Interval{Float64}} - New P node 2
--------------------------------------------------------------------------------
"""

"""
--------------------------------------------------------------------------------
Function: bisectDirX
--------------------------------------------------------------------------------
Description: Determine direction and position of cut in X dimension.
--------------------------------------------------------------------------------
Inputs:
h        function - h(z,p)
X        Vector{Interval{Float64}} - for state variables
P        Vector{Interval{Float64}} - for decision variables
x        Currently Unused
PIFlag   Bool - Flag for partial inclusion.
--------------------------------------------------------------------------------
Returns:
jmax            Int64 - The dimension to cut X in
xCut            Float64 - Cut value for X
piflagtemp      Bool - Partial Inclusion Flag
--------------------------------------------------------------------------------
"""

"""
--------------------------------------------------------------------------------
Function: bisectDirP
--------------------------------------------------------------------------------
Description: Determines the direction to branch P in based on absolute/relative
width and tolerance for P dimension size.
--------------------------------------------------------------------------------
Inputs:
P         Vector{Interval{Float64}} - P Box to bisect
ptol      Float64 - tolerance for P dimension size
pbisect   Bool - bisect based on relative width?
Pstart    Vector{Interval{Float64}} - initial P size
--------------------------------------------------------------------------------
Returns:
jmax - The dimension to branch P in.
--------------------------------------------------------------------------------
"""

"""
--------------------------------------------------------------------------------
Function:  partIncProc
--------------------------------------------------------------------------------
Description:
Partial inclusion check (uses Newton type operator test with Miranda)
--------------------------------------------------------------------------------
Inputs:
h         function - h(z,p)
hj        function - Jacobian of h(z,p) w.r.t. z
X0        Vector{Interval{Float64}} - state variables
P         Vector{Interval{Float64}} - decision variables
exDim     Int64 - Dimension flag
Iflag     Bool - Inclusion flag
Eflag     Bool - Exclusion flag
PIflag    Bool - Partial Inclusion Flag
--------------------------------------------------------------------------------
Returns:
X         IntervalBox - State variables
Iflag     Bool - Inclusion flag
Eflag     Bool - Exclusion flag
PIflag    Bool - Partial Inclusion Flag
--------------------------------------------------------------------------------
"""
=#

#=
function XP_Bisection(h,hj,CurrNode,P,NextNode,NextNodeP,Ptemp,PIflag,Iflag,Eflag,ptol,pbisect,Pstart)
  #println("CurrNode Top: ", CurrNode)
  pPart = false
  xCut = 0.15
  pCutL = 0.0
  pCutR = 0.0
  IflagL = copy(Iflag)
  IflagR = copy(Iflag)
  EflagL = copy(Eflag)
  EflagR = copy(Eflag)
  XLeft = copy(CurrNode)
  XRight = copy(NextNode)
  PLeft = copy(P)
  PRight = copy(NextNodeP)
  Pmid = copy(Ptemp)
  j = 0
  jp = 0
  bp = 0
  j, xCut,PIflag = bisectDirX(h,CurrNode,P,xCut,PIflag)
  #println("bisectDirX j: ", j)
  #println("bisectDirX xCut: ", xCut)
  #println("bisectDirX PIflag: ", PIflag)
  PIflagL = PIflag
  PIflagR = PIflag
  #println("PIflag:  ", PIflag)
  if (PIflag == false)
    XLeft[j] = Interval(CurrNode[j].lo,xCut)
    XRight[j] = Interval(xCut,CurrNode[j].hi)
    CurrNode = copy(XLeft)
    NextNode = copy(XRight)
    println("bisect trace 1")
    return 1,CurrNode,P,NextNode,NextNodeP,Ptemp
  else
    jp = bisectDirP(P,ptol,pbisect,Pstart)
    if (jp != -1)
      if pPart==true
        bp,jp,pCutL,pCutR = branchP(CurrNode,P,jp,pCutL,pCutR)
      else
        bp = 1
      end
      #println("jp")
      if (bp == 1)
        PLeft[jp] = Interval(P[jp].lo,mid(P[jp]))
        PRight[jp] = Interval(mid(P[jp]),NextNodeP[jp].hi)
        P = copy(PLeft)
        NextNodeP = copy(PRight)
        #println("bisect trace 2")
        return 2,CurrNode,P,NextNode,NextNodeP,Ptemp
      elseif (bp == 3)
        PLeft[jp] = Interval(P[jp].lo,pCutL)
        Pmid[jp] = Interval(pCutL,pCutR)
        PRight[jp] = Interval(pCutR,P[jp].hi)
        P = copy(PLeft)
        Ptemp = copy(Pmid)
        NextNodeP = copy(PRight)
        #println("bisect trace 3")
        return 3,CurrNode,P,NextNode,NextNodeP,Ptemp
      end
    end
  end
  #println("jp = :",jp)
  #println("bp = :",bp)
  if (false) # originally just ((jp == -1) && (bp != 1))
    xCut = xCut*CurrNode[j].lo + (1.0-xCut)*CurreNode[j].hi
    XLeft[j] = Interval(CurrNode[j].lo,xCut)
    XRight[j] = Interval(xCut,CurrNode[j].hi)
    k,rhoL = spectralR(hj,XLeft,P,rhoL)
    k,rhoR = spectralR(hj,XRight,P,rhoR)
    if (rhoL)
      XLeft,Iflag,Eflag,PIflag = partIncProc(h,hj,XLeft,P,j,Iflag,Eflag,PIflag)
    end
    if (rhoR)
      XRight,Iflag,Eflag,PIflag = partIncProc(h,hj,XRight,P,j,Iflag,Eflag,PIflag)
    end
    if ((~PIflagL) && (~PIflagR))
      PIflag = false
    end
  end

  CurrNode = copy(XLeft)
  NextNode = copy(XRight)
  #println("bisect trace 4")
  return bp,CurrNode,P,NextNode,NextNodeP,Ptemp
end

function bisectDirX(h,X,P,x,PIFlag)

  piflagtemp = copy(PIFlag)
  v = 1
  iterm = 4
  jmax = 1
  k = -1
  maxL = 0.0
  disc = 20
  Xold = copy(X)
  Xnew = copy(X)
  xCut = 0.0
  maxL,jmax = findmax(diam.(X))
  if (v == 1)
    rad = diam(X[jmax])/2.0
    for l=0:(disc/2)
      xCut = X[jmax].lo+rad/(iterm) + l*(diam(X[jmax])-2.0*rad/(iterm))/(disc-1)
      Xnew[jmax] = xCut
      k = Miranda(h,Xnew,P)
      if (k>=0)
        return jmax, xCut, piflagtemp
      end
      xCut = X[jmax].hi-rad/(iterm) - l*(diam(X[jmax])-2.0*rad/(iterm))/(disc-1)
      Xnew[jmax] = xCut
      k = Miranda(h,Xnew,P)
      if (k>=0)
        return jmax, xCut, piflagtemp
      else
        Xnew[jmax] = Xold[jmax]
      end
    end
  end

  if (k == -1)
    v = 2
  end

  if (v == 2)
    iterm = iterm*4
    Xnew = copy(Xold)
    disc = 4*disc
    for i = 1:length(X)
        rad = diam(X[i])/2.0
        for l = 0:(disc/2)
          xCut = X[i].lo+rad/(iterm) + l*(diam(X[i])-2.0*rad/(iterm))/(disc-1)
          Xnew[i] = xCut
          k = Miranda(h,Xnew,P)
          if (k>=0)
            return i, xCut, piflagtemp
          end
          xCut = X[i].hi-rad/(iterm) - l*(diam(X[i])-2.0*rad/(iterm))/(disc-1)
          Xnew[i] = xCut
          k = Miranda(h,Xnew,P)
          if (k>=0)
            return i, xCut, piflagtemp
          else
            Xnew[i] = Xold[i]
          end
        end
    end
    if (k == -1)
      piflagtemp = true
    end
  end
  return jmax, xCut, piflagtemp
end

function bisectDirP(P,ptol,pbisect,Pstart)
  jmax = 1
  maxL = 0.0
  if (~pbisect)
    for i=1:length(P)
      if (maxL < diam(P[i]))
        maxL = diam(P[i])
        jmax = i
      end
    end
    if (maxL < ptol)
      return -1
    end
  else
    for i=1:length(P)
      if (maxL < diam(P[i])/diam(Pstart[i]))
        maxL = diam(P[i])/diam(Pstart[i])
        jmax = i
      end
    end
    if (maxL < ptol)
      return -1
    end
  end
  return jmax
end

# Unfinish bisection based on spectral radius
function branchP(X,P,jmax,pCutL,pCutR)

  Plocal = P
  maxD = 0.0
  try
    Pf = [P[j] for j=1:length(P)]
    p = [abs(Pf[j].lo+Pf[j].hi)/2.0 for j=1:length(P)]
    Xf = [X[j] for j=1:length(X)]
    rhok = rho
    fpsecR(Xf,Pf,rho)
  catch
  end

end

function partIncProc(h,hj,X0,P,exDim,Iflag,Eflag,PIflag)

  nx = length(X)
  S1 = Interval(0.0)
  S2 = Interval(0.0)
  k = 0
  incl = -1
  side = 0
  inclusion = [false for i=1:length(X)]
  exclusion = false
  Iflag = false
  Eflag = false

  X = copy(X0)
  Xpart = copy(X0)
  Xpart[exDim] = xCut
  if (X[exDim].lo == xCut)
    side = 1
  end

  x_mid = mid.(X)
  H = h(x_mid,P)
  J = hj(X,P)
  Y = Preconditioner(hj,X,P,jac="User")
  B = Y*H
  M = Y*J
  Xi = X

  for i=1:nx
    S1 = S2 = Interval(0.0)
    for j=1:nx
      if (j<i)
        S1 += M[i,j]*(X[j]-x_mid[j])
      elseif (j>i)
        S2 += M[i,j]*(X[j]-x_mid[j])
      end
    end
    N[i] = x_mid[i] - (B[i]+S1+S2)/M[i,i]
    if N[i] ⊆ X[i]
      inclusion[i] = true
    else
      inclusion[i] = false
    end
    if ~isdisjoint(N[i],X[i])
      if (i == exDim)
        wL = X[i].lo - N[i].lo
        wU = N[i].hi - X[i].hi
        if (side == 0)
          X[i] = Interval(X[i].lo,X[i].hi+0.25*wU)
        else
          X[i] = Interval(X[i].lo-wL*0.25,X[i].hi)
        end
      else
        X[i] = N[i] ∩ X[i]
      end
    else
      exclusion = true
      break
    end
  end
  k += 1

  # checks if all components are included
  if (Iflag == false)
    for i=1:nx
      if (inclusion[i] == true)
        Iflag = true
        continue
      elseif (inclusion[i] == false)
        Iflag = false
        break
      end
    end
  end

  Xpart = X
  if (side == 0)
    Xpart[exDim] = X[exDim].hi
    if (Xpart[exDim]<=xCut)
      PIflag = false
      return X,Iflag,Eflag,PIflag #FIX ME
    end
  else
    Xpart[exDim] = X[exDim].lo
    if (Xpart[exDim]>=xCut)
      PIflag = false
      return X,Iflag,Eflag,PIflag #FIX ME
    end
  end
  incl = Miranda(h,Xpart,P)

  return X,Iflag,Eflag,PIflag #FIX ME
end
=#


text:
### A Pluto.jl notebook ###
# v0.16.3

using Markdown
using InteractiveUtils

# ╔═╡ 6446adac-a3e2-4961-bea5-620118149b50
begin
    import Pkg
    Pkg.add(PackageSpec(url="https://github.com/Pocket-titan/DarkMode"))
    import DarkMode
    DarkMode.enable()
    # OR DarkMode.Toolbox(theme="default")
end

# ╔═╡ 81bfcb8c-3345-11ec-21ec-bf37d08f624f
begin
	using DifferentialEquations
	using Plots
	using PlutoUI
	#using DarkMode; DarkMode.Toolbox(theme="default")
end;

# ╔═╡ 369a25d1-6c20-4b22-9f66-a845364b3a0d
function model!(du, u, p, t)
	ω0, ωF, ωM, ϕF, ϕM, μ, ν, αF, αM, ψF, ψM = p
	du[1] = ω0
	du[2] = ωF + μ*sin(u[3]-u[2]-ϕF) + αF*sin(u[1]-u[2]-ψF)
	du[3] = ωM + ν*sin(u[2]-u[3]-ϕM) + αM*sin(u[1]-u[3]-ψM)
end;

# ╔═╡ ecf11f46-b4fa-46c5-8c10-b833adcbd72b
begin
	θ  = rand(3).%(2π) # Initial phases (Sunlight, fibroblast, macrophage)
	tMax = 10.0 	   # Total time
	ω0 = 2π            # Sunlight phase rate of change /day
	ωF = 2π            # Fibroblast phase intrinsic rate of change /day
	ωM = 2π            # Macrophage phase intrinsic rate of change /day
	ϕF = 3π/5          # Fibroblast phase offset in coupling to macrophage
	ϕM = 2π/5 		   # Macrophage phase offset in coupling to fibroblast
	μ  = 1.0           # Amplitude of fibroblast phase coupling to macrophage
	ν  = 0.0           # Amplitude of macrophage phase coupling to fibroblast
	αF = 10.0          # Amplitude of fibroblast phase coupling to sunlight
	αM = 10.0          # Amplitude of macrophage phase coupling to sunlight
	ψF = π/8           # Fibroblast phase offset from sunlight
	ψM = π/4          # Macrophage phase offset from sunlight
	p = [ω0, ωF, ωM, ϕF, ϕM, μ, ν, αF, αM, ψF, ψM]

	prob = ODEProblem(model!,θ,(0.0,tMax),p)
	sol = solve(prob, Tsit5(),saveat=0.01)

	anim = @animate for u in sol.u
		scatter(sin.(u),cos.(u),xlims=(-1.1,1.1),ylims=(-1.1,1.1),aspect_ratio=:equal,marker=[:star8,:circle,:circle],color=[:red,:green,:blue],ms=10,series_annotations=["Light","Fib","Mac"],legend=:none)
	end
	gif(anim,"test.gif",fps=3)

end

# ╔═╡ Cell order:
# ╠═6446adac-a3e2-4961-bea5-620118149b50
# ╠═81bfcb8c-3345-11ec-21ec-bf37d08f624f
# ╠═369a25d1-6c20-4b22-9f66-a845364b3a0d
# ╠═ecf11f46-b4fa-46c5-8c10-b833adcbd72b

max_stars_repo_name: chris-revell/ClockOscillators

max_issues_repo_name: chris-revell/ClockOscillators

max_forks_repo_name: chris-revell/ClockOscillators


text:
@testset "redirect" begin
    gotd()
    write(joinpath(td, "config.md"), "")
    mkpath(joinpath(td, "_css"))
    mkpath(joinpath(td, "_layout"))
    write(joinpath(td, "_layout", "head.html"), "")
    write(joinpath(td, "_layout", "foot.html"), "")
    write(joinpath(td, "_layout", "page_foot.html"), "")
    write(joinpath(td, "foo.md"), """
    {{redirect bar/foo.html}}
    {{redirect bar.html}}
    # Hello
    baz
    """)
    write(joinpath(td, "index.md"), """
    # Index
    """)
    serve(single=true)
    @test isfile(joinpath("__site", "index.html"))
    @test isfile(joinpath("__site", "foo", "index.html"))
    @test isfile(joinpath("__site", "bar.html"))
    @test isfile(joinpath("__site", "bar", "foo.html"))
    @test read(joinpath("__site", "bar", "foo.html"), String) // """
                    <!-- Generated Redirect -->
                    <!doctype html>
                    <html>
                    <head>
                      <meta http-equiv="refresh" content="0; url=/foo/index.html">
                    </head>
                    </html>"""
end


text:
# fftconv.jl
# test FFT-based convolution and
# adjoint consistency for FFT convolution methods on very small case

using SPECTrecon: plan_psf
using SPECTrecon: fft_conv!, fft_conv_adj!, fft_conv_adj2!
using SPECTrecon: fft_conv, fft_conv_adj
using LinearMapsAA: LinearMapAA
using Test: @test, @testset, @test_throws, @inferred

@testset "plan_psf" begin
    plan = plan_psf(10, 10, 5)
    show(isinteractive() ? stdout : devnull, "text/plain", plan)
    show(isinteractive() ? stdout : devnull, "text/plain", plan[1])
    @test sizeof(plan) isa Int
    @test sizeof(plan[1]) isa Int
end

@testset "fftconv" begin
    img = randn(Float32, 12, 8)
    ker = rand(Float64, 7, 7)
    ker_sym = ker .+ reverse(ker, dims=:)
    ker_sym /= sum(ker_sym)
    out = @inferred fft_conv(img, ker_sym)
    @test eltype(out) == Float64
    @test_throws String fft_conv(img, ker)
end

@testset "fftconv3" begin
    nx = 12
    nz = 10
    nx_psf = 5
    T = Float32
    plan = plan_psf(nx, nz, nx_psf; T)
    image3 = rand(T, nx, nx, nz)
    ker3 = ones(T, nx_psf, nx_psf, nx) / (nx_psf)^2
    result = similar(image3)
    fft_conv!(result, image3, ker3, plan)
    @test maximum(result) ≤ 1
    fft_conv_adj!(result, image3, ker3, plan)
    @test maximum(result) ≤ 1.5 # boundary is the sum of replicate padding
    fft_conv_adj2!(result, image3[:, 3, :], ker3, plan)
    @test maximum(result) ≤ 1.5 # boundary is the sum of replicate padding
end

@testset "adjoint-fftconv" begin
    M = 20
    N = 14
    T = Float32
    for i = 1:4 # test with different kernels
        img = randn(T, M, N)
        ker = rand(T, 5, 5)
        ker = ker .+ reverse(ker, dims=:)
        ker /= sum(ker)

        idim = (M, N)
        odim = (M, N)
        forw = img -> fft_conv(img, ker)
        back = img -> fft_conv_adj(img, ker)

        A = LinearMapAA(forw, back, (prod(odim), prod(idim)); T, odim, idim)
        @test Matrix(A') ≈ Matrix(A)'
    end
end


text:
mutable struct VerticeType
    Point :: Array{Float64,1}
    adj ::  Array{Int64,1}
    J :: Array{Int64,1}
    T :: Float64
end


text:
#this file includes the creation of the various julian meshes, returning true if all of
#them are creates succesfully with no errors

using FEniCS
using PyCall
const fenics=pyimport("fenics")

test_triangle = UnitTriangleMesh()
#test_tetrahedron = UnitTetrahedronMesh() # deprecated? https://fenicsproject.org/docs/dolfin/2019.1.0/python/_autogenerated/dolfin.cpp.generation.html
test_interval = UnitIntervalMesh(10)
test_square = UnitSquareMesh(10,10,"crossed")
test_cube = UnitCubeMesh(10,10,10)
test_box = BoxMesh(Point([0.0,0.0,0.0]),Point([1.0,1.0,1.0]),10,10,10)
test_rectangle = RectangleMesh(Point([0.0, 0.0]), Point([10.0, 4.0]), 10, 10)
test_copy  = Mesh(test_square)
test_file = Mesh("./dolphin.xml")

#the below functions simply check the creation of the objects,
#without (currently) verifying values
x1 = cell_orientations(test_square)
x2 = cells(test_square)
x3 = hmin(test_square)
x4 = hmax(test_square)
x5 = coordinates(test_square)
x6 = data(test_square)
x7 = domains(test_square)
x8 = geometry(test_square)
x9 = num_cells(test_square)
x10 = num_edges(test_square)
x11 = num_entities(test_square,1)
x12 = num_faces(test_square)
x13 = num_facets(test_square)
x14 = num_vertices(test_square)
x15 = bounding_box_tree(test_square)
x16 = init(test_square)
x17 = init(test_square,1)
x18 = topology(test_square)
x19 = rmax(test_square)
x20 = rmin(test_square)
x21 = ufl_id(test_square)
x22 = ufl_domain(test_square)

try
  pyimport("mshr")
  circle1 = Circle(Point([0.0, 0.0]),3)
  circle2 = Circle(Point([0.0, 0.0]),4)
  rectangle = Rectangle(Point([0.0, 0.0]),Point([1.0, 1.0]))
  box = Box(Point([0.0, 0.0, 0.0]),Point([1.0, 1.0,1.0]))
  cone = Cone(Point([0.0,0.0,1.0]),Point([0.0, 0.0, 0.0]),3,50)
  sphere = Sphere(Point([0.0, 0.0, 0.0]),3)
  mesh1 = generate_mesh(circle1,64)
  mesh2= generate_mesh(circle2,64)
  circle3 = circle1+circle2
  circle4 = circle2-circle1
  circle5 = circle1*circle2
  mesh3 = generate_mesh(circle3,64)
  mesh4 = generate_mesh(circle4,64)
  union = CSGUnion(circle1,circle2)
  intersection = CSGIntersection(circle1,circle2)
  difference = CSGDifference(circle1,circle2)
  scaling = CSGScaling(circle1,2)
  rotation_1 = CSGRotation(rectangle,45)
  rotation_2 = CSGRotation(rectangle,Point([11.0, 111.0]),45)
  rotation_3 = CSGRotation(sphere,Point([1.0, 1.0]),Point([1.0, 2.0]),45)
  translation = CSGTranslation(rectangle,Point([2.0,2.0]))
  x = Extrude2D(union,0.2)
catch ee
 print("mshr has not been included, so it has not been tested")
end

list_lu_solver_methods()

list_krylov_solver_methods()

list_krylov_solver_preconditioners()

list_linear_solver_methods()

info_NonLinearVariationalSolver()

info_LinearVariatonalSolver()

true


text:
using Test
using NQCDynamics
using NQCDynamics.InitialConditions.QuantisedDiatomic
using LinearAlgebra: norm
using Unitful
using UnitfulAtomic

atoms = Atoms([:N, :O])
model = DiatomicHarmonic(r₀=2)
sim = Simulation(atoms, model)
surface = QuantisedDiatomic.SurfaceParameters(sim.atoms.masses, [1, 2], Matrix{Float64}(undef, 3, 0), 10.0, [0, 0, 1.0])

@testset "separate/combine_slab_and_molecule" begin
    atom_indices = [2, 3]
    r0 = [
        1 2 3 4 5;
        6 7 8 9 10;
        11 12 13 14 15
    ]
    molecule, slab = QuantisedDiatomic.separate_slab_and_molecule(atom_indices, r0)
    @test molecule == [2 3; 7 8; 12 13]
    @test slab == [1 4 5; 6 9 10; 11 14 15]

    r = QuantisedDiatomic.combine_slab_and_molecule(atom_indices, molecule, slab)
    @test r == r0

    atom_indices = [1, 2]
    r0 = [1 2; 3 4; 5 6]
    molecule, slab = QuantisedDiatomic.separate_slab_and_molecule(atom_indices, r0)
    @test molecule == r0

    r = QuantisedDiatomic.combine_slab_and_molecule(atom_indices, molecule, slab)
    @test r == r0
end

@testset "calculate_diatomic_energy" begin
    @test QuantisedDiatomic.calculate_diatomic_energy(model, 2.0, surface) ≈ 0.0 atol=1e-3
    @test QuantisedDiatomic.calculate_diatomic_energy(model, 3.0, surface) ≈ 0.5
end

@testset "calculate_force_constant" begin
    bond_lengths = 0.5:0.01:5.0
    potential = QuantisedDiatomic.calculate_diatomic_energy.(model, bond_lengths, surface) # Calculate binding curve

    @test QuantisedDiatomic.calculate_force_constant(bond_lengths, potential)[1] ≈ 1 atol=0.1

    model2 = NQCModels.DiatomicHarmonic(r₀=3.0)
    potential = QuantisedDiatomic.calculate_diatomic_energy.(model2, bond_lengths, surface) # Calculate binding curve
    @test QuantisedDiatomic.calculate_force_constant(bond_lengths, potential)[1] ≈ 1 atol=0.1
end

@testset "subtract_centre_of_mass!" begin
    atoms = Atoms([:O, :O])
    r = rand(3, 2)
    r = QuantisedDiatomic.subtract_centre_of_mass(r, atoms.masses)
    @test r[:,1] ≈ -r[:,2]
end

@testset "apply_random_rotation!" begin
    x = fill(5.0, 3, 2)
    y = fill(-4.0, 3, 2)
    @test all(x .== 5)
    @test all(y .== -4)
    QuantisedDiatomic.apply_random_rotation!(x, y)
    @test all(x .!= 5)
    @test all(y .!= -4)
    for col in eachcol(x)
        @test norm(col) ≈ norm([5, 5, 5])
    end
    for col in eachcol(y)
        @test norm(col) ≈ norm([-4, -4, -4])
    end
end

@testset "generate and check results" begin
    numbers = [2, 5]
    configs = generate_configurations(sim, numbers[1], numbers[2]; samples=10, translational_energy=1u"eV")
    for config in configs
        ν, J = quantise_diatomic(sim, config...)
        @test ν ≈ numbers[1] rtol=1e-1
        @test J ≈ numbers[2] rtol=1e-1
    end
end

@testset "position_above_surface!" begin
    r = zeros(3, 2)
    height = 10
    QuantisedDiatomic.position_above_surface!(r, height, InfiniteCell())
    @test r == [0 0; 0 0; height height]
end

@testset "velocity_from_energy" begin
    v = QuantisedDiatomic.velocity_from_energy([100, 100], 10u"eV")
    v2 = QuantisedDiatomic.velocity_from_energy([100, 100], austrip(10u"eV"))
    @test v ≈ v2
end

@testset "apply_translational_impulse!" begin
    masses = [100, 100]
    r = rand(3, 2) .+ [0 2; 0 2; 0 2]
    v = rand(3, 2) ./ 100
    v_before = copy(v)
    e = 10u"eV"
    direction = [0, 0.5, 0.5]
    QuantisedDiatomic.apply_translational_impulse!(v, masses, e, direction)
    @test v_before[1,:] == v[1,:]
    @test all(v_before[2:3,:] .< v[2:3,:])

    νi, Ji = quantise_diatomic(sim, v, r)
    QuantisedDiatomic.apply_translational_impulse!(v, sim.atoms.masses, 1, rand(3))
    νf, Jf = quantise_diatomic(sim, v, r)
    @test νi ≈ νf
    @test Ji ≈ Jf
end


text:
@noinline function dostuff(A; B = I)
    StrideArrays.@avx for i in eachindex(A)
        A[i] += 1
    end
    C = A * B
    s = zero(eltype(C))
    StrideArrays.@avx for i in eachindex(C)
        s += C[i]
    end
    s
end
function gc_preserve_test()
    A = @StrideArray rand(8, $(2<<2));
    B = @StrideArray rand($(1<<3), 8);
    @gc_preserve dostuff(A, B = B)
end

@testset "Miscellaneous" begin
    x = @StrideArray rand(127);
    @test maximum(abs, x) == maximum(abs, Array(x))
    y = @StrideArray rand(3);
    @test maximum(abs, y) == maximum(abs, Array(y))
    @test iszero(@allocated gc_preserve_test())

    A = @StrideArray rand($(5<<1),$(1<<3));
    @test StrideArrays.size(A) === (StaticInt(10), StaticInt(8))
    A_u = view(A, StaticInt(1):StaticInt(6), :)
    A_l = view(A, StaticInt(7):StaticInt(10), :)
    @test A == @inferred(vcat(A_u, A_l))
    @test StrideArrays.size(A) === StrideArrays.size(vcat(A_u, A_l))

    # On 1.5 tests fail if you don't do this first.
    @test pointer(view(A, 1:StaticInt(6), :)) == pointer(A)
    A_u = view(A, 1:StaticInt(6), :);
    A_l = view(A, StaticInt(7):StaticInt(10), :);
    @test A == @inferred(vcat(A_u, A_l))
    @test StrideArrays.size(A) == StrideArrays.size(vcat(A_u, A_l))
    @test StrideArrays.size(A) !== StrideArrays.size(vcat(A_u, A_l))

    Aa = Array(A)
    @test sum(A) ≈ sum(Aa)
    @test maximum(A) == maximum(Aa) == maximum(abs, @. -A)
    @test mapreduce(abs2, +, A) ≈ mapreduce(abs2, +, Aa)

    @test A[1,1] == A[1] == vec(A)[1]
    A[1,1] = 3
    @test A[1,1] == A[1] == vec(A)[1] == 3
    A[2] = 4
    @test A[2,1] == A[2] == vec(A)[2] == 4
    vec(A)[3] = 5
    @test A[3,1] == A[3] == vec(A)[3] == 5

    @test ArrayInterface.stride_rank(StrideArrays.similar_layout(A')) === ArrayInterface.stride_rank(StrideArrays.similar_layout(Aa')) === (StaticInt(2),StaticInt(1))
    @test ArrayInterface.contiguous_axis(StrideArrays.similar_layout(A')) === ArrayInterface.contiguous_axis(StrideArrays.similar_layout(Aa')) === StrideArrays.StaticInt{2}()

    let B = @StrideArray rand(10,10)
        @test exp(B) ≈ exp(Array(B))
    end
end


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here

# a custom shape that could be wrapped in a function */
x         = 25.6;        # parameters like cairo_rectangle */
y         = 25.6;
xw         = 204.8;
yw        = 204.8;
aspect        = 1.0;     # aspect ratio */
corner_radius = yw / 10.0;   #* and corner curvature radius */

radius = corner_radius / aspect;
degrees = pi / 180.0;

new_sub_path(cr);
arc(cr, x + xw - radius, y + radius, radius, -90 * degrees, 0 * degrees);
arc(cr, x + xw - radius, y + yw - radius, radius, 0 * degrees, 90 * degrees);
arc(cr, x + radius, y + yw - radius, radius, 90 * degrees, 180 * degrees);
arc(cr, x + radius, y + radius, radius, 180 * degrees, 270 * degrees);
close_path(cr);

set_source_rgb(cr, 0.5, 0.5, 1);
fill_preserve(cr);
set_source_rgba(cr, 0.5, 0, 0, 0.5);
set_line_width(cr, 10.0);
stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_rounded_rectangle.png");


text:
using TrussMorph
using Test

# The testing setup
# Here we split different categories into tests
# Each category is a separate file
# By placing `@testset` on the file, it will run all of the tests there,
# and report back when the entire file is complete

# Optional: If you @testset begin ... end around the entire testing setup
# then all of your tests will run, and the failures will be counted
# and printed at the end. There is a balance between the completeness of
# the test results and the time it takes to run tests!

# Optional: tic() ... toc() around the tests as a quick way to keep track of
# large performance regressions.

@testset "TrussMorph.jl" begin
    include("parse_json_tests.jl")
    # include("stiffness_matrix_tests.jl")
    # include("viz_tests.jl")
end


text:
# Copyright © 2021 Alexander L. Hayes
# Apache 2.0 License

"""Convert (array-based) propositional datasets to a relational/ILP
representation.
"""

include("types.jl")

function _is_multiclass(y::Vector{Int64})::Bool
    return length(unique(y)) > 2
end

"""
    from_vector(X::Matrix{Int}, y::Vector{Int}, names::Union{Vector{String}, Nothing} = nothing)

Convert a classification dataset to an ILP representation.
"""
function from_vector(X::Matrix{Int64}, y::Vector{Int64}, names::Union{Vector{String}, Nothing} = nothing)

    # Make sure the X and y are valid.
    @assert size(y)[1] == size(X)[2]

    if names == nothing
        names = ["v$(i)" for i in 1:size(X)[1] + 1]
    end

    # Make sure the names vector makes sense.
    @assert length(names) - 1 == size(X)[1]

    pos, neg, facts = String[], String[], String[]

    is_multiclass = _is_multiclass(y)

    if is_multiclass

        for (i, row) in enumerate(y)
            push!(pos, "$(last(names))(id$(i),$(row)).")
        end

    else

        for (i, row) in enumerate(y)
            if Bool(row)
                push!(pos, "$(last(names))(id$(i)).")
            else
                push!(neg, "$(last(names))(id$(i)).")
            end
        end
    end

    for (i, col) in enumerate(eachcol(X))
        var = names[i]
        facts = vcat(facts, ["$(var)(id$(j),$(row))." for (j, row) in enumerate(col)])
    end

    modes = ["$(name)(+id,#var$(name))." for name in names[1:end-1]]

    if is_multiclass
        push!(modes, "$(last(names))(+id,#classlabel).")
    else
        push!(modes, "$(last(names))(+id).")
    end

    return RelationalDataset((pos, neg, facts)), modes
end

"""
    from_vector(X::Matrix{Int}, y::Vector{Float64}, names::Union{Vector{String}, Nothing} = nothing)

Convert a regression dataset to an ILP representation.
"""
function from_vector(X::Matrix{Int64}, y::Vector{Float64}, names::Union{Vector{String}, Nothing} = nothing)

    # Make sure the X and y are valid.
    @assert size(y)[1] == size(X)[2]

    if names == nothing
        names = ["v$(i)" for i in 1:size(X)[1] + 1]
    end

    # Make sure the names vector makes sense.
    @assert length(names) - 1 == size(X)[1]

    pos, neg, facts = String[], String[], String[]

    for (i, row) in enumerate(y)
        push!(pos, "regressionExample($(last(names))(id$(i)),$(row)).")
    end

    for (i, col) in enumerate(eachcol(X))
        var = names[i]
        facts = vcat(facts, ["$(var)(id$(j),$(row))." for (j, row) in enumerate(col)])
    end

    modes = ["$(name)(+id,#var$(name))." for name in names[1:end-1]]
    push!(modes, "$(last(names))(+id).")

    return RelationalDataset((pos, neg, facts)), modes
end


text:
"""
Get array of sensor ids from folder names (only tested for DS18B20 temperature sensor)
"""
function sensorids(sensordir::String="/sys/bus/w1/devices")
    sids = filter(x->startswith(x, "28-"), readdir(sensordir))
    return sids
end

"""
Get the temperature from a particular sensor based on unique id
"""
function readtemp(sensorid::String, sensordir::String="/sys/bus/w1/devices")
    sid = joinpath(sensordir, sensorid)
    p = joinpath(sid, "w1_slave")
    isfile(p) || error("Couldn't find sensor $sensorid at '$sensordir'")

    readings = readlines(p)
    if endswith(readings[1], "YES")
        m = match(r"t=(\d+)$", readings[2])
        ctemp = parse(Float64, m.captures[1]) / 1000.
        ctemp == 85. && warn("Sensor $sensorid may be improperly connected (temp = 85.000°C)")
        ftemp = ctemp * 1.8 + 32.

        return ctemp, ftemp
    else
        @warn "Sensor $sensorid found, but couldn't get temp"
    end
end


text:
# -*- coding: utf-8 -*-
# + {}
module Transformations

using SemanticModels.Parsers
using SemanticModels.ModelTools

abstract type Transformation end

struct ExprTransformation <: Transformation
    func::Symbol
    state::Symbol
    expr::Expr
end

struct ConcatTransformation
  seq::Vector{Transformation}
end

promote(::ConcatTransformation, ::Transformation) = ConcatTransformation
convert(::ConcatTransformation, t::Transformation) = ConcatTransformation([t])

function ∘(f::ConcatTransformation, g::ConcatTransformation)
    append!(g.seq, f.seq)
    return g
end

function ∘(f::Transformation, g::ConcatTransformation)
    append!(g.seq, [f])
    return g
end

function (t::ExprTransformation)(m::AbstractProblem)
    getfield(ModelTools, t.func)(m, ExpStateTransition(t.state, t.expr))
    return m
end

function (f::ConcatTransformation)(m::AbstractProblem)
    return foldl((m, t)->t(m), [m; f.seq])
end

expr = parsefile("../examples/agentbased.jl")
m = model(ExpStateModel, expr)
T = ConcatTransformation([])
T = ExprTransformation(:put!, :D, :((x...)->:D)) ∘ T
T = ExprTransformation(:replace!, :I, :((x...)->rand(Bool) ? :R : :D)) ∘ T
@show m
@show T.seq
T(m)
@show m
#M = eval(Expr(m))
#sol = M.solve()

end
# -


text:
module par_module
export f
using Roots
function f(k,b,c)
    b*[k k^2; k^2 k]+c*[1 1; 1 1]
end
end


text:
using DrWatson
@quickactivate "StatisticalRethinkingDHMC"
DrWatson.greet()


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2020- Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

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

"""
    steady_state_solve()

Fill in an awesome description here
"""
function steady_state_solve(data_dictionary::Dict{String,Any})

    # Get the initial conditions -
    initial_condition_array = data_dictionary["initial_condition_array"];

    # calculate the steady-state soln -
    steady_state_prob = SteadyStateProblem(balances, initial_condition_array, data_dictionary)
    steady_state_soln  = solve(steady_state_prob, DynamicSS(AutoTsit5(Rosenbrock23(autodiff=false));
        abstol=1e-8,reltol=1e-6,tspan=Inf))

    # return -
    return abs.(steady_state_soln.u)
end

"""
    dynamic_solve()

Fill in an awesome description here
"""
function dynamic_solve(TSTART::Float64,TSTOP::Float64,Ts::Float64,
    data_dictionary::Dict{String,Any})

    # Setup time scale -
    time_span = (TSTART,TSTOP)

    # grab the initial condition from the data dictionary -
    initial_condition_array = data_dictionary["initial_condition_array"];

    # build the ODE problem object -
    problem_object = ODEProblem(balances, initial_condition_array, time_span, data_dictionary)

    # solve -
    solution = solve(problem_object, AutoTsit5(Rosenbrock23(autodiff=false)), reltol=1e-8,abstol=1e-8,saveat=Ts)

    # pull solution apart -
    T = solution.t

    # initialize the state array -
    number_of_times_steps = length(T)
    number_of_states = length(initial_condition_array)
    X = zeros(number_of_times_steps,number_of_states)
    for step_index=1:number_of_times_steps

        # grab the solution 0
        soln_array = solution.u[step_index]
        for state_index = 1:number_of_states
            X[step_index, state_index] = soln_array[state_index]
        end
    end

    # return -
    return (T,X)
end

max_stars_repo_name: varnerlab/N-Compartment-Chip-Model

max_issues_repo_name: varnerlab/N-Compartment-Chip-Model

max_forks_repo_name: varnerlab/N-Compartment-Chip-Model


text:
# interfacing with Random standard library

using Random

using GPUArrays

mutable struct RNG <: Random.AbstractRNG
    handle::rocrand_generator
    typ::rocrand_rng_type

    function RNG(typ=ROCRAND_RNG_PSEUDO_DEFAULT)
        handle = Ref{rocrand_generator}()
        rocrand_create_generator(handle, typ)
        obj = new(handle[], typ)
        finalizer(unsafe_destroy!, obj)
        return obj
    end
end

function unsafe_destroy!(rng::RNG)
    rocrand_destroy_generator(rng.handle)
end

Base.unsafe_convert(::Type{rocrand_generator}, rng::RNG) = rng.handle

## seeding
function Random.seed!(rng::RNG, seed=Base.rand(UInt64), offset=0)
    rocrand_set_seed(rng, seed)
    rocrand_set_offset(rng, seed)
    rocrand_initialize_generator(rng)
    return
end

Random.seed!(rng::RNG, ::Nothing) = Random.seed!(rng)

## in-place

# uniform
const UniformType = Union{Type{UInt32}, Type{Cuchar}, Type{UInt16}, Type{Float16}, Type{Float32}, Type{Float64}}
const DUniformArray = ROCArray{<:Union{UInt32, Cuchar, UInt16, Float16, Float32,Float64}}
for (f,T) in ((:rocrand_generate, :UInt32), (:rocrand_generate_char,:Cuchar),
              (:rocrand_generate_short, :UInt16), (:rocrand_generate_uniform, :Float32),
              (:rocrand_generate_uniform_double, :Float64), (:rocrand_generate_uniform_half, :Float16))
    @eval begin
        function Random.rand!(rng::RNG, A::ROCArray{$(T)})
            $(f)(rng, A, length(A))
            return A
        end
    end
end

# some functions need pow2 lengths: use a padded array and copy back to the original one
function inplace_pow2(A, f)
    len = length(A)
    if len > 1 && ispow2(len)
        f(A)
    else
        padlen = max(2, nextpow(2, len))
        B = similar(A, padlen)
        f(B)
        copyto!(A, 1, B, 1, len)
        # TODO: why is unsafe_free used in CUDA?
        # here it causes "Release of dead AMDGPU.Mem.Buffer" error
        #AMDGPU.unsafe_free!(B)
    end
    A
end

# normal
const NormalType = Union{Type{Float16}, Type{Float32}, Type{Float64}}
const NormalArray = ROCArray{<:Union{Float16, Float32,Float64}}
function Random.randn!(rng::RNG, A::ROCArray{Float32}; mean=0, stddev=1)
    inplace_pow2(A, B->rocrand_generate_normal(rng, B, length(B), mean, stddev))
    return A
end
function Random.randn!(rng::RNG, A::ROCArray{Float64}; mean=0, stddev=1)
    inplace_pow2(A, B->rocrand_generate_normal_double(rng, B, length(B), mean, stddev))
    return A
end
function Random.randn!(rng::RNG, A::ROCArray{Float16}; mean=0, stddev=1)
    inplace_pow2(A, B->rocrand_generate_normal_half(rng, B, length(B), mean, stddev))
    return A
end

# log-normal
const LognormalType = Union{Type{Float16}, Type{Float32},Type{Float64}}
const LognormalArray = ROCArray{<:Union{Float16, Float32, Float64}}
function rand_logn!(rng::RNG, A::ROCArray{Float32}; mean=0, stddev=1)
    inplace_pow2(A, B->rocrand_generate_log_normal(rng, B, length(B), mean, stddev))
    return A
end
function rand_logn!(rng::RNG, A::ROCArray{Float64}; mean=0, stddev=1)
    inplace_pow2(A, B->rocrand_generate_log_normal_double(rng, B, length(B), mean, stddev))
    return A
end
function rand_logn!(rng::RNG, A::ROCArray{Float16}; mean=0, stddev=1)
    inplace_pow2(A, B->rocrand_generate_log_normal_half(rng, B, length(B), mean, stddev))
    return A
end

# poisson
const PoissonType = Union{Type{Cuint}}
const PoissonArray = ROCArray{Cuint}
function rand_poisson!(rng::RNG, A::ROCArray{Cuint}; lambda=1)
    rocrand_generate_poisson(rng, A, length(A), lambda)
    return A
end

## out of place

# some functions need pow2 lengths: construct a compatible array and return part of it
function outofplace_pow2(shape, ctor, f)
    len = prod(shape)
    if ispow2(len)
        A = ctor(shape)
        f(A)
    else
        padlen = max(2, nextpow(2, len))
        A = ctor(padlen)
        f(A)
        B = reshape(A[1:len], shape)
        return B
    end
end

# arrays
Random.rand(rng::RNG, T::UniformType, dims::Dims) =
    Random.rand!(rng, ROCArray{T}(undef, dims))
Random.randn(rng::RNG, T::NormalType, dims::Dims; kwargs...) =
    outofplace_pow2(dims, shape->ROCArray{T}(undef, dims), A->randn!(rng, A; kwargs...))
rand_logn(rng::RNG, T::LognormalType, dims::Dims; kwargs...) =
    outofplace_pow2(dims, shape->ROCArray{T}(undef, dims), A->rand_logn!(rng, A; kwargs...))
rand_poisson(rng::RNG, T::PoissonType, dims::Dims; kwargs...) =
    rand_poisson!(rng, ROCArray{T}(undef, dims); kwargs...)

# specify default types
Random.rand(rng::RNG, dims::Dims; kwargs...) = Random.rand(rng, Float32, dims; kwargs...)
Random.randn(rng::RNG, dims::Dims; kwargs...) = Random.randn(rng, Float32, dims; kwargs...)
rand_logn(rng::RNG, dims::Dims; kwargs...) = rand_logn(rng, Float32, dims; kwargs...)
rand_poisson(rng::RNG, dims::Dims; kwargs...) = rand_poisson(rng, Cuint, dims; kwargs...)

# support all dimension specifications
Random.rand(rng::RNG, T::UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(rng, T, Dims((dim1, dims...)))
Random.randn(rng::RNG, T::NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(rng, T, Dims((dim1, dims...)); kwargs...)
rand_logn(rng::RNG, T::LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn(rng, T, Dims((dim1, dims...)); kwargs...)
rand_poisson(rng::RNG, T::PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson(rng, T, Dims((dim1, dims...)); kwargs...)

# rand_logn! and rand_poisson! without specified rng
rand_logn!(A::rocRAND.LognormalArray; kwargs...) = rand_logn!(default_rng(), A; kwargs...)
rand_poisson!(A::rocRAND.PoissonArray; kwargs...) = rand_poisson!(default_rng(), A; kwargs...)

rand_logn(T::rocRAND.LognormalType, dims::Dims; kwargs...) = rand_logn(default_rng(), T, dims; kwargs...)
rand_poisson(T::rocRAND.PoissonType, dims::Dims; kwargs...) = rand_poisson(default_rng(), T, dims; kwargs...)

rand_logn(T::rocRAND.LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn(default_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_poisson(T::rocRAND.PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson(default_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_logn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn!(ROCArray{T}(undef, dim1, dims...); kwargs...)
rand_poisson(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson!(ROCArray{T}(undef, dim1, dims...); kwargs...)
rand_logn(dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn(default_rng(), Dims((dim1, dims...)); kwargs...)
rand_poisson(dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson(default_rng(), Dims((dim1, dims...)); kwargs...)
rand_logn(T::Type, dims::Dims; kwargs...) = rand_logn!(ROCArray{T}(undef, dims...); kwargs...)
rand_poisson(T::Type, dims::Dims; kwargs...) = rand_poisson!(ROCArray{T}(undef, dims...); kwargs...)
rand_logn!(A::ROCArray; kwargs...) =
    error("AMDGPU.jl does not support generating lognormally-distributed random numbers of type $(eltype(A))")
rand_poisson!(A::ROCArray; kwargs...) =
    error("AMDGPU.jl does not support generating Poisson-distributed random numbers of type $(eltype(A))")


text:
using YaoHIR: X, Y, Z, H

const naive = @theory begin
    Chain(Gate($X, locs), Chain(Gate($X, locs), C)) => C
    Chain(Gate($X, locs), Gate($X, locs)) => Chain()
    Chain(Gate($Y, locs), Chain(Gate($Y, locs), C)) => C
    Chain(Gate($Y, locs), Gate($Y, locs)) => Chain()
    Chain(Gate($Z, locs), Chain(Gate($Z, locs), C)) => C
    Chain(Gate($Z, locs), Gate($Z, locs)) => Chain()
    Chain(Gate($H, locs), Chain(Gate($H, locs), C)) => C
    Chain(Gate($H, locs), Gate($H, locs)) => Chain()
    Chain(Chain(C)) => Chain(C)
    Chain(C, Chain()) => C
end


text:
function make_table(cur_rows::Int, cur_cols::Int, cur_hidden_rows::Int=0; class_names::AbstractString="")
  cur_row_array = Array{Node}(0)

  total_rows = cur_rows + cur_hidden_rows

  for cur_row in total_rows:-1:1

    cur_cell_array = Array{Node}(0)

    for cur_col in 1:cur_cols
      push!(
        cur_cell_array,
        Node(:td)
      )
    end

    cur_class = "cs-row-$(cur_row)"

    is_hidden_row = ( cur_row > cur_rows )

    if is_hidden_row
      cur_class *= " hidden"
    else
      cur_class *= " cs-light-grey"
    end

    push!(
      cur_row_array,
      Node(
        :tr,
        cur_cell_array...,
        attributes=Dict(:class => cur_class)
      )
    )

  end

  cur_table = Node(
    :table,
    cur_row_array...,
    attributes=Dict(:class => "cs-tetris-table z-depth-1-half $(class_names)")
  )

  cur_table
end


text:
@testset "AbstractTriangular ambiguity methods" begin
    ut = UpperTriangular([1 0 0; 0 2 0; 0 0 3])
    v = [2,3,4]
    tv = v.'

    @test (tv*ut)::TransposedVector == [2 6 12]
    @test_throws Exception ut*tv

    @test (tv*ut.')::TransposedVector == [2 6 12]
    @test (ut*tv.')::Vector == [2,6,12]

    @test (ut.'*tv.')::Vector == [2,6,12]
    @test_throws Exception tv.'*ut.'

    @test_throws Exception ut.'*tv
    @test_throws Exception tv.'*ut

    @test (tv*ut')::TransposedVector == [2 6 12]
    @test (ut*tv')::Vector == [2,6,12]

    @test_throws Exception tv'*ut'
    @test (ut'*tv')::Vector == [2,6,12]

    @test_throws Exception ut'*tv
    @test_throws Exception tv'*ut

    @test (tv/ut)::TransposedVector ≈ [2/1  3/2  4/3]
    @test (tv/ut.')::TransposedVector ≈ [2/1  3/2  4/3]
    @test (tv/ut')::TransposedVector ≈ [2/1  3/2  4/3]

    @test_throws Exception ut\tv

end


text:
"""
Return a decoded JSON file.
"""
function read_json(filename::AbstractString)
    open(filename, "r") do io
        JSON.parse(io)
    end
end

"""
Return a DataFrame from a CSV file.
"""
function read_dataframe(filename::AbstractString)
    open(filename, "r") do io
        DataFrames.DataFrame(CSV.File(io))
    end
end

"""
Return the SHA 256 hash of a file.
"""
function compute_sha256(filename::AbstractString)
    return open(filename) do io
        return bytes2hex(SHA.sha256(io))
    end
end

"""
Return the key for the given value
"""
function find_key_with_value(d, value)
    for (k, v) in d
        v == value && return k
    end
    error("dict does not have value == $value")
end

function compute_file_hash(path::String, files::Vector{String})
    data = Dict("files" => [])
    for file in files
        file_path = joinpath(path, file)
        # Don't put the path in the file so that we can move results directories.
        file_info = Dict("filename" => file, "hash" => compute_sha256(file_path))
        push!(data["files"], file_info)
    end

    open(joinpath(path, HASH_FILENAME), "w") do io
        write(io, JSON.json(data))
    end
end

function compute_file_hash(path::String, file::String)
    return compute_file_hash(path, [file])
end

function read_file_hashes(path)
    data = open(joinpath(path, HASH_FILENAME), "r") do io
        JSON.parse(io)
    end

    return data["files"]
end

# this ensures that the timestamp is not double shortened
function find_variable_length(es::Dict, e_list::Array)
    return size(es[Symbol(splitext(e_list[1])[1])], 1)
end

"""
    check_file_integrity(path::String)

Checks the hash value for each file made with the file is written with the new hash_value to verify the file hasn't been tampered with since written

# Arguments

  - `path::String`: this is the folder path that contains the results and the check.sha256 file
"""
function check_file_integrity(path::String)
    matched = true
    for file_info in read_file_hashes(path)
        filename = file_info["filename"]
        @info "checking integrity of $filename"
        expected_hash = file_info["hash"]
        actual_hash = compute_sha256(joinpath(path, filename))
        if expected_hash != actual_hash
            @error "hash mismatch for file" filename expected_hash actual_hash
            matched = false
        end
    end

    if !matched
        throw(
            IS.HashMismatchError(
                "The hash value in the written files does not match the read files, results may have been tampered.",
            ),
        )
    end
end

to_namedtuple(val) = (; (x => getfield(val, x) for x in fieldnames(typeof(val)))...)

convert_for_path(x::Dates.DateTime) = replace(string(x), ":" => "-")


text:
using Resonance

profiles = String[]
for (root, dirs, files) in walkdir("/grace/sequencing/processed/mgx")
    contains(root, "metaphlan") || continue
    contains(root, "links") && continue
    filter!(f-> contains(f, r"^FG.+profile"), files)
    append!(profiles, joinpath.(Ref(root), files))
end

species = metaphlan_profiles(profiles, :species)

sns = map(samplenames(species)) do s
    replace(s, r"_S\d+_profile"=>"")
end

idx = unique([findfirst(==(x), sns) for x in sns])

species = CommunityProfile(abundances(species)[:, idx], features(species), map(samplenames(species)[idx]) do s
    MicrobiomeSample(replace(s, r"_S\d+_profile"=>""))
end)

CSV.write("data/species.csv", species)

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis