text:
using ImplicitEquations, Plots

include("project2.jl")
include("helpers.jl")

probname = "simple1"
n_seeds = 3
seed = 50
num_iter = 30

prob = PROBS[probname]

x0 = prob.x0()

# TODO: Contour plot for Rosenbrock function:
# use range of (-3.0, 3.0) for x1, x2
# z(x1, x2) = f([x1, x2])
x1 = range(-3.0,stop=3.0,length=50)
x2 = range(-3.0,stop=3.0,length=50)
funcn = zeros(50,50)
# c1 = zeros(50,50)
# c2 = zeros(50,50)
for i = 1:50
    for j = 1:50
        global funcn
        funcn[i,j] = prob.f([x1[i], x2[j]])
        # c1[i,j] = x1[i] + x2[j]^2 - 1
        # c2[i,j] = -x1[i] -x2[j]
    end
end

c1(x, y) = x + y^2 - 1
c2(x, y) = -x - y
p1 = contour(x1, x2, funcn, colorbar = false,
    c = cgrad(:viridis, rev = true), legend = false, xlims = (-3, 3), ylims = (-3, 3),
    xlabel = "x₁", ylabel = "x₂", aspectratio = :equal)
plot!(c2 ≤ 0)
# contour!(x1, x2, c1, levels = [0], c = :red, fillrange(0,0))
# contour!(x1, x2, c2, levels = [0], c = :black)
# contour!(x1, x2, c1, levels = [0], c = :gray)
# plot!(xhist_nseeds[:,1,1], xhist_nseeds[:,2,1])
# plot!(xhist_nseeds[:,1,2], xhist_nseeds[:,2,2])
# plot!(xhist_nseeds[:,1,3], xhist_nseeds[:,2,3])
# savefig(p1, "rosenbrock_trajectories.pdf")

# Plot x_history on top of function contour

# # CONVERGENCE PLOTS
# fhist_nseeds = zeros(num_iter, n_seeds)
#
# for i = 1:n_seeds
#     Random.seed!(seed+i)
#     _, fhist_nseeds[:,i] = optimize_plot_data(prob.f, prob.g, prob.x0(), prob.n, probname)
# end
#
# p = plot(1:num_iter, fhist_nseeds, xlabel = "Iterations", ylabel = "Function value", title = "Rosenbrock Convergence Plot", label = ["Seed 1" "Seed 2" "Seed 3"], lw = 2)
# savefig(p, "rosenbrock_convergence.pdf")


text:
using Compat
import Base: Forward
import SortingAlgorithms: RADIX_MASK, RADIX_SIZE, uint_mapping
using FastGroupBy
using BenchmarkTools

using SortingAlgorithms, Base.Order, Compat, IndexedTables, DataFrames
import Base: ht_keyindex, rehash!, _setindex!
import SortingAlgorithms: uint_mapping, RADIX_SIZE, RADIX_MASK
import DataFrames: DataFrame, AbstractDataFrame
import IndexedTables: IndexedTable, column
import PooledArrays.PooledArray
import CategoricalArrays.CategoricalArray
import SplitApplyCombine.groupreduce

using StatsBase

T = Int
S = Int

N  = 1_000_000;
K = 100;
by = nothing;
val = nothing;
gc();
srand(1);
by = rand(Int32(1):Int32(round(N/K)), N);
val = rand(Int32(1):Int32(5), N)

#
# @time sumby_lessmem_chain(by, val);
#
# by = nothing;
# val = nothing;
# gc();
# srand(1);
# @time val = rand(Int32(1):Int32(5), N);
# @time by = rand(Int32(1):Int32(round(N/K)), N);
# @time sumby(by, val);
#
# function abc()
#   by1, val1 = copy(by), copy(val)
#   @elapsed sumby_lessmem_chain(by1, val1)
# end
#
# function def()
#   by1, val1 = copy(by), copy(val)
#   @elapsed sumby(by1, val1)
# end

# a1 = mean([abc() for i = 1:5][2:end])
# b1 = mean([def() for i = 1:5][2:end])
# 1 - a1/b1
2+2

ext: jl


text:
using AtomicOpt
using LinearAlgebra
using Random
using Test

import AtomicOpt: LBFGSBQuad, Quadratic, boxquad, grad

"""Project into the box [bl, bu]."""
proj(x, bl, bu) = max.(min.(x, bu), bl)

@testset "Box-constrained LS via projected Newton" begin
    m, n = 10, 3
    A = randn(m, n); x0 = randn(n); b = A*x0
    bl = zeros(n); bu = ones(n)
    x = randn(n)

    obj(x) = 0.5*norm(A*x)^2 - dot(A'b,x)

    Q = Quadratic(A'*A, A'b)
    @test obj(x) ≈ Q(x)

end

@testset "Box-constrained LS via LBFGS-B" begin
    m, n = 100, 30
    A = randn(m, n); x0 = randn(n); b = A*x0
    bl = rand(n) .- 0.5; bu = rand(n) .+ 0.5

    Q = Quadratic(A'*A, A'b)
    lbfgsq = LBFGSBQuad(Q, bl, bu)
    x = boxquad(lbfgsq, factr=1.0)
    @test norm( x - proj(x-grad(Q,x), bl, bu)) ≤ 1e-6*norm(x)
end


text:
using Flux, MLDataPattern, Mill, JsonGrinder, JSON, Statistics, IterTools, ThreadTools
using JsonGrinder: suggestextractor, ExtractDict
using Mill: reflectinmodel
using StatsBase: sample

###############################################################
# start by loading all samples
###############################################################
samples = Vector{Dict}(open(JSON.parse, "data/mutagenesis/data.json"))

JSON.print(samples[3],2)

metadata = open(JSON.parse, "data/mutagenesis/meta.json")
labelkey = metadata["label"]
test_num = metadata["test_samples"]
minibatchsize = 100
iterations = 5_000
neurons = 20 		# neurons per layer

targets = map(i -> i[labelkey], samples)
foreach(i -> delete!(i, labelkey), samples)

train_indices = 1:length(samples)-test_num
test_indices = length(samples)-test_num+1:length(samples)

#####
#  Create the schema and extractor
#####
sch = JsonGrinder.schema(samples)
extractor = suggestextractor(sch)

#####
#  Convert samples to Mill structure and extract targets
#####
data = tmap(extractor, samples)
labelnames = unique(targets)

#####
#  Create the model
#####
model = reflectinmodel(sch, extractor,
	k -> Dense(k, neurons, relu),
	d -> SegmentedMeanMax(d),
	fsm = Dict("" => k -> Dense(k, length(labelnames))),
)

#####
#  Train the model
#####
function minibatch()
	idx = sample(1:length(data[train_indices]), minibatchsize, replace = false)
	reduce(catobs, data[idx]), Flux.onehotbatch(targets[idx], labelnames)
end

accuracy(x,y) = mean(labelnames[Flux.onecold(model(x).data)] .== y)

trainset = reduce(catobs, data[train_indices])
testset = reduce(catobs, data[test_indices])

cb = () -> begin
	train_acc = accuracy(trainset, targets[train_indices])
	test_acc = accuracy(testset, targets[test_indices])
	println("accuracy: train = $train_acc, test = $test_acc")
end
ps = Flux.params(model)
loss = (x,y) -> Flux.logitcrossentropy(model(x).data, y)
Flux.Optimise.train!(loss, ps, repeatedly(minibatch, iterations), ADAM(), cb = Flux.throttle(cb, 2))

###############################################################
#  Classify test set
###############################################################

probs = softmax(model(testset).data)
o = Flux.onecold(probs)
pred_classes = labelnames[o]

print(mean(pred_classes .== targets[test_indices]))
# we see the accuracy is around 79% on test set

#predicted classes for test set
print(pred_classes)
#gt classes for test set
print(targets[test_indices])
# probabilities for test set
print(probs)


text:
"""
    Polynomials functions for spline interpolation: basis functions for b-splines
    and bernstein polynomials for bezier curves.

    Most polynomials have two implementations: one in which the polynomial is
    computed over the entier parameter range, the other for evaluating the
    polynomial at a single parameter value.
"""
module Polynomials

    # ---------------------------------------------------------------------------- #
    #                           B-Spline BASIS FUNCTIONS                           #
    # ---------------------------------------------------------------------------- #
    """
        N_0(τ, k; i=0)

    Zero-th order basis function for b-splines for the i-th knot.
        for t∈[k[i], k[i+1]] = 1
        else: 0

    The additional "+1" index in the code is because in the bspline maths
    the knots are indexed starting from 0.
    """
    function N_0(τ; k::AbstractArray, i::Int=0)
        if i+2 < length(k)
            return k[i+1] .<= τ .< k[i+2]
        else
            return k[i+1] .<= τ .<= k[i+2]
        end
    end

    """
        ω(τ, j)

    Function used within N_D to compute the factors for the linear interpolation
    of the two lover level basis functions
    """
    function ω(τ, j, d, k)
        Δ = k[j+d]-k[j]
        @. Δ!=0 ? (τ - k[j])/Δ : zero(τ)
    end

    """
        N_D(τ, k; i=0, d=0)

    D-th order basis function for b-splines for the i-th knot.
    Built recursively on lower order basis functions.
    See: https://www.geometrictools.com/Documentation/BSplineCurveLeastSquaresFit.pdf
    and: https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-basis.html

    The "j=i+1" index in the code is because in the bspline maths
    the knots are indexed starting from 0.
    """
    function N_D(τ; k::AbstractArray, i::Int=0, d::Int=0)
        ω(τ, i+1, d, k) .* N_basis(τ; k=k, i=i, d=d-1) .+ (1 .- ω(τ, i+1+1, d, k)) .* N_basis(τ; k=k, i=i+1, d=d-1)
    end

    """
        N(τ, k; i=0, d=0)

    B-spline basis function for the i-th knot and d-th order.
    Calls either N_0 or N_D depending on the value of d.
    """
    N_basis(τ; k::AbstractArray, i::Int=0, d::Int=0) = (d == 0) ? N_0(τ; k, i=i) : N_D(τ; k, i=i, d=d)

    # ---------------------------------------------------------------------------- #
    #                             BERNSTEIN POLYNOMIALS                            #
    # ---------------------------------------------------------------------------- #
    """
        bernstein(t::Float64; i::Int, n::Int)

    Evaluate the bernstein polynomial at parameter value `t` given the index `i` and the number
    of polynomials `n`.
    """
    bernstein(τ; i::Int, n::Int) = @. binomial(n, i) * τ^i * (1 - τ)^(n-i)

end

max_stars_repo_name: FedeClaudi/Splines

max_issues_repo_name: FedeClaudi/Splines

max_forks_repo_name: FedeClaudi/Splines


text:
using Revise

using JDF, CSV, DataFrames, Blosc, Base.GC, BenchmarkTools, Serialization

# use 12 threads
Blosc.set_num_threads(6)
# GC.gc()
@time a = CSV.read(
    "C:/data/Performance_All/Performance_2003Q1.txt",
    delim = '|',
    header = false
);
GC.gc()

t = time()
@time metadatas = psavejdf(a, "c:/data/large8.dir.jdf");
time() - t
GC.gc()

GC.gc()
@time metadatas = savejdf(a, "c:/data/large8.jdf");
GC.gc()
serialize("c:/data/large8.meta", metadatas)

using Revise, JDF, DataFrames, Serialization
metadatas = deserialize("c:/data/large8.meta")
@time a2 = loadjdf("c:/data/large8.jdf", metadatas);
# before=Base.summarysize(a2)

@time type_compress!(a2, verbose = true)

# after=Base.summarysize(a2)
#
# before/after

GC.gc()
@time metadatas = savejdf(a2, "c:/data/large8c.jdf");
GC.gc()
serialize("c:/data/large8c.meta", metadatas)

using Serialization, JDF
metadatas = deserialize("c:/data/large8c.meta")
@time a2 = loadjdf("c:/data/large8c.jdf", metadatas);

GC.gc()

all(names(a) .== names(a2))
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)]))


text:
using ConicBenchmarkUtilities
using Mosek
using MathProgBase

# Usage: process_output.jl [--noconicsolve] RESULTFOLDERS... OUTPUT.csv

@assert length(ARGS) >= 2
NOCONIC = false
if ARGS[1] == "--noconicsolve"
    NOCONIC = true
    shift!(ARGS)
end
resultfiles = []
for i in 1:length(ARGS)-1
    append!(resultfiles,[joinpath(pwd(),ARGS[i],f) for f in readdir(joinpath(pwd(), ARGS[i]))])
end

fd = open(joinpath(pwd(), ARGS[end]), "w")

# Process into a CSV file with columns:
println(fd,"solver,instance,filename,sense,timelimit,solver_time,total_time,status,newstatus,objval,objbound,calc_objgap,objval_error,validator_status,validator_relobjdiff,iter_count,node_count,conic_count,conic_opt_count,conic_inf_count,relax_status,relax_time,conic_time,mip_time,int_viol,lin_viol,soc_viol,socrot_viol,exp_viol,psd_viol")

# From instance name to file name
function find_instance(name)
    instances = joinpath(dirname(@__FILE__),"..","cbfs")
    try
        match = "$name.*"
        f = readstring(`find $instances -name $match`)
        if length(split(f,"\n")) == 2
            return convert(String,chomp(f))
        end
    end
    return ""
end

function make_smat(svec::Vector{Float64})
    dim = Int(sqrt(1/4+2length(svec)) - 1/2)
    smat = Array{Float64}(dim,dim)
    k = 1

    for i in 1:dim, j in i:dim
        if i == j
            smat[i,j] = svec[k]
        else
            smat[i,j] = (1/sqrt(2))*svec[k]
        end
        k += 1
    end

    return smat
end

function violation_cone(subvec,cone)
    if cone == :Zero
        return :Linear, maximum(abs, subvec)
    elseif cone == :NonNeg
        return :Linear, -minimum(min.(subvec, 0.))
    elseif cone == :NonPos
        return :Linear, maximum(max.(subvec, 0.))
    elseif cone == :SOC
        return :SOC, max(0., vecnorm(subvec[2:end]) - subvec[1])
    elseif cone == :SOCRotated
        # (y,z,x) in RSOC <=> (sqrt2inv*(y+z),sqrt2inv*(-y+z),x) in SOC
        return :SOCRotated, max(0., sqrt(1/2*(subvec[1] - subvec[2])^2 + sum(abs2, subvec[3:end])) - 1/sqrt(2)*(subvec[1] + subvec[2]))
    elseif cone == :ExpPrimal
        return :Exp, max(subvec[2]*exp(subvec[1]/subvec[2]) - subvec[3], 0)
    elseif cone == :Free
        return :Linear, 0.0
    elseif cone == :SDP
        smat = make_smat(subvec)
        return :PSD, -min(eigmin(Symmetric(smat)), 0.0)
    else
        error("Unrecognized cone $cone")
    end
end

function compute_viols(dat, solution)
    c, A, b, con_cones, var_cones, vartypes, sense, objoffset = cbftompb(dat)

    if length(solution) > length(c)
        println("Solution is too long")
    end
    @assert length(solution) >= length(c)
    solution = solution[1:length(c)]
    objval = dot(c,solution)

    y = b - A*solution
    lin_viol = -Inf
    soc_viol = -Inf
    socrot_viol = -Inf
    exp_viol = -Inf
    psd_viol = -Inf
    for (cones,x) in [(var_cones,solution),(con_cones,y)]
        for (cone, idx) in cones
            t, viol = violation_cone(x[idx],cone)
            if t == :Linear
                lin_viol = max(lin_viol,viol)
            elseif t == :SOC
                soc_viol = max(soc_viol,viol)
            elseif t == :SOCRotated
                socrot_viol = max(socrot_viol,viol)
            elseif t == :Exp
                exp_viol = max(exp_viol,viol)
            elseif t == :PSD
                psd_viol = max(psd_viol,viol)
            end
        end
    end

    int_viol = 0.0
    for i in 1:length(vartypes)
        if vartypes[i] == :Int || vartypes[i] == :Bin
            int_viol = max(int_viol, abs(solution[i] - round(solution[i])))
        end
    end

    return objval, lin_viol, soc_viol, socrot_viol, exp_viol, psd_viol, int_viol
end

function validate_with_conic_solver(dat, solution)
    c, A, b, con_cones, var_cones, vartypes, sense, objoffset = cbftompb(dat)
    @assert objoffset == 0.0

    if sense == :Max
        c = -c
    end

    con_cones = copy(con_cones)
    I_eq = Int[]
    J_eq = Int[]
    b_eq = Float64[]
    intcount = 1
    for i in 1:length(vartypes)
        if vartypes[i] == :Int
            push!(I_eq,intcount)
            intcount += 1
            push!(J_eq,i)
            push!(b_eq, round(solution[i]))
        else
            @assert vartypes[i] == :Cont
        end
    end
    con_cones = vcat(con_cones, (:Zero,size(A,1)+1:size(A,1)+intcount-1))
    A = vcat(A, sparse(I_eq,J_eq,ones(length(I_eq)),length(I_eq),length(c)))
    b = vcat(b, b_eq)

    m = MathProgBase.ConicModel(MosekSolver(LOG=0))
    MathProgBase.loadproblem!(m, c, A, b, con_cones, var_cones)
    MathProgBase.optimize!(m)
    status = MathProgBase.status(m)
    if status == :Optimal
        objval = MathProgBase.getobjval(m)
        if sense == :Max
            objval *= -1
        end
    else
        objval = NaN
    end
    return status, objval
end

for (cnt, filename) in enumerate(resultfiles)
    if startswith(basename(filename), ".") || contains(filename,"META")
        continue
    end
    println("$cnt of $(length(resultfiles)): $(basename(filename))")

    solver = split(basename(filename),".")[1]
    # instance name may have '.' in it
    instance = join(split(basename(filename),".")[2:end-1],".")
    instance = string(instance,".cbf")

    timelimit = ""
    status = ""
    objval = Inf
    objbound = -Inf
    solver_time = ""
    total_time = ""
    node_count = ""
    lin_viol = ""
    soc_viol = ""
    socrot_viol = ""
    exp_viol = ""
    psd_viol = ""
    int_viol = ""
    conic_count = ""
    conic_opt_count = ""
    conic_inf_count = ""
    relax_status = ""
    conic_time = ""
    relax_time = ""
    mip_time = ""
    iter_count = ""
    solution = []

    for line in eachline(filename)
        if startswith(line, "#SOLVERNAME#")
            @assert solver == split(line)[2]
        elseif startswith(line, "#INSTANCE#")
            inst = split(line)[2]
            if endswith(inst, ".gz")
                @assert instance == inst[1:end-3]
            else
                @assert instance == inst
            end
        elseif startswith(line, "#TIMELIMIT#")
            timelimit = split(line)[2]
        elseif startswith(line, "#STATUS#")
            status = split(line)[2]
        elseif startswith(line, "#OBJVAL#")
            objval = parse(Float64,split(line)[2])
        elseif startswith(line, "#OBJBOUND#")
            objbound = parse(Float64,split(line)[2])
        elseif startswith(line, "#TIMESOLVER#")
            solver_time = split(line)[2]
        elseif startswith(line, "#TIMEALL#")
            total_time = split(line)[2]
        elseif startswith(line, "#NODECOUNT#")
            node_count = split(line)[2]
        elseif startswith(line, "#SOLUTION#")
            solutionvec = split(line)[2]
            if startswith(solutionvec,'[') && endswith(solutionvec,']') && !startswith(solutionvec, "[]")
                solution = [parse(Float64,x) for x in split(solutionvec[2:end-1],',')]
            end
        elseif startswith(line, " - Iterations           =")
            iter_count = split(line)[4]
        elseif startswith(line, " -- Conic subproblems   =")
            conic_count = split(line)[5]
        elseif startswith(line, " --- Optimal            =")
            conic_opt_count = split(line)[4]
        elseif startswith(line, " --- Infeasible         =")
            conic_inf_count = split(line)[4]
        elseif startswith(line, " - Relaxation status    =")
            relax_status = split(line)[5]
        elseif startswith(line, " -- Solve subproblems   =")
            conic_time = split(line)[5]
        elseif startswith(line, " -- Solve relaxation    =")
            relax_time = split(line)[5]
        elseif startswith(line, " -- Solve MIP models    =") || startswith(line, " -- MIP solver driving  =")
            mip_time = split(line)[6]
        end
    end

    instancefile = find_instance(instance)
    dat = readcbfdata(instancefile)
    sense = string(dat.sense)
    validator_status = ""
    validator_objval = ""
    validator_relobjdiff = ""
    objval_error = ""
    calc_objgap = ""

    # Decide status classification (more can later be excluded if objvals disagree)
    newstatus = ""

    if length(solution) > 0 && all(isfinite,solution)
        objval_sol, lin_viol, soc_viol, socrot_viol, exp_viol, psd_viol, int_viol = compute_viols(dat,solution)

        objval_error = abs(objval_sol - objval)/(abs(objval) + 1e-5)

        if startswith(solver, "BONMIN")
            if status == "Optimal"
                # Have to trust bonmin because cannot query objectivebound
                calc_objgap = 1e-5
            else
                calc_objgap = Inf
            end
        else
            calc_objgap = (objval_sol - objbound) / (abs(objval_sol) + 1e-5)
        end

        if !NOCONIC
            # use conic solver to calculate objval
            validator_status, validator_objval = validate_with_conic_solver(dat, solution)
            validator_relobjdiff = abs(objval_sol - validator_objval)/(abs(objval_sol) + 1e-5)
        end

        if status == "UserLimit" || status == "Suboptimal"
            if parse(Float64, solver_time) > 0.9*parse(Float64, timelimit)
                newstatus = "lim"
            else
                newstatus = "err"
            end
        elseif lin_viol > 1e-6 || soc_viol > 1e-5 || socrot_viol > 1e-5 || exp_viol > 1e-5 || psd_viol > 1e-4 || int_viol > 1e-6
            newstatus = "excl"
        elseif objval_error > 1e-8
            newstatus = "excl"
        elseif !NOCONIC && (validator_status == "Infeasible" || (isfinite(validator_relobjdiff) && validator_relobjdiff > 1e-5))
            newstatus = "excl"
        elseif status == "Optimal"
            if -1e-3 < calc_objgap < 1.05e-5
                newstatus = "conv"
            else
                newstatus = "excl"
            end
        end
    elseif status == "UserLimit"
        newstatus = "lim"
    end

    if status == "Infeasible"
        newstatus = "excl"
    elseif status == "KilledTime" || status == "KilledMemory"
        newstatus = "lim"
    elseif status == "" || status == "Error" || status == "FailedMIP" || status == "UnboundedOA" || status == "FailedOA"
        newstatus = "err"
    end

    if newstatus == ""
        error("new status not handled for $solver $instance with solver status $status")
    end

    println(fd,"$solver,$instance,$(basename(filename)),$sense,$timelimit,$solver_time,$total_time,$status,$newstatus,$objval,$objbound,$calc_objgap,$objval_error,$validator_status,$validator_relobjdiff,$iter_count,$node_count,$conic_count,$conic_opt_count,$conic_inf_count,$relax_status,$relax_time,$conic_time,$mip_time,$int_viol,$lin_viol,$soc_viol,$socrot_viol,$exp_viol,$psd_viol")
end

close(fd)

max_stars_repo_name: chriscoey/PajaritoSupplement

max_issues_repo_name: chriscoey/PajaritoSupplement

max_forks_repo_name: chriscoey/PajaritoSupplement


text:
"""
Optimizer Adam algorithm following

[Kingma et. al., 2014](http://arxiv.org/abs/1412.6980)

Fields:
    - x: the current state
    - t: iteration counter
    - stepsize: s.stepsize(s.t) gives the stepsize to be used in iteration s.t
    - mt: geometrically averaged position
    - vt: geometrically averged componentwise square
    - beta1: geometric factor for averaged position
    - beta2: geometric factor for averaged componentwise square
    - minmas: used in preconditioning step ./(sqrt(vt)+s.minmass
Constructor:
    - construct using `AdamState(x::Array{Float64};stepsize=.001,beta1=.9,beta2=.999,minmass=1e-8)`

"""
type AdamState <: SamplerState
    x::Array{Float64}
    t::Int
    mt::Array{Float64}
    vt::Array{Float64}

    stepsize::Float64
    beta1::Float64
    beta2::Float64
    minmass::Float64

    function AdamState(x::Array{Float64};
        stepsize=.001,beta1=.9,beta2=.999,minmass=1e-8)
        nparams = size(x)
        new(x,0,zeros(nparams),zeros(nparams),
        stepsize,beta1,beta2,minmass)
    end
end

function sample!(s::AdamState, grad::Function)
    gt = grad(s.x)
    s.t += 1
    s.mt = s.beta1*s.mt + (1-s.beta1)*gt
    s.vt = s.beta2*s.vt + (1-s.beta2)*(gt.*gt)
    mt = s.mt/(1.0-s.beta1^s.t) # debiasing
    vt = s.vt/(1.0-s.beta2^s.t) # debiasing
    alphat = s.stepsize./(sqrt(vt)+s.minmass)
    s.x[:] += alphat .* mt
    s
end


text:
amaxmin(x::T, y::T) where {T} = ifelse( abs(x) < abs(y), (y,x), (x,y) )
function amaxmin(x::T, y::T, z::T) where {T}
     y, z = amaxmin(y, z)
     x, z = amaxmin(x, z)
     x, y = amaxmin(x, y)
     return x, y, z
end
function amaxmin(a::T, b::T, c::T, d::T) where {T}
    a, b = amaxmin(a, b)
    c, d = amaxmin(c, d)
    a, c = amaxmin(a, c)
    b, d = amaxmin(b, d)
    b, c = amaxmin(b, c)

    return a, b, c, d
end

"""
    two_sum(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_sum(a::T, b::T) where {T<:Real}
    hi = a + b
    v  = hi - a
    lo = (a - (hi - v)) + (b - v)
    return hi, lo
end

"""
    one_sum(a, b)

Computes `fl(a+b)`.
"""
@inline function one_sum(a::T, b::T) where {T<:Real}
    return a + b
end

"""
   three_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `md = err(a+b+c), lo = err(md)`.
"""
function three_sum(a::T, b::T, c::T) where {T}
    hi, md, lo = amaxmin(a, b, c)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    md, lo = two_hilo_sum(md, lo)
    return hi, md, lo
end

"""
    two_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
"""
@inline function two_sum(a::T, b::T, c::T) where {T}
    hi, md, lo = amaxmin(a, b, c)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    md += lo
    return hi, md
end

"""
    one_sum(a, b, c)

Computes `fl(a+b+c)`
"""
@inline function one_sum(a::T,b::T,c::T) where {T}
    md, lo = two_sum(b, c)
    hi, md = two_sum(a, md)
    md, lo = two_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    md = md + lo
    hi = hi + md
    return hi
end

"""
    four_sum(a, b, c, d)

Computes `hi = fl(a+b+c+d)` and `hm = err(a+b+c+d), ml = err(hm), lo = err(ml)`.
"""
function four_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_sum(a,  b)
    t2, t3 = two_sum(c,  d)
    hi, t4 = two_sum(t0, t2)
    t5, lo = two_sum(t1, t3)
    hm, ml = two_sum(t4, t5)
    ml, lo = two_hilo_sum(ml, lo)
    hm, ml = two_hilo_sum(hm, ml)
    hi, hm = two_hilo_sum(hi,hm)
    return hi, hm, ml, lo
end

"""
    three_sum(a, b, c, d)

Computes `hi = fl(a+b+c+d)` and `md = err(a+b+c+d), lo = err(md)`.
"""
function three_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_sum(a,  b)
    t2, t3 = two_sum(c,  d)
    hi, t4 = two_sum(t0, t2)
    t5, t6 = two_sum(t1, t3)
    md, lo = two_sum(t4, t5)
    lo = lo + t6
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    return hi, md, lo
end

"""
    two_sum(a, b, c, d)

Computes `hi = fl(a+b+c+d)` and `lo = err(a+b+c+d)`.
"""
function two_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_sum(a,  b)
    t2, t3 = two_sum(c,  d)
    hi, t4 = two_sum(t0, t2)
    t5, t6 = two_sum(t1, t3)
    lo, t7 = two_sum(t4, t5)
    t7 = t7 + t6
    lo, t7 = two_hilo_sum(lo, t7)
    hi, lo = two_hilo_sum(hi, lo)
    return hi, lo
end

"""
    one_sum(a, b, c, d)

Computes `fl(a+b+c+d)`.
"""
function one_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_sum(a,  b)
    t2, t3 = two_sum(c,  d)
    hi, t4 = two_sum(t0, t2)
    t5, t6 = two_sum(t1, t3)
    lo, t7 = two_sum(t4, t5)
    t7 = t7 + t6
    lo, t7 = two_hilo_sum(lo, t7)
    hi = hi + lo
    return hi
end

function vec_sum(x0::T, x1::T, x2::T, x3::T) where {T}
    s3 = x3
    s2, e3 = two_sum(x2, s3)
    s1, e2 = two_sum(x1, s2)
    s0, e1 = two_sum(x0, s1)
    return s0,e1,e2,e3
end

function vsum_errbranch(x::NTuple{4,T}) where {T}
    y = zeros(T, 4)
    r = zeros(T, 4)
    e = zeros(T, 4)
    j = 1
    e[1] = x[1]
    for i = 1:2
        r[i], t = two_sum(e[i], x[i+1])
        if t !== zero(T)
            y[j] = r[i]
            e[i+1] = t
            j += 1
        else
            e[i+1] = r[i]
        end
    end
    y[j], y[j+1] = two_sum(e[3], x[4])
    return y
end

function foursum(x1::T, x2::T, x3::T, x4::T) where {T}
    a1, a2 = two_sum(x1, x2)
    b1, b2 = two_sum(x3, x4)
    c1, c2 = two_sum(a1, b1)
    d1, d2 = two_sum(a2, b2)
    e1to4 = vec_sum(c1,c2,d1,d2)
    y = vsum_errbranch(e1to4)
    return (y...,)
end

"""
    five_sum(a, b, c, d, e)

Computes `s = fl(a+b+c+d+e)` and
    `e1 = err(a+b+c+d), e2 = err(e1), e3 = err(e2), e4 = err(e3)`.
"""
function five_sum(v::T, w::T, x::T, y::T, z::T) where {T}
    t0, t4 = two_sum(y, z)
    t0, t3 = two_sum(x, t0)
    t0, t2 = two_sum(w, t0)
    a, t1  = two_sum(v, t0)
    t0, t3 = two_sum(t3, t4)
    t0, t2 = two_sum(t2, t0)
    b, t1  = two_sum(t1, t0)
    t0, t2 = two_sum(t2, t3)
    c, t1  = two_sum(t1, t0)
    d, e   = two_sum(t1, t2)
    return a, b, c, d, e
end

"""
    two_diff(a, b)

Computes `s = fl(a-b)` and `e = err(a-b)`.
"""
@inline function two_diff(a::T, b::T) where {T}
    hi = a - b
    v  = hi - a
    lo = (a - (hi - v)) - (b + v)
    return hi, lo
end

"""
    three_diff(a, b, c)

Computes `s = fl(a-b-c)` and `e1 = err(a-b-c), e2 = err(e1)`.
"""
function three_diff(a::T,b::T,c::T) where {T}
    s, t = two_diff(-b, c)
    x, u = two_sum(a, s)
    y, z = two_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    four_diff(a, b, c, d)

Computes `hi = fl(a-b-c-d)` and `hm = err(a-b-c-d), ml = err(hm), lo = err(ml)`.
"""
function four_diff(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_diff(a ,  b)
    t0, t2 = two_diff(t0,  c)
    hi, t3 = two_diff(t0,  d)
    t0, t1 = two_sum(t1, t2)
    hm, t2 = two_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    two_square(a)

Computes `hi = fl(a*a)` and `lo = fl(err(a*a))`.
"""
@inline function two_square(a::T) where {T}
    hi = a * a
    lo = fma(a, a, -hi)
    hi, lo
end

"""
    two_prod(a, b)

Computes `hi = fl(a*b)` and `lo = fl(err(a*b))`.
"""
@inline function two_prod(a::T, b::T) where {T}
    hi = a * b
    lo = fma(a, b, -hi)
    hi, lo
end

@inline max_min(a,b) = abs(a) < abs(b) ? (b,a) : (a,b)

"""
    three_prod(a, b, c)

Computes `hi = fl(a*b*c)` and `md = err(a*b*c), lo = err(md)`.
"""
function three_prod(a::T, b::T, c::T) where {T}
    a, b = max_min(a,b)
    a, c = max_min(a,c)
    b, c = max_min(b,c)
    abhi, ablo = two_prod(a, b)
    hi, abhiclo = two_prod(abhi, c)
    ablochi, abloclo = two_prod(ablo, c)
    md, lo, tmp  = three_sum(ablochi, abhiclo, abloclo)
    hi, md = two_sum(hi, md)
    return hi, md, lo
end

#=
   three_fma algorithm from
   Sylvie Boldo and Jean-Michel Muller
   Some Functions Computable with a Fused-mac
=#

"""
   three_fma(a, b, c)

Computes `hi = fl(fma(a,b,c))` and `md = fl(err(fma(a,b,c))), lo = fl(err(md))`.
"""
function three_fma(a::T, b::T, c::T) where {T}
    hi = fma(a, b, c)
    hi0, lo0 = two_prod(a, b)
    hi1, lo1 = two_sum(c, lo0)
    hi2, lo2 = two_sum(hi0, hi1)
    y = ((hi2 - hi) + lo2)
    md, lo = two_hilo_sum(y, lo1)
    return hi, md, lo
end

"""
   two_fma(a, b, c)

Computes `hi = fl(fma(a,b,c))` and `lo = fl(err(fma(a,b,c)))`.
"""
function two_fma(a::T, b::T, c::T) where {T}
    hi = fma(a, b, c)
    hi0, lo0 = two_prod(a, b)
    hi1, lo1 = two_sum(c, lo0)
    hi2, lo2 = two_sum(hi0, hi1)
    lo = ((hi2 - hi) + lo2) + lo1
    return hi, lo
end

"""
   two_muladd(a, b, c)

Computes `hi = fl(muladd(a,b,c))` and `lo = fl(err(muladd(a,b,c)))`.
"""
function two_muladd(a::T, b::T, c::T) where {T}
    hi = fma(a, b, c)
    lo = fma(a, b, c-hi)
    return hi, lo
end

# with arguments sorted by magnitude

"""
    two_hilo_sum(a, b)

*unchecked* requirement `|a| ≥ |b|`

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_hilo_sum(a::T, b::T) where {T}
    hi = a + b
    lo = b - (hi - a)
    return hi, lo
end

"""
    two_lohi_sum(a, b)

*unchecked* requirement `|b| ≥ |a|`

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_lohi_sum(a::T, b::T) where {T}
    hi = b + a
    lo = a - (hi - b)
    return hi, lo
end

"""
    two_hilo_diff(a, b)

*unchecked* requirement `|a| ≥ |b|`

Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_hilo_diff(a::T, b::T) where {T}
    hi = a - b
    lo = (a - hi) - b
    hi, lo
end

"""
    two_lohi_diff(a, b)

*unchecked* requirement `|b| ≥ |a|`

Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_lohi_diff(a::T, b::T) where {T}
    hi = b - a
    lo = (b - hi) - a
    hi, lo
end

"""
    two_hilo_sum(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
"""
function two_hilo_sum(a::T, b::T, c::T) where {T}
    lo, t = two_hilo_sum(b, c)
    hi, lo = two_sum(a, lo)
    lo += t
    hi, lo = two_hilo_sum(hi, lo)
    return hi, lo
end

"""
    three_hilo_sum(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`

Computes `x = fl(a+b+c)` and `y = err(a+b+c), z = err(y)`.
"""
function three_hilo_sum(a::T, b::T, c::T) where {T}
    md, lo = two_hilo_sum(b, c)
    hi, md = two_sum(a, md)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    return hi,md,lo
end

"""
    three_lohi_sum(a, b, c)

*unchecked* requirement `|c| ≥ |b| ≥ |a|`

Computes `x = fl(a+b+c)` and `y = err(a+b+c), z = err(y)`.
"""
function three_lohi_sum(a::T,b::T,c::T) where {T}
    md, lo = two_lohi_sum(a, b)
    hi, md = two_sum(c, md)
    md, lo = two_lohi_sum(lo, md)
    hi, md = two_lohi_sum(md, hi)
    return hi,md,lo
end

"""
    three_hilo_diff(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`

Computes `x = fl(a-b-c)` and `y = err(a-b-c), z = err(y)`.
"""
function three_hilo_diff(a::T,b::T,c::T) where {T}
    md, lo = two_hilo_sum(-b, -c)
    hi, md = two_sum(a, md)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    return hi,md,lo
end

"""
    three_lohi_diff(a, b, c)

*unchecked* requirement `|c| ≥ |b| ≥ |a|`

Computes `x = fl(a-b-c)` and `y = err(a-b-c), z = err(y)`.
"""
function three_lohi_diff(c::T,b::T,a::T) where {T}
    three_lohi_sum(c, -b, -a)
end

"""
    four_hilo_sum(a, b, c, d)

*unchecked* requirement `|a| ≥ |b| ≥ |c| ≥ |d|`

Computes `hi = fl(a+b+c+d)` and `hm = err(a+b+c+d), ml = err(hm), lo = err(ml)`.
"""
function four_hilo_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_hilo_sum(a ,  b)
    t0, t2 = two_hilo_sum(t0,  c)
    hi, t3 = two_hilo_sum(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    four_lohi_sum(a, b, c, d)

*unchecked* requirement `|d| ≥ |c| ≥ |b| ≥ |a|`

Computes `hi = fl(a+b+c+d)` and `hm = err(a+b+c+d), ml = err(hm), lo = err(ml)`.
"""
function four_lohi_sum(d::T,c::T,b::T,a::T) where {T}
    t0, t1 = two_hilo_sum(a ,  b)
    t0, t2 = two_hilo_sum(t0,  c)
    hi, t3 = two_hilo_sum(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3)
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    four_hilo_diff(a, b, c, d)

*unchecked* requirement `|a| ≥ |b| ≥ |c| ≥ |d|`

Computes `hi = fl(a-b-c-d)` and `hm = err(a-b-c-d), ml = err(hm), lo = err(ml)`.
"""
function four_hilo_diff(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_hilo_diff(a,  b)
    t0, t2 = two_hilo_diff(t0,  c)
    hi, t3 = two_hilo_diff(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    four_hilo_diff(a, b, c, d)

*unchecked* requirement `|d| ≥ |c| ≥ |b| ≥ |a|`

Computes `hi = fl(a-b-c-d)` and `hm = err(a-b-c-d), ml = err(hm), lo = err(ml)`.
"""
function four_lohi_diff(d::T,c::T,b::T,a::T) where {T}
    t0, t1 = two_hilo_diff(a,  b)
    t0, t2 = two_hilo_diff(t0,  c)
    hi, t3 = two_hilo_diff(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end


text:
# Copyright(c) Rafael Dias
# MIT Licensed
include("propriedades-secoes.jl")

module AnaliseEstruturas
    include("propriedades-secoes.jl")
    mutable struct Elemento
        coordenadaInicial::PropriedadesSecoes.Coordenada
        coordenadaFinal::PropriedadesSecoes.Coordenada
    end

    mutable struct Estrutura
        elementos::Array{Elemento}
    end

    function comprimento(elemento::Elemento)::Float32
        return sqrt(((elemento.coordenadaFinal.x - elemento.coordenadaInicial.x) ^ 2) + ((elemento.coordenadaFinal.y - elemento.coordenadaInicial.y) ^ 2))
    end

    function buildElemento(inicioX::Float32, inicioY::Float32, terminoX::Float32, terminoY::Float32)::Elemento
        return Elemento(PropriedadesSecoes.Coordenada(inicioX, inicioY), PropriedadesSecoes.Coordenada(terminoX, terminoY))
    end

    function metodoDeslocamentos(elemento::Elemento)::Array{Float32}
        @show comprimento(elemento)

        return [0f0]
    end

    function metodoForcas(elemento::Elemento)::Array{Float32}
        @show elemento

        return [0f0]
    end
end

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
@testset "Line" begin
    p1 = SVector(1.0, 1.0, 1.0); p2 = SVector(2.0, 2.0, 0.0);
    l1 = Line(p1, p2)

    @test l1[0] ≈ p1
    @test l1[end/2] ≈ SVector(1.5, 1.5, 0.5)
    @test l1[end] ≈ p2

    @test all(map(isapprox, l1[[0, length(l1)]], [p1, p2]))

    p3 = SVector(1.0, 1.0, 1.0); p4 = SVector(1.0, 1.0, 2.0);
    l2 = Line(p3, p4)

    l3 = Line(p1, p2) # New instance of l1
    @test isequal(l1, l3)
    @test !isequal(l1, l2)
    @test hash(l1) == hash(l3)

    @test distance(l2, SVector(1.0,1.0,1.0)) == 0
    @test distance(SVector(1,1,1.5), l2) == 0
    @test distance(SVector(3,2,3), l2) == sqrt(6)
    @test distance(l2, SVector(3,2,0)) == sqrt(6)

    @test distance(l2, SVector(1.0,1.0)) == 0
    @test distance(SVector(1,1), l2) == 0
    @test distance(SVector(3,2), l2) == sqrt(5)
    @test distance(l2, SVector(3,2)) == sqrt(5)

    l1 = Line(SVector(0.0, 0.0), SVector(1.0, 0.0))
    @test convert2d(l1) == l1
    @test convert3d(l1) == Line(SVector(0.0, 0.0, 0.0), SVector(1.0, 0.0, 0.0))
    @test convert3d(l1, 1.0) == Line(SVector(0.0, 0.0, 1.0), SVector(1.0, 0.0, 1.0))
    @test convert3d(Line(SVector(0.0, 0.0, 1.0), SVector(1.0, 0.0, 1.0)), 2.0) == Line(SVector(0.0, 0.0, 2.0), SVector(1.0, 0.0, 2.0))

    # boundingbox
    @test boundingbox(Line(SVector(1,2,3), SVector(4,5,6))) == BoundingBox(1,2,3,4,5,6)

    @testset "2D Line distances" begin
        # Intersecting lines
        @test distance(Line(SVector(0.0, 0.0), SVector(1.0, 1.0)), Line(SVector(0.0, 1.0), SVector(1.0, 0.0))) === 0.0
        @test distance(Line(SVector(0.3f0, 0.3f0), SVector(0.6f0, 0.6f0)), Line(SVector(0.3f0, 0.6f0), SVector(0.6f0, 0.3f0))) === 0f0

        # Intersecting at vertices
        @test distance(Line(SVector(0.0, 0.0), SVector(1.0, 1.0)), Line(SVector(1.0, 1.0), SVector(1.0, 0.0))) === 0.0
        @test closest_points(Line(SVector(0.0, 0.0), SVector(1.0, 1.0)), Line(SVector(1.0, 1.0), SVector(1.0, 0.0)))[1] === SVector(1.0, 1.0)
        @test closest_points(Line(SVector(0.0, 0.0), SVector(1.0, 1.0)), Line(SVector(1.0, 1.0), SVector(1.0, 0.0)))[2] === SVector(1.0, 1.0)

        @test distance(Line(SVector(0.3f0, 0.3f0), SVector(0.6f0, 0.6f0)), Line(SVector(0.6f0, 0.6f0), SVector(0.6f0, 0.9f0))) === 0f0

        # Parallel lines
        @test distance(Line(SVector(0.0, 0.0), SVector(1.0, 0.0)), Line(SVector(0.0, 0.0), SVector(1.0, 0.0))) === 0.0
        @test distance(Line(SVector(0.3f0, 0.3f0), SVector(0.6f0, 0.6f0)), Line(SVector(0.3f0, 0.3f0), SVector(0.6f0, 0.6f0))) === 0f0
        @test distance(Line(SVector(0.3f0, 0.3f0), SVector(0.6f0, 0.6f0)), Line(SVector(0.4f0, 0.4f0), SVector(0.5f0, 0.5f0))) === 0f0
        @test distance(Line(SVector(0.4f0, 0.4f0), SVector(0.8f0, 0.8f0)), Line(SVector(0.3f0, 0.3f0), SVector(0.6f0, 0.6f0))) === 0f0

        @test distance(Line(SVector(0.0, 0.5), SVector(1.0, 0.5)), Line(SVector(0.0, 0.5), SVector(1.0, 0.5))) === 0.0
        @test distance(Line(SVector(0.0, 0.5), SVector(1.0, 0.5)), Line(SVector(1.0, 0.5), SVector(2.0, 0.5))) === 0.0
        @test distance(Line(SVector(0.0, 0.5), SVector(1.0, 0.5)), Line(SVector(2.0, 0.5), SVector(3.0, 0.5))) === 1.0
        @test distance(Line(SVector(0.0, 0.0), SVector(1.0, 0.0)), Line(SVector(2.0, 0.5), SVector(3.0, 0.5))) ≈ sqrt(1.25)

        # Non-intersecting lines
        @test distance(Line(SVector(0.0, 0.0), SVector(1.0, 0.0)), Line(SVector(1.0, 1.0), SVector(1.0, 2.0))) === 1.0
        @test closest_points(Line(SVector(0.0, 0.0), SVector(1.0, 0.0)), Line(SVector(1.0, 1.0), SVector(1.0, 2.0)))[1] === SVector(1.0, 0.0)
        @test closest_points(Line(SVector(0.0, 0.0), SVector(1.0, 0.0)), Line(SVector(1.0, 1.0), SVector(1.0, 2.0)))[2] === SVector(1.0, 1.0)

    end

    @testset "3D Line distances" begin
        l1 = Line(SVector(0.0, 0.0, 0.0), SVector(1.0, 1.0, 0.0))
        l2 = Line(SVector(0.0, 1.0, 0.5), SVector(1.0, 0.0, 0.5))
        @test distance(l1, l2) ≈ 0.5
        @test closest_points(l1, l2)[1] ≈ SVector(0.5, 0.5, 0.0)
        @test closest_points(l1, l2)[2] ≈ SVector(0.5, 0.5, 0.5)
    end

    # TODO test transformations
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export ConvolutionLayer
function compute_sliding_window_output_shape(input_shape, kernel, pad, stride)
    width, height, channels = input_shape
    width_out = div((width - kernel + 2 * pad), stride) + 1
    height_out = div((height - kernel + 2 * pad), stride) + 1
    width_out, height_out
end

function ConvolutionLayer(name::Symbol, net::Net,
                          input_ensemble::AbstractEnsemble, num_filters::Int,
                          kernel::Int, stride::Int, pad::Int; weight_init=xavier, bias_init::Real=0)
    @assert(ndims(input_ensemble) == 3,
            "Input ensemble to convolution layer must be 3-dimensional")
    width, height, channels = size(input_ensemble)
    # compute the output shape as a function of the input
    width_out, height_out = compute_sliding_window_output_shape(
        size(input_ensemble), kernel, pad, stride)

    # instantiate weights
    weights = weight_init(Float32, kernel*kernel*channels, num_filters)
    ∇weights = zeros(Float32, kernel*kernel*channels, num_filters)

    # instantiate shared bias
    bias = Array(Float32, 1, num_filters)
    fill!(bias, bias_init)
    ∇bias = zeros(Float32, 1, num_filters)

    # allocate neurons
    neurons = Array(WeightedNeuron, width_out, height_out, num_filters)

    # Neurons in the same output channel `o` share filters
    for o = 1:num_filters, h = 1:height_out, w = 1:width_out
        neurons[w, h, o] = WeightedNeuron(view(weights, :, o), view(∇weights, :, o),
                                          view(bias, :, o), view(∇bias, :, o))
    end

    conv = Ensemble(net, name, neurons, [Param(name, :weights, 1.0f0, 1.0f0),
                                         Param(name, :bias, 2.0f0, 0.0f0)])
    add_connections(net, input_ensemble, conv, function (x, y, _)
        in_x = (x-1)*stride - pad
        in_y = (y-1)*stride - pad
        return (in_x+1:in_x+kernel, in_y+1:in_y+kernel, 1:channels)
    end; padding=pad)
    conv
end


text:
using DataFrames
using CSV
using JLD
using Plots
pyplot()

"""Data"""

# Read data from CSV file
df = CSV.read(pwd()*"/data/SNLS80mV.csv", ignoreemptylines=true)
df = select(df, [:V1, :V2])

# Number of samples
num_samples = size(df, 1)

"""Visualization"""

# Check full signals
scatter(1:num_samples, df[:,1], markerstrokewidth=0, markersize=1, size=(1000,400))
xlabel!("samples")
title!("input")
savefig(pwd()*"/data/input.png")

scatter(1:num_samples, df[:,2], markerstrokewidth=0, markersize=1, size=(1000,400))
xlabel!("samples")
title!("output")
savefig(pwd()*"/data/output.png")

# Zooms
zoom = 1:100;

# Check full signals
plot(zoom, df[zoom,1], markerstrokewidth=0, markersize=1, size=(1000,400))
xlabel!("samples")
title!("input")
savefig(pwd()*"/data/input_zoom.png")

plot(zoom, df[zoom,2], markerstrokewidth=0, markersize=1, size=(1000,400))
xlabel!("samples")
title!("output")
savefig(pwd()*"/data/output_zoom.png")

# """Saving"""
#
# #Save the data-frame
# fid = jldopen(pwd()*"/data/SNLS80mV.jld", "w");
# write(fid, "Silverbox", df);
# close(fid);

max_stars_repo_name: wmkouw/NSI-Silverbox

max_issues_repo_name: wmkouw/NSI-Silverbox

max_forks_repo_name: wmkouw/NSI-Silverbox


text:
function op!(Op::ITensor,
             ::OpName"Id",
             ::SiteType"Generic",
             s::Index)
  for n=1:dim(s)
    Op[n,n] = 1.0
  end
end

op!(Op::ITensor,
    ::OpName"F",
    st::SiteType"Generic",
    s::Index) = op!(Op,OpName("Id"),st,s)


text:
"""
FBeam!(FV, U, bdata)
Function evaluation for PTC example.
F(u) = -u'' - lambda sin(u)
"""
function FBeam!(FV, U, bdata)
    D2 = bdata.D2
    lambda = bdata.lambda
    su = lambda * sin.(U)
    FV .= (D2 * U - su)
end

"""
BeamJ!(FP,FV,U,bdata)

Jacobian for the beam problem
F(u) = -u'' - lambda sin(u)
so
F'(u) w  = D2 w - lambda cos(u) w
"""
function BeamJ!(FP, FV, U, bdata)
    D2 = bdata.D2
    lambda = bdata.lambda
    cu = lambda * cos.(U)
    n = length(U)
    zr = zeros(n - 1)
    FP .= D2 - Tridiagonal(zr, cu, zr)
end

"""
BeamtdJ!(FP, FV, U, bdata)

Jacobian evaluation for the time-dependent beam problem.
If F^n(w) = w - u_n + dt F(w) = 0 then
F^n(w)' = I + dt F'(w)
"""
function BeamtdJ!(FP, FV, U, bdata)
    FP .= BeamJ!(FP, FV, U, bdata)
    dt = bdata.dt
    FP .= I + dt * FP
end

"""
FBeamtd!(FV, U, bdata)

Function evaluation for the time-dependent beam problem.
The implicit Euler step for u_t = - F(u) is
u_{n+1} = u_n - dt F(u_{n+1})
so the nonlinear equation is
F^n(w) = w - u_n + dt F(w) = 0
"""
function FBeamtd!(FV, U, bdata)
    un = bdata.UN
    dt = bdata.dt
    FV .= FBeam!(FV, U, bdata)
    dU = U - un
    FV .= dU + dt * FV
    #axpby!(1.0,dU,dt,FV)
end

"""
beaminit(n,dt,lambda=20.0)

Set up the beam problem with n interior grid points.

dt is only needed for the temporal integration examples
"""

function beaminit(n, dt, lambda = 20.0)
    #
    # deltaval is a place to store the current pseudo-time step. I need this
    # for preconditioning. You MUST have this in your pdata because ptcsoli
    # writes to it:  pdata.deltaval[1]=delta
    # So it has to be there in exactly this way.
    #
    deltaval = zeros(1)
    D2 = Lap1d(n)
    dx = 1.0 / (n + 1)
    x = collect(dx:dx:1.0-dx)
    UN = zeros(size(x))
    bdata = (D2 = D2, x = x, dx = dx, dt = dt, lambda = lambda, UN = UN, deltaval = deltaval)
end

"""
Lap1d(n)

returns -d^2/dx^2 on [0,1] zero BC
"""
function Lap1d(n)
    dx = 1 / (n + 1)
    d = 2.0 * ones(n)
    sup = -ones(n - 1)
    slo = -ones(n - 1)
    D2 = Tridiagonal(slo, d, sup)
    D2 = D2 / (dx * dx)
    return D2
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct ParameterMapping{
    VFV<:AbstractVector{<:Real},
    VFI<:AbstractVector{Int},
    VVI<:AbstractVector{Int}
}
    fixed_idxs::VFI
    fixed_values::VFV
    variable_idxs::VVI
end

function ParameterMapping(bounds::HyperRectBounds)
    vol = bounds.vol
    lo_hi_cmp = vol.lo .≈ vol.hi
    fixed_idxs = findall(x -> x == true, lo_hi_cmp)
    fixed_values = vol.lo[fixed_idxs]
    variable_idxs = findall(x -> x == false, lo_hi_cmp)
    ParameterMapping(fixed_idxs, fixed_values, variable_idxs)
end

function map_params!(
    mapped_params::AbstractVector{<:Real},
    parmap::ParameterMapping,
    params::AbstractVector{<:Real}
)
    mapped_params[parmap.fixed_idxs] = parmap.fixed_values
    mapped_params[parmap.variable_idxs] = params
    mapped_params
end

function map_params!(
    mapped_params::VectorOfSimilarVectors{<:Real},
    parmap::ParameterMapping,
    params::VectorOfSimilarVectors{<:Real}
)
    flatview(mapped_params)[parmap.fixed_idxs, :] = parmap.fixed_values
    flatview(mapped_params)[parmap.variable_idxs, :] = flatview(params)
    mapped_params
end

function map_params(parmap::ParameterMapping, params::AbstractVector{<:Real})
    mapped_params = similar(params, size(parmap.fixed_idxs, 1) + size(parmap.variable_idxs, 1))
    map_params!(mapped_params, parmap, params)
end

function map_params(parmap::ParameterMapping, params::VectorOfSimilarVectors{<:Real})
    mapped_params = VectorOfSimilarVectors(similar(flatview(params), size(parmap.fixed_paramsidxs, 1) + size(parmap.variable_idxs, 1), size(flatview(params), 2)))
    map_params!(mapped_params, parmap, params)
end

function invmap_params!(
    params::AbstractVector{<:Real},
    parmap::ParameterMapping,
    mapped_params::AbstractVector{<:Real}
)
    params[:] = view(mapped_params, parmap.variable_idxs)
    params
end

function invmap_params!(
    params::VectorOfSimilarVectors{<:Real},
    parmap::ParameterMapping,
    mapped_params::VectorOfSimilarVectors{<:Real}
)
    flatview(params)[:, :] = view(flatview(mapped_params), parmap.variable_idxs, :)
    params
end

function invmap_params(parmap::ParameterMapping, mapped_params::AbstractVector{<:Real})
    params = similar(mapped_params, size(parmap.variable_idxs, 1))
    invmap_params!(params, parmap, mapped_params)
end

function invmap_params(parmap::ParameterMapping, mapped_params::VectorOfSimilarVectors{<:Real})
    params = VectorOfSimilarVectors(similar(flatview(mapped_params), size(parmap.variable_idxs, 1), size(flatview(mapped_params), 2)))
    invmap_params!(params, parmap, mapped_params)
end

function invmap_params(parmap::ParameterMapping, bounds::HyperRectBounds)
    idxs = parmap.variable_idxs
    vol = HyperRectVolume(bounds.vol.lo[idxs], bounds.vol.hi[idxs])
    bt = bounds.bt[idxs]
    HyperRectBounds(vol, bt)
end

function invmap_params(parmap::ParameterMapping, bounds::NoParamBounds)
    NoParamBounds(size(parmap.variable_idxs, 1))
end


text:
function myfunc()
    A = rand(200, 200, 400)
    maximum(A)
end

myfunc()

using Profile

@profile myfunc()

using ProfileView

Profile.clear()  # in case we have any previous profiling data
@profile myfunc()
ProfileView.view()


text:
# index arithmetic over a string is a bit trickier when using all symbols
# we can use `prevind` and `nextind` to make sure it works properly
@testset "Inter Md 2" begin
    st = raw"""
        ~~~
        this⊙ then ⊙ ⊙ and
        ~~~
        finally ⊙⊙𝛴⊙ and
        ~~~
        escape ∀⊙∀
        ~~~
        done
        """

    inter_md, = explore_md_steps(st)[:inter_md]
    @test inter_md // """
                     ##FDINSERT##
                    finally ⊙⊙𝛴⊙ and
                     ##FDINSERT##
                    done
                    """
end

@testset "MD>HTML" begin
    st = raw"""
        text A1 \newcommand{\com}{blah}text A2 \com and
        ~~~
        escape B1
        ~~~
        \newcommand{\comb}[ 1]{\mathrm{#1}} text C1 $\comb{b}$ text C2
        \newcommand{\comc}[ 2]{part1:#1 and part2:#2} then \comc{AA}{BB}.
        """

    steps = explore_md_steps(st)
    lxdefs, tokens, braces, blocks, lxcoms = steps[:latex]
    b2insert, = steps[:b2insert]
    inter_md, mblocks = steps[:inter_md]
    inter_html, = steps[:inter_html]

    @test isapproxstr(inter_md, """
            text A1 text A2  ##FDINSERT## and
             ##FDINSERT##
             text C1  ##FDINSERT## text C2
             then  ##FDINSERT##.""")

    @test isapproxstr(inter_html, """
text A1 text A2  ##FDINSERT## and  ##FDINSERT##  text C1  ##FDINSERT## text C2  then  ##FDINSERT##.
""")

    hstring = F.convert_inter_html(inter_html, b2insert, lxdefs)
    @test isapproxstr(hstring, raw"""

text A1 text A2 blah and
                    escape B1
                    text C1 \(\mathrm{ b}\) text C2  then part1: AA and part2: BB.
""")
end

@testset "headers" begin
    set_curpath("index.md")
    h = """
        # Title
        and then
        ## Subtitle cool!
        done
        """ |> seval
    @test h // """

Title

and then

Subtitle cool!

done
"""
end


text:
# Subtyping AbstractVector{ElemT} make Base think that RepIterator implements indexing e.g. for copy!
abstract type AbstractRepIterator{N, T, ElemT, PT} end

# For a RepIterator with only one representation, we fallback to the indexing interface, the index being the iterator state
# A representation can overwrite this if it can do something more efficient or if it simply does not support indexing
const SingleRepIterator{N, T, ElemT, RepT} = AbstractRepIterator{N, T, ElemT, Tuple{RepT}}
# With MapRepIterator, we could have RepT<:Rep{N', T'} with N' != N or T' != T, with eachindex, we need to replace everything with N', T'
Base.eachindex(it::SingleRepIterator{<:Any, <:Any, ElemT, RepT}) where {N, T, ElemT, RepT<:Rep{N, T}} = Indices{N, T, similar_type(ElemT, FullDim{N}(), T)}(it.ps[1])
Base.start(it::SingleRepIterator{<:Any, <:Any, ElemT, RepT})     where {N, T, ElemT, RepT<:Rep{N, T}} = start(eachindex(it))::Index{N, T, similar_type(ElemT, FullDim{N}(), T)}
Base.done(it::SingleRepIterator, idx::Index) = done(eachindex(it), idx)
function Base.next(it::SingleRepIterator, idx::Index)
    #mapitem(it, 1, get(it.ps[1], idx)), nextindex(it.ps[1], idx)
    x = get(it.ps[1], idx)
    a = mapitem(it, 1, x)
    b = nextindex(it.ps[1], idx)
    a, b
end

# If there are multiple representations, we need to iterate.
# Builds a SingleRepIterator{ElemT} from p
function checknext(it::AbstractRepIterator, i::Int, state)
    while i <= length(it.ps) && (i == 0 || done(repit(it, i), state))
        i += 1
        if i <= length(it.ps)
            state = start(repit(it, i))
        end
    end
    i > length(it.ps) ? (i, nothing) : (i, state)
end
Base.start(it::AbstractRepIterator) = checknext(it, 0, nothing)
Base.done(it::AbstractRepIterator, state) = state[1] > length(it.ps)

function Base.next(it::AbstractRepIterator, state)
    i, substate = state
    item, newsubstate = next(repit(it, i), substate)
    newstate = checknext(it, i, newsubstate)
    item, newstate
end

# RepIterator
struct RepIterator{N, T, ElemT, PT<:Tuple{Vararg{Rep{N, T}}}} <: AbstractRepIterator{N, T, ElemT, PT}
    ps::PT
    function RepIterator{N, T, ElemT}(ps::PT) where {N, T, ElemT, PT<:Tuple{Vararg{Rep{N, T}}}}
        new{N, T, ElemT, PT}(ps)
    end
end
repit(it::RepIterator{N, T, ElemT}, i::Int) where {N, T, ElemT} = RepIterator{N, T, ElemT}((it.ps[i],))
mapitem(it::RepIterator, i, item) = item

# MapRepIterator
struct MapRepIterator{N, T, ElemT, PT} <: AbstractRepIterator{N, T, ElemT, PT}
    ps::PT
    f::Function
    function MapRepIterator{N, T, ElemT}(ps::PT, f::Function) where {N, T, ElemT, PT<:Tuple}
        new{N, T, ElemT, PT}(ps, f)
    end
end
# N, T are the output dim and coeftype, it can be different that the ones of is.ps[i]
repit(it::MapRepIterator{N, T, ElemT}, i::Int) where {N, T, ElemT} = MapRepIterator{N, T, ElemT}((it.ps[i],), (j, item) -> it.f(i, item)) # j should be 1
mapitem(it::MapRepIterator, i, item) = it.f(i, item)

function Base.first(it::AbstractRepIterator)
    next(it, start(it))[1]
end
FullDim(it::AbstractRepIterator{N}) where N = FullDim{N}()
Base.eltype(it::AbstractRepIterator{N, T, ElemT}) where {N, T, ElemT} = ElemT
function typed_map(f, d::FullDim{N}, ::Type{T}, it::RepIterator{Nin, Tin, ElemT}) where {N, T, Nin, Tin, ElemT}
    MapRepIterator{N, T, similar_type(ElemT, d, T)}(it.ps, f)
end

function RepIterator{N, T}(it::RepIterator) where {N, T}
    typed_map((i,x) -> similar_type(typeof(x), T)(x), FullDim{N}(), T, it)
end

# FIXME the variables need to be defined outside of the local scope of for
#       for Julia to know them inside the """ ... """ of the docstrings
singular = HorV = HorVRep = horvrep = singularlin = plural = plurallin = lenp = isnotemptyp = repexem = listexem = :_

for (isVrep, elt, loop_singular) in [#(true, :SymPoint, :sympoint),
                                     (true, :AbstractPoint, :point),
                                     (true, :Line, :line), (true, :Ray, :ray),
                                     (false, :HyperPlane, :hyperplane), (false, :HalfSpace, :halfspace)]
    global singular = loop_singular
    if isVrep
        vectortype = :vvectortype
        global HorV = :V
        global HorVRep = :VRep
        global horvrep = :vrep
    else
        vectortype = :hvectortype
        global HorV = :H
        global HorVRep = :HRep
        global horvrep = :hrep
    end
    typename = :(AbstractRepIterator{N, T, <:$elt})
    singularstr = string(singular)
    elemtype = Symbol(singularstr * "type")
    donep = Symbol("done" * singularstr)
    startp = Symbol("start" * singularstr)
    nextp = Symbol("next" * singularstr)
    pluralstr = singularstr * "s"
    global plural = Symbol(pluralstr)
    global lenp = Symbol("n" * pluralstr)
    global isnotemptyp = Symbol("has" * pluralstr)
    mapit = Symbol("map" * pluralstr)
    inc = Symbol("incident" * pluralstr)
    incidx = Symbol("incident" * singularstr * "indices")

    @eval begin
        export $plural, $lenp, $isnotemptyp, $startp, $donep, $nextp, $elemtype
        export $inc, $incidx

        """
            $plural($horvrep::$HorVRep)

        Returns an iterator over the $plural of the $HorV-representation `$horvrep`.
        """
        function $plural end

        """
            incident$plural(p::Polyhedron, idx)

        Returns the list of $plural incident to idx for the polyhedron `p`.
        """
        $inc(p::Polyhedron{N, T}, idx) where {N, T} = get(p, IncidentElements{N, T, $elemtype(p)}(p, idx))

        """
            incident$(singular)indices(p::Polyhedron, idx)

        Returns the list of the indices of $plural incident to idx for the polyhedron `p`.
        """
        $incidx(p::Polyhedron{N, T}, idx) where {N, T} = get(p, IncidentIndices{N, T, $elemtype(p)}(p, idx))

        """
            $lenp($horvrep::$HorVRep)

        Returns the number of $plural of the $HorV-representation `$horvrep`.
        """
        $lenp($horvrep::$HorVRep{N, T}) where {N, T} = length(Indices{N, T, $elemtype($horvrep)}($horvrep))

        """
            $isnotemptyp($horvrep::$HorVRep)

        Returns whether the $HorV-representation `$horvrep` has any $singular.
        """
        $isnotemptyp($horvrep::$HorVRep{N, T}) where {N, T} = !isempty(Indices{N, T, $elemtype($horvrep)}($horvrep))

        $elemtype(p::Polyhedron) = $elemtype($horvrep(p))
        if $singularstr == "point"
            $elemtype(p::$HorVRep) = $vectortype(typeof(p))
        else
            $elemtype(p::$HorVRep{N, T}) where {N, T} = $elt{N, T, $vectortype(typeof(p))}
        end

        function $plural(p::$HorVRep{N, T}...) where {N, T}
            ElemT = promote_type($elemtype.(p)...)
            RepIterator{N, T, ElemT}(p)
        end

        function $mapit(f::Function, d::FullDim{N}, ::Type{T}, p::$HorVRep...) where {N, T}
            ElemT = promote_type(similar_type.($elemtype.(p), d, T)...)
            MapRepIterator{N, T, ElemT}(p, f)
        end

        Base.length(it::$typename) where {N, T} = sum($lenp, it.ps)
        Base.isempty(it::$typename) where {N, T} = !any($isnotemptyp.(it.ps))
    end
end

# Combines an element type with its linear version.
# e.g. combines rays with the lines by splitting lines in two points.
struct AllRepIterator{N, T, ElemT, LinElemT, PT}
    itlin::RepIterator{N, T, LinElemT, PT}
    it::RepIterator{N, T, ElemT, PT}
end

Base.eltype(it::AllRepIterator{N, T, ElemT}) where {N, T, ElemT} = ElemT
Base.length(it::AllRepIterator) = 2length(it.itlin) + length(it.it)
Base.isempty(it::AllRepIterator) = isempty(it.itlin) && isempty(it.it)

function checknext(it::AllRepIterator, i, state)
    while i <= 3 && ((i <= 2 && done(it.itlin, state)) || (i == 3 && done(it.it, state)))
        i += 1
        if i <= 2
            state = start(it.itlin)
        elseif i == 3
            state = start(it.it)
        end # Otherwise we leave state as it is to be type stable
    end
    i, state
end

splitlin(h::HyperPlane, i) = (i == 1 ? HalfSpace(h.a, h.β) : HalfSpace(-h.a, -h.β))
#splitlin(s::SymPoint, i) = (i == 1 ? coord(s) : -coord(s))
splitlin(l::Line, i) = (i == 1 ? Ray(coord(l)) : Ray(-coord(l)))

Base.start(it::AllRepIterator) = checknext(it, 1, start(it.itlin))
Base.done(::AllRepIterator, state) = state[1] > 3
function Base.next(it::AllRepIterator, istate)
    i, state = istate
    if i <= 2
        @assert i >= 1
        itemlin, newstate = next(it.itlin, state)
        item = splitlin(itemlin, i)
    else
        @assert i == 3
        item, newstate = next(it.it, state)
    end
    newistate = checknext(it, i, newstate)
    (item, newistate)
end

for (isVrep, loop_singularlin,
     loop_singular, loop_repexem,
     loop_listexem) in [#(true, :sympoint, :point, "convexhull(SymPoint([1, 0]), [0, 1])", "[1, 0], [-1, 0], [0, 1]"),
                   (true, :line, :ray, "Line([1, 0]) + Ray([0, 1])",
                    "Ray([1, 0]), Ray([-1, 0]), Ray([0, 1])"),
                   (false, :hyperplane, :halfspace,
                    "HyperPlane([1, 0], 1) ∩ HalfSpace([0, 1], 1)",
                    "HalfSpace([1, 0]), HalfSpace([-1, 0]), HalfSpace([0, 1])")]
    global singularlin = loop_singularlin
    global singular = loop_singular
    global repexem = loop_repexem
    global listexem = loop_listexem
    if isVrep
        global HorV = :V
        global HorVRep = :VRep
        global horvrep = :vrep
    else
        global HorV = :H
        global HorVRep = :HRep
        global horvrep = :hrep
    end
    pluralstrlin = string(singularlin) * "s"
    global plurallin = Symbol(pluralstrlin)
    lenplin = Symbol("n" * pluralstrlin)
    isnotemptyplin = Symbol("has" * pluralstrlin)
    pluralstr = string(singular) * "s"
    global plural = Symbol(pluralstr)
    global lenp = Symbol("n" * pluralstr)
    global isnotemptyp = Symbol("has" * pluralstr)
    allpluralstr = "all" * pluralstr
    allplural = Symbol(allpluralstr)
    alllenp = Symbol("n" * allpluralstr)
    allisnotemptyp = Symbol("has" * allpluralstr)
    @eval begin
        export $allplural, $alllenp, $allisnotemptyp

        """
            all$plural($horvrep::$HorVRep)

        Returns an iterator over the $plural and $plurallin in the $HorV-representation `$horvrep` splitting $plurallin in two $plural.

        ### Examples

        ```julia
        $horvrep = $repexem
        collect(all$plural($horvrep)) # Returns [$listexem]
        ```
        """
        $allplural(p::$HorVRep...) = AllRepIterator($plurallin(p...), $plural(p...))

        """
            nall$plural($horvrep::$HorVRep)

        Returns the number of $plural plus twice the number of $plurallin in the $HorV-representation `$horvrep`, i.e. `length(all$plural($horvrep))`
        """
        $alllenp($horvrep::$HorVRep) = 2 * $lenplin($horvrep) + $lenp($horvrep)

        """
            hasall$plural($horvrep::$HorVRep)

        Returns whether the $HorV-representation `$horvrep` contains any $singular or $singularlin.
        """
        $allisnotemptyp($horvrep::$HorVRep) = $isnotemptyplin($horvrep) || $isnotemptyp($horvrep)
    end
end

const ElemIt{ElemT} = Union{AllRepIterator{<:Any, <:Any, ElemT}, AbstractRepIterator{<:Any, <:Any, ElemT}, AbstractVector{ElemT}}
const HyperPlaneIt{N, T} = ElemIt{<:HyperPlane{N, T}}
const HalfSpaceIt{N, T} = ElemIt{<:HalfSpace{N, T}}
const HIt{N, T} = Union{HyperPlaneIt{N, T}, HalfSpaceIt{N, T}}

#const SymPointIt{N, T} = ElemIt{<:SymPoint{N, T}}
const PointIt{N, T} = ElemIt{<:AbstractPoint{N, T}}
const PIt{N, T} = PointIt{N, T}
#const PIt{N, T} = Union{SymPointIt{N, T}, PointIt{N, T}}
const LineIt{N, T} = ElemIt{<:Line{N, T}}
const RayIt{N, T} = ElemIt{<:Ray{N, T}}
const RIt{N, T} = Union{LineIt{N, T}, RayIt{N, T}}
const VIt{N, T} = Union{PIt{N, T}, RIt{N, T}}

const It{N, T} = Union{HIt{N, T}, VIt{N, T}}

function fillvits(::FullDim{N}, points::ElemIt{AT}, lines::ElemIt{Line{N, T, AT}}=Line{N, T, AT}[], rays::ElemIt{Ray{N, T, AT}}=Ray{N, T, AT}[]) where {N, T, AT<:AbstractPoint{N, T}}
    if isempty(points) && !(isempty(lines) && isempty(rays))
        vconsistencyerror()
    end
    return points, lines, rays
end
function fillvits(::FullDim{N}, lines::ElemIt{Line{N, T, AT}}, rays::ElemIt{Ray{N, T, AT}}=Ray{N, T, AT}[]) where {N, T, AT}
    if isempty(lines) && isempty(rays)
        ps = AT[]
    else
        ps = [origin(AT, FullDim{N}())]
    end
    ps, lines, rays
end

hreps(p::HRep{N, T}...) where {N, T} = hyperplanes(p...), halfspaces(p...)
hreps(p::HAffineSpace{N, T}...) where {N, T} = tuple(hyperplanes(p...))

hmap(f, d::FullDim, ::Type{T}, p::HRep...) where T = maphyperplanes(f, d, T, p...), maphalfspaces(f, d, T, p...)
hmap(f, d::FullDim, ::Type{T}, p::HAffineSpace...) where T = tuple(maphyperplanes(f, d, T, p...))

hconvert(RepT::Type{<:HRep{N, T}}, p::HRep{N, T}) where {N, T} = constructpolyhedron(RepT, (p,), hreps(p)...)
hconvert(RepT::Type{<:HRep{N, T}}, p::HRep{N})    where {N, T} = constructpolyhedron(RepT, (p,), RepIterator{N, T}.(hreps(p))...)

vreps(p...) = preps(p...)..., rreps(p...)...
preps(p::VRep...) = tuple(points(p...))
preps(p::VCone...) = tuple()
rreps(p::VRep...) = lines(p...), rays(p...)
rreps(p::VLinearSpace...) = tuple(lines(p...))
rreps(p::VPolytope...) = tuple()

vmap(f, d::FullDim, ::Type{T}, p::VRep...) where T = pmap(f, d, T, p...)..., rmap(f, d, T, p...)...
pmap(f, d::FullDim, ::Type{T}, p::VRep...) where T = tuple(mappoints(f, d, T, p...))
pmap(f, d::FullDim, ::Type, p::VCone...) = tuple()
rmap(f, d::FullDim, ::Type{T}, p::VRep...) where T = maplines(f, d, T, p...), maprays(f, d, T, p...)
rmap(f, d::FullDim, ::Type{T}, p::VLinearSpace...) where T = tuple(maplines(f, d, T, p...))
rmap(f, d::FullDim, ::Type, p::VPolytope...) = tuple()

vconvert(RepT::Type{<:VRep{N, T}}, p::VRep{N, T}) where {N, T} = constructpolyhedron(RepT, (p,), vreps(p)...)
vconvert(RepT::Type{<:VRep{N, T}}, p::VRep{N})    where {N, T} = constructpolyhedron(RepT, (p,), RepIterator{N, T}.(vreps(p))...)


text:
############################################################
## joAbstractSAparallelLinearOperator - outer constructors
############################################################

#export

# outer constructors
export joSAdistributedLinOp
"""
    julia> joSAdistributedLinOp(A,psin; [,fclean] [,rclean])

Create a linear operator working on 2D SharedArray in multi-vector (over 2nd dimension) mode.

# Signature

    function joSAdistributedLinOp(A::joAbstractLinearOperator{DDT,RDT},psin::joPAsetup,
        fclean::Bool=false,rclean::Bool=false) where {DDT<:Number,RDT<:Number,INT<:Integer}

# Arguments

- `A`: joAbstractLinearOperator type
- `psin`: parallel setup from jpPAsetup
- `fclean`: close SharedArray after forward operation
- `rclean`: close SharedArray after forward operation in transpose/adjoint mode

# Examples

- `ps=joPAsetup((12,30))`: define parallel setup for operator with 12 columns and multi-vector with 30 columns
- `joSAdistributedLinOp(A,ps)`: operator that will apply A to distributed multi-vector with 30 columns

# Notes

- nvc of joPAsetup must be >= then # of workers in the WorkerPool
- DT in joPAsetup will be overwritten by DDT of the operator

"""
function joSAdistributedLinOp(A::joAbstractLinearOperator{DDT,RDT},psin::joPAsetup,
        fclean::Bool=false,rclean::Bool=false) where {DDT<:Number,RDT<:Number,INT<:Integer}

    psin.chunks[1]==1 || throw(joSAdistributedLinearOperatorException("joSAdistributedLinearOperator: invalid joPAsetup - must not be distributed in 1st dimension"))
    A.n==psin.dims[1]  || throw(joSAdistributedLinearOperatorException("joSAdistributedLinearOperator: invalid joPAsetup - 1st dimension off joPAsetup does not match size(A,2)"))

    wpool=WorkerPool(psin.procs)
    nvc=psin.dims[2]
    parts=joPAsetup_etc.parts(psin,2)
    idst=joPAsetup(wpool,(A.n,nvc),2,parts=parts,DT=DDT)
    odst=joPAsetup(wpool,(A.m,nvc),2,parts=parts,DT=RDT)
    return joSAdistributedLinearOperator{DDT,RDT,2}("joSAdistributedLinearOperator($(A.name))",A.m,A.n,nvc,
               v1->A*v1,
               v2->transpose(A)*v2,
               v3->adjoint(A)*v3,
               v4->conj(A)*v4,
               @joNF, @joNF, @joNF, @joNF,
               idst,odst,fclean,rclean)
end
"""
    julia> joSAdistributedLinOp(A,nvc; [parts] [,fclean] [,rclean])

Create a linear operator working on 2D SharedArray in multi-vector (over 2nd dimension) mode.

# Signature

    function joSAdistributedLinOp(wpool::WorkerPool,A::joAbstractLinearOperator{DDT,RDT},nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),nvc);
        fclean::Bool=false,rclean::Bool=false) where {DDT<:Number,RDT<:Number,INT<:Integer}
    joSAdistributedLinOp(A::joAbstractLinearOperator{ADDT,ARDT},nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),nvc);
        fclean::Bool=false,rclean::Bool=false) where {ADDT<:Number,ARDT<:Number,INT<:Integer}

# Arguments

- `A`: joAbstractLinearOperator type
- `nvc`: number of columns in multi-vector
- `parts`: custom partitioning of 2nd diemnsion
- `fclean`: close SharedArray after forward operation
- `rclean`: close SharedArray after forward operation in transpose/adjoint mode

# Examples

- `joSAdistributedLinOp(A,30)`: operator that will apply A to distributed multivector with 30 columns

# Notes

- nvc must be >= then # of workers in the WorkerPool

"""
function joSAdistributedLinOp(wpool::WorkerPool,A::joAbstractLinearOperator{DDT,RDT},nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),nvc);
        fclean::Bool=false,rclean::Bool=false) where {DDT<:Number,RDT<:Number,INT<:Integer}

    length(parts)==nworkers(wpool) || throw(joSAdistributedLinearOperatorException("joSAdistributedLinearOperator: lenght(parts) does not much nworkers()"))
    sum(parts)==nvc || throw(joSAdistributedLinearOperatorException("joSAdistributedLinearOperator: sum(parts) does not much nvc"))

    idst=joPAsetup(wpool,(A.n,nvc),2,parts=parts,DT=DDT)
    odst=joPAsetup(wpool,(A.m,nvc),2,parts=parts,DT=RDT)
    return joSAdistributedLinearOperator{DDT,RDT,2}("joSAdistributedLinearOperator($(A.name))",A.m,A.n,nvc,
               v1->A*v1,
               v2->transpose(A)*v2,
               v3->adjoint(A)*v3,
               v4->conj(A)*v4,
               @joNF, @joNF, @joNF, @joNF,
               idst,odst,fclean,rclean)
end
joSAdistributedLinOp(A::joAbstractLinearOperator{ADDT,ARDT},nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),nvc);
        fclean::Bool=false,rclean::Bool=false) where {ADDT<:Number,ARDT<:Number,INT<:Integer} =
        joSAdistributedLinOp(WorkerPool(workers()),A,nvc,fclean=fclean,rclean=rclean)

############################################################
## joAddSolver{Any,All}

# exported in joAbstractOperator/constructors.jl

# joAddSolver - outer constructor for adding solver to operator
#joAddSolverAny(A::joSAdistributedLinOp{DDT,RDT},slvr::Function) where {DDT,RDT} =
#end
#joAddSolverAll(A::joSAdistributedLinOp{DDT,RDT},
#    slvr::Function,slvr_T::Function,slvr_A::Function,slvr_C::Function) where {DDT,RDT} =
#end


text:
@testset "wrapping_integers_cmp.jl" begin
    @testset begin
        @test WrappingInt32(3) != WrappingInt32(1)
        N_REPS = 4096
        for i in 1:N_REPS
            n = rand(typemin(UInt32):typemax(UInt32))
            diff = rand(UInt16(1):typemax(UInt16))
            m = n + diff
            @test WrappingInt32(m) != WrappingInt32(n)
        end
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
include("data.jl")
# splitticks(x) = ( [ e[1] for e in x], [ e[2] for e in x])

mordat = [ (alter=a, jahr=j, monat=monat, geschlecht=g,
            # interpolate for monat:
            N=population(popdata; alter=a, jahr=j, monat=monat, geschlecht=g),
            # interpolate with deaths for monat:
            # N=population(popdata, deathdat; alter=a, jahr=j, monat=monat, geschlecht=g),
            D=deaths(deathdat; alter=a, jahr=j, monat=monat, geschlecht=g))
           # TODO: integrate as population_deaths(popdata; alter=a, geschlecht=g, jahr=j, deaths=(monat=monat,) )
           for a in DeutschlandDaten.Altersgruppen_Monate
               for g in ["Männlich", "Weiblich"]
                   for j in 2000:2021
                       for monat in 1:12
                           ] |> DataFrame;

filter!(x-> x.D !==missing && x.D>0, mordat)
@transform! mordat :date = Date.(:jahr, :monat, 1)
# @transform! mordat :Ndiff = :N .- :Ns

@transform! mordat :group = collect(zip(:geschlecht, :alter))

# us = unstack(d[:,Symbol[:row,:col,:value]], :col, :value)

# @transform! groupby(mordat, :date) :stackedP = cumsum(:P_AG)

# P_AG: Wahrscheinlichkeit P(A=alter, G=geschlecht | T=date)
# PD: Wahrscheinlichkeit P(D = 1 | A=alter, G=geschlecht, T=date)
dataprep!(mordat)

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
using DataDeps

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

dd = DataDep("cubed_sphere_32_grid",
    "Conformal cubed sphere grid with 32×32 grid points on each face",
    "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/cubed_sphere_grids/cubed_sphere_32_grid.jld2"
)

DataDeps.register(dd)


text:
@enum PFibTile begin
    α₁=1
    α₂=2
    β=3
end
L = [α₁, α₂, β]
Square = Tuple{PFibTile, PFibTile, PFibTile, PFibTile}
transitions = Dict{PFibTile, Vector{PFibTile}}([α₁ => [α₂,β], α₂ => [α₁,β], β => [α₁]])
H₊ = [[α₁, α₂], [α₁, β], [α₂, β], [β, α₁]]
H₀ = [[ℓ, ℓ] for ℓ in L]
S₊₊ = [ [transitions[h[1]][end], transitions[h[2]][1], h[1], h[2] ] for h in H₊ ]
S₊₀ = [[α₂, β, α₁, α₁], [α₁, β, α₂, α₂]]
S₊ = S₊₊ ∪ S₊₀
S₀ = [ [ℓ, ℓ, κ, κ] for κ in L for ℓ in transitions[κ]]
P = [ [ℓ, κ] for κ in L for ℓ in transitions[κ]]

function complete_up(h)
    return [square for square in S₀ ∪ S₊ if square[[1,2]] == h]
end

Path = Vector{PFibTile}
BBisection = Tuple{Path, Path}
Bisection = Vector{BBisection}

S = [ ([ℓ₁, κ₁], [ℓ₂, κ₂]) for κ₁ in L for ℓ₁ in transitions[κ₁] for κ₂ in L for ℓ₂ in transitions[κ₂]
    if (([κ₁, κ₂] ∈ H₊) || (ℓ₁ < ℓ₂ && κ₁ == κ₂))]
function promote_bisection(B :: Bisection, n)
    m = length(B[1][1])
    for m = m:n-1
        new_B = Bisection()
        for O in B
            h = [O[1][end], O[2][end]]
            for square in complete_up(h)
                push!(new_B, ([O[1];square[3]], [O[2];square[4]]))
            end
        end
        B = new_B
    end
    return B
end

function commutes_from(O :: BBisection, n)
    m = length(O[1])
    for i=n:m-1
        if [O[1][i], O[2][i], O[1][i+1], O[2][i+1]] ∉ S₊ ∪ S₀
            return false
        end
    end
    return true
end

function Base.:∘(B₁ :: Bisection, B₂ :: Bisection)
    if isempty(B₁) || isempty(B₂)
        return Bisection()
    end
    m = length(B₁[1][1])
    n = length(B₂[1][1])
    N = max(m,n)
    B₁ = promote_bisection(B₁, N+2)
    B₂ = promote_bisection(B₂, N+2)
    B = Bisection()
    for O₁ in B₁
        for O₂ in B₂
            if O₁[2] == O₂[1]
                push!(B, (O₁[1], O₂[2]))
            end
        end
    end
    return B
end

part = []
T = filter(!isempty, [[O₁] ∘ [O₂] ∘ [O₃] for O₁ in S for O₂ in S for O₃ in S])
promote_bisection([S[1]] ∘ [S[5]],5)


text:
"""
    TicToc -- timing measurement

Low-overhead time interval measurement via minimally invasive macros.

"""
module TicToc

using Printf

export @tic, @toc, tictoc

# disable to reduce overhead
const tictoc_track_memory = true

if tictoc_track_memory

mutable struct TimingInfo
    ncalls::Int
    time::UInt64
    allocd::Int64
    gctime::UInt64
    curr::UInt64
    mem::Base.GC_Num
end
TimingInfo() = TimingInfo(0, 0, 0, 0, 0,
                          Base.GC_Num(0, 0, 0, 0, 0, 0, 0,
                                      0, 0, 0, 0, 0, 0, 0))

else # !tictoc_track_memory

mutable struct TimingInfo
    ncalls::Int
    time::UInt64
    curr::UInt64
end
TimingInfo() = TimingInfo(0, 0, 0)

end # if tictoc_track_memory

const timing_infos = TimingInfo[]
const timing_info_names = Symbol[]
const atexit_function_registered = Ref(false)

# `@tic` helper
function _tic(nm)
    exti = Symbol("tictoc__", nm)
    global timing_info_names
    if exti in timing_info_names
        throw(ArgumentError("$(nm) already used in @tic"))
    end
    push!(timing_info_names, exti)
    @static if tictoc_track_memory
        quote
            global $(exti)
            $(exti).curr = time_ns()
            $(exti).mem = Base.gc_num()
        end
    else
        quote
            global $(exti)
            $(exti).curr = time_ns()
        end
    end
end

"""
    @tic nm

Indicate the start of the interval `nm`.
"""
macro tic(args...)
    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in @tic"))
    end
    ex = args[1]
    if !isa(ex, Symbol)
        throw(ArgumentError("need a name argument to @tic"))
    end
    return _tic(ex)
end

# `@toc` helper
function _toc(nm)
    exti = Symbol("tictoc__", nm)
    @static if tictoc_track_memory
        quote
            global $(exti)
            $(exti).time += time_ns() - $(exti).curr
            $(exti).ncalls += 1
            local diff = Base.GC_Diff(Base.gc_num(), $(exti).mem)
            $(exti).allocd += diff.allocd
            $(exti).gctime += diff.total_time
        end
    else
        quote
            global $(exti)
            $(exti).time += time_ns() - $(exti).curr
            $(exti).ncalls += 1
        end
    end
end

"""
    @toc nm

Indicate the end of the interval `nm`.
"""
macro toc(args...)
    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in @toc"))
    end
    ex = args[1]
    if !isa(ex, Symbol)
        throw(ArgumentError("need a name argument to @toc"))
    end
    return _toc(ex)
end

"""
    print_timing_info()

`atexit()` function, writes all information about every interval to
`stdout`.
"""
function print_timing_info()
    println("TicToc timing information")
    @static if tictoc_track_memory
        println("name,ncalls,tottime(ns),allocbytes,gctime")
    else
        println("name,ncalls,tottime(ns)")
    end
    for i = 1:length(timing_info_names)
        @static if tictoc_track_memory
            s = @sprintf("%s,%d,%d,%d,%d", timing_info_names[i],
                         timing_infos[i].ncalls, timing_infos[i].time,
                         timing_infos[i].allocd, timing_infos[i].gctime)
        else
            s = @sprintf("%s,%d,%d", timing_info_names[i],
                         timing_infos[i].ncalls, timing_infos[i].time)
        end
        println(s)
    end
end

"""
    tictoc()

Call at program start (only once!) to set up the globals used by the
macros and to register the at-exit callback.
"""
function tictoc()
    global timing_info_names
    for nm in timing_info_names
        exti = Symbol(nm)
        isdefined(@__MODULE__, exti) && continue
        expr = quote
            const $exti = $TimingInfo()
        end
        eval(Expr(:toplevel, expr))
        push!(timing_infos, getfield(@__MODULE__, exti))
    end
    if parse(Int, get(ENV, "TICTOC_PRINT_RESULTS", "0")) == 1
        if !atexit_function_registered[]
            atexit(print_timing_info)
            atexit_function_registered[] = true
        end
    end
    return length(timing_info_names)
end

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Catalyst_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Catalyst")
JLLWrappers.@generate_main_file("Catalyst", UUID("513fda4f-8b2e-59b2-a185-097ac7c26d11"))
end  # module Catalyst_jll


text:
export variables, nvariables, exponents, degree, isconstant, powers, constantmonomial, mapexponents

"""
    monomialtype(p::AbstractPolynomialLike)

Return the type of the monomials of `p`.

    termtype(::Type{PT}) where PT<:AbstractPolynomialLike

Returns the type of the monomials of a polynomial of type `PT`.
"""
monomialtype(::Union{M, Type{M}}) where M<:AbstractMonomial = M
monomialtype(::Union{PT, Type{PT}}) where PT <: APL = monomialtype(termtype(PT))
monomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}) where PT <: APL = monomialtype(PT)

"""
    variables(p::AbstractPolynomialLike)

Returns the tuple of the variables of `p` in decreasing order. It could contain variables of zero degree, see the example section.

### Examples

Calling `variables(x^2*y)` should return `(x, y)` and calling `variables(x)` should return `(x,)`.
Note that the variables of `m` does not necessarily have nonzero degree.
For instance, `variables([x^2*y, y*z][1])` is usually `(x, y, z)` since the two monomials have been promoted to a common type.
"""
function variables end

"""
    nvariables(p::AbstractPolynomialLike)

Returns the number of variables in `p`, i.e. `length(variables(p))`. It could be more than the number of variables with nonzero degree (see the Examples section of [`variables`](@ref)).

### Examples

Calling `nvariables(x^2*y)` should return at least 2 and calling `nvariables(x)` should return at least 1.
"""
nvariables(::Union{AbstractVariable, Type{<:AbstractVariable}}) = 1

"""
    exponents(t::AbstractTermLike)

Returns the exponent of the variables in the monomial of the term `t`.

### Examples

Calling `exponents(x^2*y)` should return `(2, 1)`.
"""
exponents(t::AbstractTerm) = exponents(monomial(t))
exponents(v::AbstractVariable) = (1,)

"""
    degree(t::AbstractTermLike)

Returns the *total degree* of the monomial of the term `t`, i.e. `sum(exponents(t))`.

    degree(t::AbstractTermLike, v::AbstractVariable)

Returns the exponent of the variable `v` in the monomial of the term `t`.

### Examples

Calling `degree(x^2*y)` should return 3 which is ``2 + 1``.
Calling `degree(x^2*y, x)` should return 2 and calling `degree(x^2*y, y)` should return 1.

"""
degree(t::AbstractTermLike) = sum(exponents(t))

degree(t::AbstractTermLike, var::AbstractVariable) = degree(monomial(t), var)
degree(v::AbstractVariable, var::AbstractVariable) = (v == var ? 1 : 0)
#_deg(v::AbstractVariable) = 0
#_deg(v::AbstractVariable, power, powers...) = v == power[1] ? power[2] : _deg(v, powers...)
#degree(m::AbstractMonomial, v::AbstractVariable) = _deg(v, powers(t)...)

function degree(m::AbstractMonomial, v::AbstractVariable)
    i = findfirst(isequal(v), variables(m))
    if i === nothing || iszero(i)
        0
    else
        exponents(m)[i]
    end
end

"""
    isconstant(t::AbstractTermLike)

Returns whether the monomial of `t` is constant.
"""
isconstant(t::AbstractTermLike) = all(iszero.(exponents(t)))
isconstant(v::AbstractVariable) = false

"""
    powers(t::AbstractTermLike)

Returns an tuple of the powers of the monomial of `t`.

### Examples

Calling `powers(3x^4*y) should return `((x, 4), (y, 1))`.
"""
powers(t::AbstractTermLike) = tuplezip(variables(t), exponents(t))

"""
    constantmonomial(p::AbstractPolynomialLike)

Returns a constant monomial of the monomial type of `p` with the same variables as `p`.

    constantmonomial(::Type{PT}) where {PT<:AbstractPolynomialLike}

Returns a constant monomial of the monomial type of a polynomial of type `PT`.
"""
function constantmonomial end

"""
    mapexponents(f, m1::AbstractMonomialLike, m2::AbstractMonomialLike)

If ``m_1 = \\prod x^{\\alpha_i}`` and ``m_2 = \\prod x^{\\beta_i}`` then it returns the monomial ``m = \\prod x^{f(\\alpha_i, \\beta_i)}``.

### Examples

The multiplication `m1 * m2` is equivalent to `mapexponents(+, m1, m2)`, the unsafe division `_div(m1, m2)` is equivalent to `mapexponents(-, m1, m2)`, `gcd(m1, m2)` is equivalent to `mapexponents(min, m1, m2)`, `lcm(m1, m2)` is equivalent to `mapexponents(max, m1, m2)`.
"""
mapexponents(f, m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(f, monomial(m1), monomial(m2))

Base.one(::Type{TT}) where {TT<:AbstractMonomialLike} = constantmonomial(TT)
Base.one(t::AbstractMonomialLike) = constantmonomial(t)
# See https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/issues/82
# By default, Base do oneunit(v::VT) = VT(one(v)).
# This tries to convert a monomial to a variable which does not work.
# The issue here is there is no way to represent the multiplicative identity
# using the variable type. The best we can do is return a monomial even
# if it does not exactly match the definition of oneunit.
Base.oneunit(v::AbstractVariable) = one(v)
Base.oneunit(VT::Type{<: AbstractVariable}) = one(VT)


text:
using .RationalFunctionFields: modular_reduction, rational_reconstruction, CRT,
                                CRT_reconstruction

@testset "reduction + reconstruction tests" begin
    modulo = 2^31 - 1

    FF = Sing.N_ZpField(modulo)

    xs  = [ Sing.QQ(rand(1:10000), rand(1:10000)) for _ in 1:10 ]
    xrs = modular_reduction(xs, FF)
    ys  = rational_reconstruction(xrs, BigInt(modulo))

    @test xs == ys

end

@testset "complex reduction tests" begin
    R, (x, y) = AA.PolynomialRing(Sing.QQ, ["x", "y"])

    modulo = 2^31 - 1
    FF = AA.GF(2^31 - 1)

    fs = [
        R(0),
        R(1),
        7x^2*y + 5x + 3y,
        (x + y)^5
    ]

    frs = modular_reduction(fs, FF)

    fr = frs[3]
    @test AA.base_ring(fr) == FF

    frs = rational_reconstruction(frs, BigInt(modulo))

    @test fs == frs
end

@testset "CRT tests" begin

    moduli = [11]
    xs = [3, 4, 5, 10]

    for x in xs
        rems = [mod(x, mi) for mi in moduli]
        @test CRT(moduli, rems) == x
    end

    moduli = Primes.primes(5, 20)
    m = prod(moduli)

    xs = [ rand(1:(m-1)) for _ in 1:100 ]

    for x in xs
        rems = [mod(x, mi) for mi in moduli]

        @test CRT(moduli, rems) == x
    end

    moduli = Primes.primes(2^5, 2^6)
    m = prod(moduli)

    xs = [ rand(1:(m-1)) for _ in 1:100 ]

    for x in xs
        rems = [mod(x, mi) for mi in moduli]

        @test CRT(moduli, rems) == x
    end

end

@testset "CRT basis reconstruction" begin
    moduli = [7, 5]

    FF1 = Sing.N_ZpField(moduli[1])
    FF2 = Sing.N_ZpField(moduli[2])

    R1, (x, y) = AA.PolynomialRing(FF1, ["x1", "y1"])
    R11, (a, b) = AA.PolynomialRing(AA.FractionField(R1), ["a1", "b1"])

    fs1 = [
        4*a^2*b^2 + (x + 6y)//(x + y)*a + (3x + 2x*y + 1*y)//(y)*b + 5y^2,
        1a + 3b + 2a*b + 2x//(1y)
    ]

    R2, (x, y) = AA.PolynomialRing(FF2, ["x2", "y2"])
    R22, (a, b) = AA.PolynomialRing(AA.FractionField(R2), ["a2", "b2"])

    fs2 = [
        2*a^2*b^2 + (1x + 3y)//(x + y)*a + (3x + 1x*y + 4*y)//(y)*b + 5y^2,
        1a + 3b + 1a*b + 4x//(y)
    ]

    basis = CRT_reconstruction([fs1, fs2], moduli)

    yoverx = AA.parent(basis[1])
    basepolyrings = AA.base_ring(yoverx)
    ground = AA.ZZ

    x, y = AA.gens(AA.base_ring(basepolyrings))
    a, b = AA.gens(yoverx)

    # to check when number is in the denominator
    ans = [
           32*a^2*b^2 + (x + 13*y)//(x + y)*a + (16*x*y + 3*x + 29*y)//(y)*b + 5*y^2,
           1a + 3b + 16a*b + 9x//(y)
    ]

    @test basis == ans

end


text:
@testset "Inference" begin
    i = AnalyticVI()
    @test length(i) == 1
end


text:
using ClosedLoopReachability
using Test

@time @testset "Toy model (black-box network)" begin include("black_box_toy_model.jl") end


text:
#!/usr/bin/env julia

"""
This simple program is a loop that prints out each key pressed. It terminates
on Ctrl+D.
"""

using KeyboardDriven

function main()
    keyboarddriven()
    while true
        k = getkeypress()
        println(k)
        if k == ctrl | CharKey('d')
            break
        end
    end
end

main()


text:
using Documenter, FlowFarm

makedocs(;
    modules=[FlowFarm],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/byuflowlab/FlowFarm.jl/blob/{commit}{path}#L{line}",
    sitename="FlowFarm.jl",
    authors="PJ Stanley <stanley_andrewpj@yahoo.com>",
    assets=String[],
)


text:
using Surrogates

#1D
n = 10
lb = 0.0
ub = 10.0
x = sample(n,lb,ub,SobolSample())
f = x -> 2*x
y = f.(x)
my_varfid = VariableFidelitySurrogate(x,y,lb,ub)
val = my_varfid(3.0)
add_point!(my_varfid,3.0,6.0)
val = my_varfid(3.0)

my_varfid_change_struct = VariableFidelitySurrogate(x,y,lb,ub, num_high_fidel = 2,
                                                    low_fid_structure = InverseDistanceStructure(p = 1.0),
                                                    high_fid_structure = RadialBasisStructure(radial_function = linearRadial(), scale_factor=1.0,sparse = false))
#ND
n = 10
lb = [0.0,0.0]
ub = [5.0,5.0]
x = sample(n,ub,lb,SobolSample())
f = x -> x[1]*x[2]
y = f.(x)
my_varfidND = VariableFidelitySurrogate(x,y,lb,ub)
val = my_varfidND((2.0,2.0))
add_point!(my_varfidND,(3.0,3.0),9.0)
my_varfidND_change_struct = VariableFidelitySurrogate(x,y,lb,ub, num_high_fidel = 2,
                                                    low_fid_structure = InverseDistanceStructure(p = 1.0),
                                                    high_fid_structure = RadialBasisStructure(radial_function = linearRadial(), scale_factor=1.0,sparse = false))


text:
"""
    servedocs_callback!(docwatcher, filepath, path2makejl, literate, foldername)

Custom callback used in [`servedocs`](@ref) triggered when the file corresponding to `filepath`
is changed. If that file is `docs/make.jl`, the callback will check whether any new files have
been added in the `docs/src` folder and add them to the watched files, it will also remove any
file that may have been deleted or renamed.
Otherwise, if the modified file is in `docs/src` or is `docs/make.jl`, a pass of Documenter is
triggered to regenerate the documents, subsequently the LiveServer will render the produced pages
in `docs/build`.
`foldername` can be set other than "docs" if needed.
`buildfoldername` can be set other than "build" if needed.
"""
function servedocs_callback!(dw::SimpleWatcher, fp::AbstractString, makejl::AbstractString,
                             literate::String="", skip_dirs::Vector{String}=String[],
                             foldername::String="docs", buildfoldername="build")
    # ignore things happening in build (generated files etc)
    startswith(fp, joinpath(foldername, buildfoldername)) && return nothing
    if !isempty(skip_dirs)
        for dir in skip_dirs
            startswith(fp, dir) && return nothing
        end
    end
    # if the file that was changed is the `make.jl` file, assume that maybe new files are # referenced and so refresh the vector of watched files as a result.
    if fp == makejl
        # it's easier to start from scratch (takes negligible time)
        empty!(dw.watchedfiles)
        scan_docs!(dw, literate, foldername)
    end
    fext = splitext(fp)[2]
    P1 = fext ∈ (".md", ".jl", ".css")
    # the second condition is for CSS files, we want to track it but not the output
    # if we track the output then there's an infinite loop being triggered (see docstring)
    if P1
        Main.include(makejl)
        file_changed_callback(fp)
    end
    return nothing
end

"""
    scan_docs!(dw::SimpleWatcher, literate="", foldername="docs")

Scans the `docs/` folder in order to recover the path to all files that have to be watched and add
those files to `dw.watchedfiles`. The function returns the path to `docs/make.jl`. A list of
folders and file paths can also be given for files that should be watched in addition to the
content of `docs/src`. `foldername` can be changed if it's different than docs.
"""
function scan_docs!(dw::SimpleWatcher, literate::String="", foldername::String="docs")
    src = joinpath(foldername, "src")
    if !(isdir(foldername) && isdir(src))
        @error "I didn't find a $foldername/ or $foldername/src/ folder."
    end
    makejl = joinpath(foldername, "make.jl")
    push!(dw.watchedfiles, WatchedFile(makejl))
    if isdir(foldername)
        # add all files in `docs/src` to watched files
        for (root, _, files) ∈ walkdir(joinpath(foldername, "src")), file ∈ files
            push!(dw.watchedfiles, WatchedFile(joinpath(root, file)))
        end
    end
    if !isempty(literate)
        isdir(literate) || @error "I didn't find the provided literate folder $literate."
        for (root, _, files) ∈ walkdir(literate), file ∈ files
            push!(dw.watchedfiles, WatchedFile(joinpath(root, file)))
        end
    end

    # When using literate.jl, we should only watch the source file otherwise we would double
    # trigger: first when the script.jl is modified then again when the script.md is created
    # which would cause an infinite loop if both `script.jl` and `script.md` are watched.
    # So here we remove from the watchlist all files.md that have a files.jl with the same path.
    remove = Int[]
    if isempty(literate)
        # assumption is that the scripts are in `docs/src/...` and that the generated markdown
        # goes in exactly the same spot so for instance:
        # docs
        # └── src
        #     ├── index.jl
        #     └── index.md
        for wf ∈ dw.watchedfiles
            spath = splitext(wf.path)
            spath[2] == ".jl" || continue
            k = findfirst(e -> splitext(e.path) == (spath[1], ".md"), dw.watchedfiles)
            k === nothing || push!(remove, k)
        end
    else
        # assumption is that the scripts are in `literate/` and that the generated markdown goes
        # in `docs/src` with the same relative paths so for instance:
        # docs
        # ├── lit
        # │   └── index.jl
        # └── src
        #     └── index.md
        for (root, _, files) ∈ walkdir(literate), file ∈ files
            spath = splitext(joinpath(root, file))
            spath[2] == ".jl" || continue
            path = replace(spath[1], Regex("^$literate") => joinpath(foldername, "src"))
            k = findfirst(e -> splitext(e.path) == (path, ".md"), dw.watchedfiles)
            k === nothing || push!(remove, k)
        end
    end
    deleteat!(dw.watchedfiles, remove)
    return makejl
end

"""
    servedocs(; verbose=false, literate="", doc_env=false, foldername="docs")

Can be used when developing a package to run the `docs/make.jl` file from Documenter.jl and
then serve the `docs/build` folder with LiveServer.jl. This function assumes you are in the
directory `[MyPackage].jl` with a subfolder `docs`.

* `verbose=false` is a boolean switch to make the server print information about file changes and
connections.
* `literate=""` is the path to the folder containing the literate scripts, if left empty, it will be
assumed that they are in `docs/src`.
* `doc_env=false` is a boolean switch to make the server start by activating the doc environment or not (i.e. the `Project.toml` in `docs/`).
* `skip_dir=""` is a subpath of `docs/` where modifications should not trigger the generation of the docs, this is useful for instance if you're using Weave and Weave generates some files in `docs/src/examples` in which case you should give `skip_dir=joinpath("docs","src","examples")`.
* `skip_dirs=[]` same as `skip_dir`  but for a vector of such dirs. Takes precedence over `skip_dir`.
* `foldername="docs"` specify the name of the content folder if different than "docs".
* `buildfoldername="build"` specify the name of the build folder if different than "build".
* `host="127.0.0.1"` where the server will start.
* `port` is an integer between 8000 (default) and 9000.
* `launch_browser=false` specifies whether to launch the ambient browser at the localhost URL or not.
"""
function servedocs(; verbose::Bool=false, literate::String="",
                     doc_env::Bool=false, skip_dir::String="",
                     skip_dirs::Vector{String}=String[],
                     foldername::String="docs",
                     buildfoldername::String="build",
                     host::String="127.0.0.1", port::Int=8000,
                     launch_browser::Bool = false)
    # Custom file watcher: it's the standard `SimpleWatcher` but with a custom callback.
    docwatcher = SimpleWatcher()

    if isempty(skip_dirs) && !isempty(skip_dir)
        skip_dirs = [skip_dir]
    end

    set_callback!(docwatcher,
                  fp->servedocs_callback!(docwatcher, fp, makejl, literate, skip_dirs, foldername, buildfoldername))

    # Retrieve files to watch
    makejl = scan_docs!(docwatcher, literate, foldername)

    if doc_env
        Pkg.activate("$foldername/Project.toml")
    end
    # trigger a first pass of Documenter (& possibly Literate)
    Main.include(abspath(makejl))

    # note the `docs/build` exists here given that if we're here it means the documenter
    # pass did not error and therefore that a docs/build has been generated.
    serve(docwatcher, host=host, port=port, dir=joinpath(foldername, buildfoldername), verbose=verbose, launch_browser=launch_browser)
    if doc_env
        Pkg.activate()
    end
    return nothing
end

#
# Miscellaneous utils
#

"""
    setverbose(b)

Set the verbosity of LiveServer to either `true` (showing messages upon events) or `false` (default).
"""
setverbose(b::Bool) = (VERBOSE.x = b)

"""
    example()

Simple function to copy an example website folder to the current working directory that can be
watched by the LiveServer to get an idea of how things work.

### Example

```julia
LiveServer.example()
cd("example")
serve()
```
"""
function example(; basedir="")
    isempty(basedir) && (basedir = pwd())
    cp(joinpath(dirname(dirname(pathof(LiveServer))), "example"), joinpath(basedir, "example"))
end

#
# Generate example repo for servedocs
#

const INDEX_MD = raw"""
    # Test
    A link to the [other page](/man/pg1/)
    """
const PG1_JL  = raw"""
    # # Test literate
    # We can include some code like so:
    f(x) = x^5
    f(5)
    """
const MAKE_JL = raw"""
    using Documenter, Literate
    src = joinpath(@__DIR__, "src")
    lit = joinpath(@__DIR__, "literate")
    for (root, _, files) ∈ walkdir(lit), file ∈ files
        splitext(file)[2] == ".jl" || continue
        ipath = joinpath(root, file)
        opath = splitdir(replace(ipath, lit=>src))[1]
        Literate.markdown(ipath, opath)
    end
    makedocs(
        sitename = "testlit",
        pages = ["Home" => "index.md",
                 "Other page" => "man/pg1.md"]
        )
    """

"""
servedocs_literate_example(dir="servedocs_literate_example")

Generates a folder with the right structure for servedocs+literate example.
You can then `cd` to that folder and use servedocs:

```
julia> using LiveServer
julia> LiveServer.servedocs_literate_example()
julia> cd("servedocs_literate_example")
julia> servedocs(literate=joinpath("docs","literate"))
```
"""
function servedocs_literate_example(dirname="servedocs_literate_example")
    isdir(dirname) && rm(dirname, recursive=true)
    mkdir(dirname)
    # folder structure
    src  = joinpath(dirname, "src")
    mkdir(src)
    write(joinpath(src, "$dirname.jl"), "module $dirname\n foo()=1\n end")
    docs = joinpath(dirname, "docs")
    mkdir(docs)
    src = joinpath(docs, "src")
    lit = joinpath(docs, "literate")
    man = joinpath(lit, "man")
    mkdir(src)
    mkdir(lit)
    mkdir(man)
    write(joinpath(src, "index.md"), INDEX_MD)
    write(joinpath(man, "pg1.jl"), PG1_JL)
    write(joinpath(docs, "make.jl"), MAKE_JL)
    return
end


text:
#deps = ["ImagineWorker", "ImagineInterface"]

#for dep in deps
#    if Base.find_in_path(dep) == nothing
#        Pkg.clone("https://github.com/HolyLab/$dep.git")
#        Pkg.build(dep)
#    end
#end