text:
using AutoHashEquals

part1(input) = play(input, part = 1)
part2(input) = play(input, part = 2)

@auto_hash_equals mutable struct State
    mana::Int
    hero_hp::Int
    boss_hp::Int
    shield_timer::Int
    poison_timer::Int
    recharge_timer::Int
end

function play(input; part)
    boss_hp, boss_damage = parse.(Int, last.(split.(eachline(input), ": ")))
    hero_hp = 50
    hero_mana = 500

    queue = [(0, State(hero_mana, hero_hp, boss_hp, 0, 0, 0))]
    explored_states = Set{State}()
    while !isempty(queue)
        _, i = findmin(first, queue)
        spent_mana, state = queue[i]
        deleteat!(queue, i)

        if part == 2
            state.hero_hp -= 1
            state.hero_hp < 0 && continue
        end
        apply_effects!(state)
        state.boss_hp <= 0 && return spent_mana
        for spell! in (magic_missile!, drain!, shield!, poison!, recharge!)
            next_state = deepcopy(state)
            success, cost = spell!(next_state)
            success || continue
            next_state.mana -= cost
            next_state.boss_hp <= 0 && return spent_mana + cost
            apply_effects!(next_state)
            next_state.boss_hp <= 0 && return spent_mana + cost
            boss_hit!(next_state, boss_damage)
            next_state.hero_hp <= 0 && continue
            next_state in explored_states && continue
            push!(queue, (spent_mana + cost, next_state))
            push!(explored_states, next_state)
        end
    end
    @assert false
end

function apply_effects!(state)
    if state.poison_timer > 0
        state.boss_hp -= 3
    end
    if state.recharge_timer > 0
        state.mana += 101
    end
    state.shield_timer -= 1
    state.poison_timer -= 1
    state.recharge_timer -= 1
end

function magic_missile!(state)
    state.mana < 53 && return false, 53
    state.boss_hp -= 4
    return true, 53
end

function drain!(state)
    state.mana < 73 && return false, 73
    state.boss_hp -= 2
    state.hero_hp += 2
    return true, 73
end

function shield!(state)
    state.mana < 113 && return false, 113
    state.shield_timer > 0 && return false, 113
    state.shield_timer = 6
    return true, 113
end

function poison!(state)
    state.mana < 173 && return false, 173
    state.poison_timer > 0 && return false, 173
    state.poison_timer = 6
    return true, 173
end

function recharge!(state)
    state.mana < 229 && return false, 229
    state.recharge_timer > 0 && return false, 229
    state.recharge_timer = 5
    return true, 229
end

function boss_hit!(state, damage)
    armor = state.shield_timer >= 0 ? 7 : 0
    state.hero_hp -= max(1, damage - armor)
end


text:
function Base.:(*)(A::AbstractMatrix, B::OneHotLike{<:Any, L}) where L
  _isonehot(B) || return invoke(*, Tuple{AbstractMatrix, AbstractMatrix}, A, B)
  size(A, 2) == L || throw(DimensionMismatch("Matrix column must correspond with OneHot size: $(size(A, 2)) != $L"))
  return A[:, onecold(B)]
end

function Base.:(*)(A::AbstractMatrix, B::OneHotLike{<:Any, L, 1}) where L
  _isonehot(B) || return invoke(*, Tuple{AbstractMatrix, AbstractMatrix}, A, B)
  size(A, 2) == L || throw(DimensionMismatch("Matrix column must correspond with OneHot size: $(size(A, 2)) != $L"))
  return NNlib.gather(A, _indices(B))
end

function Base.:(*)(A::AbstractMatrix, B::Adjoint{Bool, <:OneHotMatrix})
  B_dim = length(_indices(parent(B)))
  size(A, 2) == B_dim || throw(DimensionMismatch("Matrix column must correspond with OneHot size: $(size(A, 2)) != $B_dim"))
  return NNlib.scatter(+, A, _indices(parent(B)), dstsize=(size(A,1), size(B,2)))
end

for wrapper in [:Adjoint, :Transpose]
  @eval begin
    function Base.:*(A::$wrapper{<:Any, <:AbstractMatrix{T}}, b::OneHotVector{<:Any, L}) where {L, T}
      size(A, 2) == L ||
          throw(DimensionMismatch("Matrix column must correspond with OneHot size: $(size(A, 2)) != $L"))

      return A[:, onecold(b)]
    end

    function Base.:*(A::$wrapper{<:Number, <:AbstractVector{T}}, b::OneHotVector{<:Any, L}) where {L, T}
      size(A, 2) == L ||
          throw(DimensionMismatch("Matrix column must correspond with OneHot size: $(size(A, 2)) != $L"))

      return A[onecold(b)]
    end
  end
end

ext: jl


text:
using OpenQuantumBase, Test

H = DenseHamiltonian([(s) -> 1 - s, (s) -> s], [σx, σz])
dH(s) = (-real(σx) + real(σz))
coupling = ConstantCouplings(["Z"])

@test_logs (:warn, "The projection method only works with real Hamitonians. Convert the complex Hamiltonian to real one.") project_to_lowlevel(H, [0.0, 0.5, 1.0], coupling, dH)

t_obj = project_to_lowlevel(H, [0.0, 0.5, 1.0], coupling, dH)

@test t_obj.ev ≈ [[-1.0, 1.0], [-0.707107, 0.707107], [-1.0, 1.0]] atol = 1e-6
@test t_obj.dθ ≈ [[0.5], [1.0], [0.5]]
@test get_dθ(t_obj) ≈ -[0.5, 1.0, 0.5]

@test t_obj.op ./ 2 ./ π ≈ [
    [[0 -1.0; -1.0 0]],
    [[-0.707107 -0.707107; -0.707107 0.707107]],
    [[-1.0 0.0; 0.0 1.0]],
] atol = 1e-6

# An exact solvable example
H = DenseHamiltonian([(s)->cos(π*s/2), (s)->sin(π*s/2)], [σx, σz])
dH(s) = π*(-sin(π*s/2)*σx+cos(π*s/2)*σz)/2
coupling = ConstantCouplings(["Z"])

t_obj = project_to_lowlevel(H, 0:0.01:1, coupling, dH)
@test all((x)->isapprox(x, [-1, 1]), t_obj.ev)
@test all((x)->isapprox(x[1], π/4), t_obj.dθ)
t_obj = project_to_lowlevel(H, 0:0.01:1, coupling, dH, direction=:backward)
@test all((x)->isapprox(x, [-1, 1]), t_obj.ev)
@test all((x)->isapprox(x[1], π/4), t_obj.dθ)

#TODO: the following tests randomly fail on Julia 1.7, need to find a permanent fix

#= H = SparseHamiltonian(
    [(s) -> 1 - s, (s) -> s],
    [-standard_driver(2, sp=true) / 2, (spσz ⊗ spσi - 0.1spσz ⊗ spσz) / 2],
)
dH(s) = standard_driver(2, sp=true) / 2 + (spσz ⊗ spσi - 0.1spσz ⊗ spσz) / 2
coupling = ConstantCouplings(["ZI", "IZ"])

t_obj = project_to_lowlevel(H, [0.0, 0.5, 1.0], coupling, dH)

@test t_obj.ev ≈ [
    [-1.0, 0.0],
    [-0.6044361719689455, -0.10443617196894575],
    [-0.55, -0.45],
] atol = 1e-6 =#


text:
module TestCompositeDataFrames

using Base.Test
using DataArrays, DataFrames
using DataFramesMeta

df = CompositeDataFrame(A = [1, 2, 3], B = [2, 1, 2])
x = [2, 1, 0]

@test  df.A == df[:A]
@test  size(df[1:2,:]) == (2,2)
@test  size(df[[1, 2],:]) == (2,2)
@test  size(df[[1, 2], [:A]]) == (2,1)
@test  size(df[:, [:A]]) == (3,1)
@test  size(df[:, 1:2]) == (3,2)
@test  size(df[:, [1]]) == (3,1)
@test  size(df[df.A .< 3, :]) == (2,2)
@test  size([df df]) == (3,4)

@test  @with(df, :A + 1)   ==  df[:A] + 1
@test  @with(df, :A + :B)  ==  df[:A] + df[:B]
@test  @with(df, :A + x)   ==  df[:A] + x
x = @with df begin
    res = 0.0
    for i in 1:length(:A)
        res += :A[i] * :B[i]
    end
    res
end
@test  x == sum(df[:A] .* df[:B])
@test  @with(df, df[:A .> 1, ^([:B, :A])]) == df[df[:A] .> 1, [:B, :A]]
@test  @with(df, DataFrame(a = :A * 2, b = :A + :B)) == DataFrame(a = df[:A] * 2, b = df[:A] + df[:B])

@test  @ix(df, :A .> 1)           == df[df[:A] .> 1,:]
@test  @ix(df, :B .> 1)           == df[df[:B] .> 1,:]
@test  @ix(df, :A .> x)           == df[df[:A] .> x,:]
@test  @ix(df, :B .> x)           == df[df[:B] .> x,:]
@test  @ix(df, :A .> :B)          == df[df[:A] .> df[:B],:]
@test  @ix(df, :A .> 1, [:B, :A]) == df[df[:A] .> 1, [:B, :A]]

@test  @where(df, :A .> 1)         == df[df[:A] .> 1,:]
@test  @where(df, :B .> 1)         == df[df[:B] .> 1,:]
@test  @where(df, :A .> x)         == df[df[:A] .> x,:]
@test  @where(df, :B .> x)         == df[df[:B] .> x,:]
@test  @where(df, :A .> :B)        == df[df[:A] .> df[:B],:]

df2 = @transform(df, C = :A + 1)
@test df2.C == df.A + 1
@test df2[:C] == df2[:A] + 1

df3 = @select(df, :B, C = :A + 1, :A)
@test df3.C == df.A + 1

end # module


text:
using Cronbach
using Test

using LinearAlgebra: diag
using Statistics: cov

@testset "Cronbach.jl" begin
    include("test_tauequivalent.jl")
end


text:
using PDMats
import Statistics: mean, var
using LinearAlgebra

"""
   `calc_periodogram(times, y_obs, σ_ys; options...)`
Computes a periodogram for times, observations and uncertainties.

Inputs:
- times: obseration times
- y_obs: observed values at each time
- σ_ys:  measurement uncertainty at each time (may be either a Vector or positive definite matrix)
Optional parameters:

- period_min (2.0)
- period_max (4*span of times)
- num_periods (4000)

Output (as NamedTuple):
- period_grid: orbital periods searched
- periodogram: χ² for best fit given each orbital period
- period_best_fit: period searched resulting to minimum χ²
- coeff_best_fit: best-fit coefficients for period_best_fit
- phase_best_fit: orbital phase for period_best_fit at time = 0
- predict: predictions of best-fit model at each time
- rms: root mean square of predictions (of best-fit model) minus observations
"""
function calc_periodogram  end

function calc_periodogram(t::Vector, y_obs::Vector, σ_ys::Vector; kwargs...)
 	calc_periodogram(t, y_obs, PDiagMat(σ_ys.^2); kwargs...)
end

function calc_periodogram(t::Vector, y_obs::Vector, covar_mat::MT;
	    		period_min::Real = 2.0, period_max::Real = 4*(maximum(t)-minimum(t)), num_periods::Integer = 4000) where {
						T3<:Real, MT<:AbstractPDMat{T3} }
	period_grid =  1.0 ./ range(1.0/period_max, stop=1.0/period_min, length=num_periods)
	periodogram = map(p->-0.5*calc_χ²_general_linear_least_squares(calc_design_matrix_circ(p,t),covar_mat,y_obs),period_grid)
	period_fit = period_grid[argmax(periodogram)]
	design_matrix_fit = calc_design_matrix_circ(period_fit,t)
	coeff_fit = fit_general_linear_least_squares(design_matrix_fit,covar_mat,y_obs)
	phase_fit = atan(coeff_fit[1],coeff_fit[2])
	pred = design_matrix_fit * coeff_fit
	rms = sqrt(mean((y_obs.-pred).^2))
	return (;period_grid=period_grid, periodogram=periodogram, period_best_fit = period_fit, coeff_best_fit=coeff_fit, phase_best_fit=phase_fit, predict=pred, rms=rms )
end

# Internal routines called by calc_periodogram[]
function calc_design_matrix_circ!(result::AM, period::Real, times::AV) where { R1<:Real, AM<:AbstractMatrix{R1}, AV<:AbstractVector{R1} }
	n = length(times)
	@assert size(result,1) == n
	@assert 2<= size(result,2) <=3
	for i in 1:n
		( result[i,1], result[i,2] ) = sincos(2π/period .* times[i])
	end
	if size(result,2) == 3
		result[:,3] .= one(eltype(result))
	end
	return result
end

function calc_design_matrix_circ(period::Real, times::AV) where { R1<:Real, AV<:AbstractVector{R1} }
	n = length(times)
	dm = Array{promote_type(typeof(period),eltype(times))}(undef,n,3)
	calc_design_matrix_circ!(dm,period,times)
	return dm
end

function fit_general_linear_least_squares( design_mat::ADM, covar_mat::APD, obs::AA ) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	Xt_inv_covar_X = Xt_invA_X(covar_mat,design_mat)
	X_inv_covar_y =  design_mat' * (covar_mat \ obs)
	AB_hat = Xt_inv_covar_X \ X_inv_covar_y                            # standard GLS
end

function predict_general_linear_least_squares( design_mat::ADM, covar_mat::APD, obs::AA ) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	param = fit_general_linear_least_squares(design_mat,covar_mat,obs)
	design_mat * param
end

function calc_χ²_general_linear_least_squares( design_mat::ADM, covar_mat::APD, obs::AA ) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	pred = predict_general_linear_least_squares(design_mat,covar_mat,obs)
	invquad(covar_mat,obs-pred)
end


text:
mutable struct cfgStruct
    meterperSeg::Float64
    chn2SegBegIdx::UInt16
    reflectorFactor::Float64
    attenutionFactor::Float64
    scanrate::UInt16
    calSamples::UInt16
    SegNumber::UInt16
    function cfgStruct()
        new(0)
    end #cfgStruct
end

mutable struct  CfgDataStruct
        SpatialResolution::Float64 ;
        Refractivity::Float64 ;
        OffsetM::Float64 ;
        ATT::Float64 ;
        ScanRate::Float64 ;
        EnergyCnt ::UInt
        LcLimitHigh::Float64 ;
        LcLimitLow::Float64 ;
        StartCh1::UInt
        LenCh1::UInt
        StartCh2::UInt
        LenCh2::UInt
        function CfgDataStruct()
            new(0)
        end #CfgDataStruct

end         #mutable struct  CfgDataStruct

mutable struct RNG
    AlrmProcStartID::Int
    AlrmProcEndtID1::Int
    startID1::Int
    endID1::Int

    chnid::Int
    ProtectSpaceWidth::Int
    function RNG()
        new(0)
    end
end

mutable struct spaceCumSumData
    m1::Matrix{UInt32 }  ##振动数据
    curRow::Int
    maxTimeCnt::Int
    # rectSum::Array{ UInt64}
    dt::String

    ################################
    startID::Int #保存数据的起点
    endID::Int  #  保存数据的终点

    AlrmProcStartID::Int
    AlrmProcEndtID::Int

    defaultSpaceCnt::Int
    defaultTimeCnt::Int
    CfgInfoObj

    rectSum::Array{ UInt64}
    WindAlarmSegs
    SpaceSegs
    isInit
    function  spaceCumSumData()
        new()
    end
end # spaceCumSumData

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using NPZ
using Plots
using Plots.PlotMeasures
using ArgParse
using StatsBase
using HypothesisTests
using LsqFit
using Compose
using LaTeXStrings

function plot_minimal_energies(file::String)
        D = npzread(file)

        p = plot(size = (300, 200))

        x = 0.
        try
            x = D["annealing_times"]
            xlabel!(L"\textrm{annelaing} \ \textrm{time} \ \mu s")
        catch
            x = D["betas"]
            xlabel!(L"$\beta_{MH}$")
            plot!(xscale = (:log))
        end

        y =  [D["ground"] for _ in 1:length(x)]
        α = D["alpha"]
        Z = D["minimum_estimated"]

        plot!(p, x, y, label = "true ground", markershape = :auto, color = "green", legend=false)

        plot!(p, x, D["minimum_from_data"], markershape = :auto, color = "red", label = "minimum from D-Wave samples")
        plot!(p, x, Z, label = "minimum estimated", markershape = :auto,  color = "black")

        plot!(p, x, D["bootstrap_q95"], color = "gray", label = "95 percentile of estimated min.", line = (:dash))
        ylabel!("energy")

        str = "minimal_energies"
        file1 = replace(file, ".npz" => str*".pdf")
        savefig(p, file1)

end

function plot_bootstrad_std(file::String, l::Int)
        D = npzread(file)

        p = plot(size = (400, 250))
        x = 0.

        y = D["bootstrap_std"]
        y1 = D["squared_error"]

        try
            x = D["annealing_times"]
            xlabel!(L"\textrm{annelaing} \ \textrm{time} \ \mu s")
        catch
            x = D["betas"]
            xlabel!(L"$ \beta_{MH} $")
            #plot!(p, ylims = (0, 1.1*maximum(y1[1:10])))
            plot!(p, xaxis = (:log))
            #vline!([x[l]], style = :dot, linewidth = 2., color = "green", label = "model limit")
        end

        α = D["alpha"]

        plot!(p, x, y, label = "bootsrap", line = (:green, 1.5), marker = (:dot, :green), legend=(:topright))
        plot!(p, x, y1, label = "Eq.(12)", line = (:red, 1.5), marker = (:dot, :red) , legend=(:topright))

        ylabel!(L" \textrm{std}(E_{0})")

        str = "_bootstrad_std"
        file1 = replace(file, ".npz" => str*".pdf")
        savefig(p, file1)

end

function plot_skewness(file::String, l::Int )

    D = npzread(file)

    try
        p = plot(size = (300, 200))
        x = D["betas"]
        xlabel!(L"$\beta_{MH}$")

        z = D["eta"]

        @. model(x, p) = p[1]*x^(p[2]/2)
        p0 = [0.5, 0.5]

        fit = curve_fit(model, x[1:l], z[1:l], p0)
        f(x, a = fit.param[1], b = fit.param[2]) = a*x^(b/2)

        plot!(p, x, z, markershape = :circle, legend=(:topleft), color = "red", label = "epmpirical data", xaxis = (:log))
        #plot!(p, ylims = (0., 2.5))
        plot!(p, x, [f(i) for i in x], linewidth = 2, style = :dot, color = "black", label = "fit with α = $(round(fit.param[2], digits=2))")

        vline!([(x[l] + x[l+1])/2], style = :dot, linewidth = 2., color = "green", label = "model limit")

        ylabel!(L"$\eta(H)$")

        str = "_etas"
        file1 = replace(file, ".npz" => str*".pdf")
        savefig(p, file1)

    catch
        0
    end

end

function plot_betas(file::String, l::Int)

    D = npzread(file)
    x = 0.
    p = plot(size = (300, 200))

    Z = D["estimated_betas"]

    try
        x = D["annealing_times"]
    catch
        x = D["betas"]
        #plot!(p, x, x, style = :dot, linewidth = 2., color = "black", label = "1 to 1")
        #plot!(xlims = (0.04, 2.), ylims = (0, 1.))
        vline!([(x[l]+x[l+1])/2], style = :dot, linewidth = 2., color = "green", label = "model limit  β = $((x[l] + x[l+1])/2)")

        @. model(x, p) = p[1]+p[2]*x
        p0 = [0.5, 0.5]

        fit = curve_fit(model, x[1:l], Z[1:l], p0)
        f(x, a = fit.param[1], b = fit.param[2]) = a+x*b

        plot!(p, x, [f(i) for i in x], linewidth = 2, style = :dot, color = "blue", label = "linear model, slope $(round(fit.param[2], digits = 2))")
    end

    plot!(p, x, Z, legend=(:topright), color = "red", label = L"\beta")
    plot!(p, x, Z, label = false, marker = (:dot, :red),  color = "red")

    try
        x = D["annealing_times"]
        xlabel!(L"\textrm{annelaing} \ \textrm{time} \ \mu s")
    catch
        xlabel!(L"\beta_{MH}")
    end
    ylabel!(L"\beta \  \textrm{ estimated}")

    str = "_betas"
    file1 = replace(file, ".npz" => str*".pdf")
    savefig(p, file1)

end

function plot_p_values(file::String, l::Int)

    D = npzread(file)
    ZZ = ( D["minimum_from_data"] .- D["ground"])/abs(D["ground"])
    li = maximum(ZZ)
    p = plot(size = (400, 250))
    p1 = twinx()
    α = D["alpha"]

    x = 0.
    try
        x = D["annealing_times"]
        xlabel!(L"\textrm{annelaing} \ \textrm{time} \ \mu s")
    catch
        x = D["betas"]
        plot!(p, xaxis = (:log))
        vline!([(x[l]+x[l+1])/2], style = :dot, linewidth = 2., color = "green", label = "model limit")

        xlabel!(L"$\beta_{MH}$")
    end

    plot!(p, legend=(0.15, 0.45), ylims = (-0.01, 1.05*li))
    plot!(p1, ylims = (-0.01, 1.01), legend=(0.67, 0.90))
    plot!(p, x, ZZ, markershape = :square, linewidth = 2, markersise = 10., color = "blue", label = L"$H_s$", ylabel = L"$(H_{\min} - H_0)/|H_0|$", right_margin=12mm)

    Z1 = D["p_values_10"]
    plot!(p1, x, Z1, markershape = :circ, label = "p-val., α = 0.10", color = "orange", ylabel = "p - value", right_margin=12mm)

    Z = D["p_values"]
    plot!(p1, x, Z, markershape = :diamond, label = "p-val., α = $α", color = "red")

    Z2 = D["p_values_39"]
    plot!(p1, x, Z2, markershape = :star, label = "p-val., α = 0.39", color = "brown")

    str = "_p_values"
    file1 = replace(file, ".npz" => str*".pdf")
    savefig(file1)
end

function main(file::String)
    β_lim_ind = 0
    if occursin("artificial_trains_case1", file)
        β_lim_ind = 10
    elseif occursin("artificial_trains_short", file)
        β_lim_ind = 6
    end
    plot_minimal_energies(file)
    plot_bootstrad_std(file, β_lim_ind)
    plot_betas(file, β_lim_ind)
    plot_p_values(file, β_lim_ind)
    plot_skewness(file, β_lim_ind)

end

s = ArgParseSettings("description")
  @add_arg_table! s begin
  "file"
  arg_type = String
  help = "file with data to be plotted"
end

file = parse_args(s)["file"]
main(file)


text:
using IterTools, ProgressBars
using JLD
using LaTeXStrings
using PyPlot

include("synthetic.jl")
qrng_names = vcat(collect(keys(qrngs)), ["RR", "SO"])

SEED    = 42
Random.seed!(SEED)

RUNS    = 1
EPOCHS  = 1000
T       = n * EPOCHS

# synthetic d-dimensional least squares problem
wgen = rand(d);

function gen_syn_data_all()
    # generate the training set consisting of n examples
    # sampled iid uniformly from the underlying distribution,
    # and the same training set will be used for all runs,
    gaussian_samples = gen_syn_data(n, "IID Uniform", d+1, SEED)
    return gaussian_samples
end

function grad_errors(T::Int64, wtau::Vector{Float64}, fullgrad::Vector{Float64},
             samples::Vector{Any}, ordering::Vector{Int64})
    errs = zeros(T)
    acc = zeros(d)
    for m in 1:T
        (xi, yi) = samples[ordering[m]]
        acc += (xi'wtau - yi) * xi
        errs[m] = norm(acc/m - fullgrad)^2
    end
    return errs
end

println("generating data...")
samples_offline = gen_syn_data_all()

println("begin computing gradient errors...")
result = Dict()

w0 = ones(d)

X = hcat([pair[1] for pair in samples_offline]...)'
y = [pair[2] for pair in samples_offline]
full_grad = X'*(X*w0-y)/n

for qrng_name in qrng_names
    err_list = zeros(RUNS, T)

    for run in 1:RUNS
        if qrng_name in collect(keys(qrngs))
            distrib = qrngs[qrng_name]
            # generate QMC indices
            sample_indices = [Int(ceil(i)) for i in QuasiMonteCarlo.sample(T, 0, n, distrib)]
        else
            if qrng_name == "SO"
                Random.seed!(SEED + run)
                perm = randperm!(collect(1:n))
                sample_indices = repeat(perm, EPOCHS)

            elseif qrng_name == "RR"
                sample_indices = zeros(Int64, T)
                for epoch in 1:EPOCHS
                    Random.seed!(SEED + run + epoch)
                    perm = randperm!(collect(1:n))
                    sample_indices[(epoch-1)*n+1 : epoch*n] = perm
                end
            else
                ArgumentError("SGD variant does not exist.")
            end
        end

        err_list[run, :] = grad_errors(T, w0, full_grad, samples_offline, sample_indices)
    end
    result[qrng_name] = err_list
end

println("completed")
println("plotting...")

figure(figsize = (6,6))
colors = ["#1f77bf", "#ff7f0e", "#17becf", "#7f7f7f"]
for (q, qrng_name) in enumerate(["IID Uniform", "Sobol", "RR", "SO"])

    err_list = result[qrng_name]
    err_mean = mean(err_list, dims=1)'

    # use moving average to smooth out the curves as RR and SO are highly oscillatory
    # towards the end
    using RollingFunctions
    ws = 2000
    err_mean = runmean(reshape(err_mean, T), ws)

    loglog(collect(1000:1000:T), err_mean[1000:1000:T], label=qrng_name, color=colors[q])
    ylim(bottom=1e-8)
end

legend(loc="lower left", fontsize=18)

xlabel(L"m", fontsize=20)
ylabel(L"\Vert \frac{1}{m}\Sigma_{t=0}^{m-1}\nabla R_n(w_0; x_t) - \nabla R_n(w_0)\Vert^2", fontsize=20)
title("Offline", fontsize=24)
xticks(fontsize=20)
yticks(fontsize=20)

tight_layout()
plot_fn = joinpath(pwd(),"figs", "qmc-graderr-offline.pdf")
savefig(plot_fn)
println("plot saved at $plot_fn")

max_stars_repo_name: EugeneLYC/qmc-ordering

max_issues_repo_name: EugeneLYC/qmc-ordering

max_forks_repo_name: EugeneLYC/qmc-ordering


text:
#=
try
    Pkg.clone("https://github.com/sisl/ProbabilisticGraphicalModels.jl.git")
catch e
    println("Exception when cloning ProbabilisticGraphicalModels.jl: $(e)")
end
=#


text:
# testing nearal ODE

using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))

# dudt2 = FastChain((x,p) -> x.^3,
#             FastDense(2,50,tanh),
#             FastDense(50,2))
# n_ode = NeuralODE(dudt2,tspan,Tsit5(),saveat=t)

dudt2 = Chain(x -> x.^3,
             Dense(2,50,tanh),
             Dense(50,2))
n_ode = NeuralODE(dudt2,tspan,Tsit5(),saveat=t)

function predict_n_ode(p)
  n_ode(u0,p)
end

function loss_n_ode(p)
    pred = predict_n_ode(p)
    loss = sum(abs2,ode_data .- pred)
    loss,pred
end

loss_n_ode(n_ode.p) # n_ode.p stores the initial parameters of the neural ODE

cb = function (p,l,pred;doplot=false) #callback function to observe training
  display(l)
  # plot current prediction against data
  if doplot
    pl = scatter(t,ode_data[1,:],label="data")
    scatter!(pl,t,pred[1,:],label="prediction")
    display(plot(pl))
  end
  return false
end

# Display the ODE with the initial parameter values.
cb(n_ode.p,loss_n_ode(n_ode.p)..., doplot=true)

res1 = DiffEqFlux.sciml_train(loss_n_ode, n_ode.p, ADAM(0.05), cb = cb, maxiters = 300)
cb(res1.minimizer,loss_n_ode(res1.minimizer)...;doplot=true)
res2 = DiffEqFlux.sciml_train(loss_n_ode, res1.minimizer, LBFGS(), cb = cb)
cb(res2.minimizer,loss_n_ode(res2.minimizer)...;doplot=true)

## Try GPU

# u0 = Float32[2.; 0.] |> gpu
# dudt = Chain(Dense(2,50,tanh),Dense(50,2)) |> gpu
#
# # p,re = DiffEqFlux.destructure(dudt)
# p,re = Flux.destructure(dudt)
#
# dudt_(u,p,t) = re(p)(u)
# prob = ODEProblem(ODEfunc, u0,tspan, p)
#
# # Runs on a GPU
# sol = solve(prob,Tsit5(),saveat=0.1)
#
# n_ode = NeuralODE(gpu(dudt2),tspan,Tsit5(),saveat=0.1)

#


text:
"""
Package: Forecast

Store results from the function `ar`

# Arguments
`varnames`        List of variable names
`order`           Order of Autoregressive Model
`ndims`           Number of dimensions
`Φ`               Collection of d by d matrices of coefficients
`coefficients`    Alias for Φ
`Φ0`              Constant
`constant`        Alias for Φ0
`Σ2`              ML variance/covariance Matrix
`variance`        Alias for Σ2
`Σ`               Variables Standard deviation
`stdev`           Alias for Σ
`x`               Original dataset
`fitted`          Fitted values
`residuals`       Prediction Error
`ic::Dict`        Collection of Information Criteria
`stats::Dict`     Collection of Statistics
`Φse`             Parameters Standard Error
`pse`             Alias for Φse
`Φ0se`            Constant Standard Error
`p0se`            Alias for Φ0se
`Φpv`             p-value for Parameters
`ppv`             Alias for Φpv
`Φ0pv`            p-value for Constant
`p0pv`            Alias for Φ0pv
`call::String`    Method called to generate AR
"""
mutable struct AR{T}
    varnames::Vector{String}
    order::Integer
    ndims::Integer
    Φ::Union{T,Array{T}}
    Φ0::Union{T,Vector{T}}
    Σ2::Union{T,Matrix{T}}
    Σ::Union{T,Vector{T}}
    x::Union{Array{T},DataFrame}
    fitted::Union{T,Array{T}}
    residuals::Union{T,Array{T}}
    ic::Dict
    stats::Dict
    Φse::Union{T,Array{T}}
    Φ0se::Union{T,Vector{T}}
    Φpv::Union{T,Array{T}}
    Φ0pv::Union{T,Vector{T}}
    call::String
end

function Base.show(io::IO, xar::AR)

    m,p = xar.ndims, xar.order

    Φ = expand(xar.Φ,(m,m,p))
    Φ0 = expand(xar.Φ0,(m,))
    Φpv = expand(xar.Φpv,(m,m,p))
    Φ0pv = expand(xar.Φ0pv,(m,))
    Σ2 = expand(xar.Σ2,(m,m))

    printstyled("Multivariate Autoregressive Model\n",bold=true,color=:underline)
    print("\n    ",xar.call,"\n")

    printstyled("\nResiduals Summary\n",bold=true,color=:underline)
    xar_sum = summarize(xar.residuals; varnames = xar.varnames)
    pretty_table(xar_sum.quantiles, nosubheader = true, show_row_number=false)
    pretty_table(xar_sum.moments, nosubheader = true, show_row_number=false)

    printstyled("\nCoefficients\n\n",bold=true,color=:underline)
    printstyled("Φ0\n",bold=true,color=:underline)
    pretty_table(string.(round.(Φ0,digits=3)) .* sigf.(Φ0pv),
                 tf = tf_matrix, noheader=true)

    for i in 1:p
        printstyled("Φ",i,"\n",bold=true,color=:underline)
        Φi = Φ[:,:,i]
        Φpvi = Φpv[:,:,i]
        pretty_table(string.(round.(Φi,digits = 3)) .* sigf.(Φpvi), tf = tf_matrix, noheader=true)
    end
    print("Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘^’ 0.1 ‘ ’ 1  and ‘+’ if fixed\n")

    # printstyled("\nCoefficients' Std. Error\n",bold=true,color=:underline)
    # printstyled("\nΦ0se\n",bold=true,color=:underline)
    # pretty_table(xar.Φ0se, tf = tf_matrix, noheader=true)
    # for i in 1:p
    #     printstyled("Φ",i,"se\n",bold=true,color=:underline)
    #     pretty_table(xar.Φse[:,:,i], tf = tf_matrix, noheader=true)
    # end

    printstyled("\nΣ2 Variance/Covariance Matrix\n",bold=true,color=:underline)
    pretty_table(Σ2, tf = tf_matrix, noheader=true)

    printstyled("\nInformation Criteria\n",bold=true,color=:underline)
    pretty_table(DataFrame(xar.ic), nosubheader = true, show_row_number=false)

    printstyled("\nStatistics\n",bold=true,color=:underline)
    pretty_table(DataFrame(xar.stats), nosubheader = true, show_row_number=false)

end

function sigf(pv::Real)::String

    if pv < 0.001 return  " ***"   end
    if pv < 0.01  return  " ** "   end
    if pv < 0.05  return  " *  "   end
    if pv < 0.1   return  " ^  "   end
    if pv < 1     return  "    "   end
    if pv == 1    return  " +  "   end

end


text:
################################################################################

mutable struct GraphCrawler
	root_node::Node
	available_dice::Vector{Die.Face}
	graphs::Dict{String,QuellerGraph}

	current::Node
	jump_stack::Vector{Node}

	options::Vector{CMD.Command}

	queller::QuellerState

	msg_buf::String

	function GraphCrawler(startgraph::String, graphs, available_dice)
		gc = new()

		gc.root_node = root(graphs[startgraph])
		gc.available_dice = available_dice
		gc.graphs = graphs

		return initialize!(gc)
	end
end

function initialize!(gc::GraphCrawler)
	gc.current = gc.root_node

	gc.jump_stack = Vector{Node}()
	sizehint!(gc.jump_stack, 10)

	gc.options = Vector{CMD.Command}()
	sizehint!(gc.options, 30)

	gc.queller = QuellerState(deepcopy(gc.available_dice))

	gc.msg_buf = ""

	autocrawl!(gc)

	return gc
end

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

function autocrawl!(gc)
	# crawls the graph until it encounters node that requires interaction
	# store all node messages in a buffer

	gc.msg_buf = ""

	while !at_end(gc)
		gc.msg_buf = add2msgbuf(gc.msg_buf, gc.current, gc.queller)
		gc.current isa InteractiveNode && break
		gc.current = autonext!(gc, gc.current)
	end
end

autonext!(gc, node::NonInteractiveNode) = getnext(node)
autonext!(gc, node::StateInteractionNode) = getnext!(node,gc.queller)

function autonext!(gc, node::JumpToGraph)
	push!(gc.jump_stack, node)
	return root(gc.graphs[node.jump_graph])
end

function autonext!(gc, ::ReturnFromGraph)
	node = pop!(gc.jump_stack)
	return getnext(node)
end

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

get_available_dice(gc) = gc.queller.available_dice
at_end(gc) = gc.current isa End
getinteraction(gc) = (gc.msg_buf, getopt(gc.current))

function proceed!(gc, opt)
	push!(gc.options, opt)
	gc.current = gc.current isa GetAvailableDice ?
		getnext!(gc.current, gc.queller, opt) : getnext(gc.current, opt)
	autocrawl!(gc)
end

function undo!(gc)
	isempty(gc.options) && return false # nothing to undo

	# undo latest options
	options = gc.options
	pop!(options)

	# reset state and replay remaining options
	initialize!(gc)
	for opt in options
		proceed!(gc, opt)
	end

	return true
end

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

add2msgbuf(buf,node::UseActiveDie,queller) = buf*"\n"*getmsg(node,queller)*"\n"
add2msgbuf(buf,node::InteractiveNode,queller) = buf*"\n"*getmsg(node)
add2msgbuf(buf,node::Node,queller) = buf

function add2msgbuf(buf, node::Union{BinaryCondition,MultipleChoice}, queller)
	die = queller.active_die
	!isnothing(die) && (buf = buf*"\n"*"Using $(Die.article(die.use_as)) $(string(die.use_as)):\n")
	return buf*"\n"*getmsg(node)
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
"""
$(TYPEDEF)
Encapsulation structure for a DFGNode (Variable or Factor) in Graphs.jl graph.
"""
mutable struct GraphsNode
    index::Int
    dfgNode::DFGNode
end
const FGType = Graphs.GenericIncidenceList{GraphsNode,Graphs.Edge{GraphsNode},Dict{Int,GraphsNode},Dict{Int,Array{Graphs.Edge{GraphsNode},1}}}

mutable struct GraphsDFG{T <: AbstractParams} <: AbstractDFG{T}
    g::FGType
    description::String
    userId::String
    robotId::String
    sessionId::String
    userData::Dict{Symbol, String}
    robotData::Dict{Symbol, String}
    sessionData::Dict{Symbol, String}
    nodeCounter::Int64
    labelDict::Dict{Symbol, Int64}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
    # GraphsDFG{T}(x...) where T <: AbstractParams = new{T}(x...)
end

function GraphsDFG( g::FGType=Graphs.incdict(GraphsNode,is_directed=false),
                    d::String="Graphs.jl implementation",
                    n::Int64=0,
                    l::Dict{Symbol, Int64}=Dict{Symbol, Int64}(),
                    a::Vector{Symbol}=Symbol[];
                    userId::String = "DefaultUser",
                    robotId::String = "DefaultRobot",
                    sessionId::String = "Session_$(string(uuid4())[1:6])",
                    userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                    robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                    sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                    params::T=NoSolverParams()) where T <: AbstractParams
    # Validate the userId, robotId, and sessionId
    Base.depwarn("GraphsDFG is voetstoots and will no longer be maintaind, use LightDFG", :GraphsDFG)
    !isValidLabel(userId) && error("'$userId' is not a valid User ID")
    !isValidLabel(robotId) && error("'$robotId' is not a valid Robot ID")
    !isValidLabel(sessionId) && error("'$sessionId' is not a valid Session ID")
    return GraphsDFG{T}(g, d, userId, robotId, sessionId, userData, robotData, sessionData, n, l, a, params)
end

function GraphsDFG{T}(  g::FGType=Graphs.incdict(GraphsNode,is_directed=false),
                        d::String="Graphs.jl implementation",
                        n::Int64=0,
                        l::Dict{Symbol, Int64}=Dict{Symbol, Int64}(),
                        a::Vector{Symbol}=Symbol[];
                        userId::String = "DefaultUser",
                        robotId::String = "DefaultRobot",
                        sessionId::String = "Session_$(string(uuid4())[1:6])",
                        userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                        robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                        sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                        params::T=T()) where T <: AbstractParams
    # Validate the userId, robotId, and sessionId
    Base.depwarn("GraphsDFG is voetstoots and will no longer be maintaind, use LightDFG", :GraphsDFG)
    !isValidLabel(userId) && error("'$userId' is not a valid User ID")
    !isValidLabel(robotId) && error("'$robotId' is not a valid Robot ID")
    !isValidLabel(sessionId) && error("'$sessionId' is not a valid Session ID")
    return GraphsDFG{T}(g, d, userId, robotId, sessionId, userData, robotData, sessionData, n, l, a, params)
end

GraphsDFG(description::String,
          userId::String,
          robotId::String,
          sessionId::String,
          userData::Dict{Symbol, String},
          robotData::Dict{Symbol, String},
          sessionData::Dict{Symbol, String},
          solverParams::AbstractParams) =
          GraphsDFG(Graphs.incdict(GraphsNode,is_directed=false), description, userId, robotId, sessionId, userData, robotData, sessionData, 0, Dict{Symbol, Int64}(), Symbol[], solverParams)


text:
using RecipesBase
using ColorTypes
import PlotUtils: cgrad

#using Compat
#using Compat: range

const mygreen = RGBA{Float64}(151/255,180/255,118/255,1)
const mygreen2 = RGBA{Float64}(113/255,161/255,103/255,1)
const myblue = RGBA{Float64}(74/255,144/255,226/255,1)

@recipe function plot(w::Fields.Nodes{T,NX,NY}) where {T,NX,NY}
  grid --> :none
  ratio := 1
  linewidth --> 1
  legend --> :none
  framestyle --> :frame
  levels --> range(minimum(w.data),stop=maximum(w.data),length=16)
  @series begin
    seriestype --> :contour
    transpose(w.data)
  end
end

@recipe function plot(x::AbstractArray{S,1},y::AbstractArray{S,1},w::Fields.Nodes{T,NX,NY};trim=0) where {S,T,NX,NY}
      grid --> :none
      ratio := 1
      linewidth --> 1
      legend --> :none
      framestyle --> :frame
      levels --> range(minimum(w.data),stop=maximum(w.data),length=16)
      @series begin
        seriestype --> :contour
        x[1+trim:end-trim],y[1+trim:end-trim],transpose(w.data[1+trim:end-trim,1+trim:end-trim])
      end
end

@recipe function plot(q::T) where {T <: Union{Fields.Edges,Fields.NodePair}}
    #wx = Fields.Nodes(Dual,(NX,NY))
    #wy = Fields.Nodes(Dual,(NX,NY))
    #shift!((wx,wy),q)
    layout := (1,2)
    seriestype --> :contour
    grid --> :none
    ratio := 1
    linewidth --> 1
    legend --> :none
    framestyle --> :frame
    @series begin
      subplot := 1
      #levels --> range(minimum(wx.data),stop=maximum(wx.data),length=16)
      levels --> range(minimum(q.u),stop=maximum(q.u),length=16)
      #transpose(wx.data)
      transpose(q.u)
    end

    @series begin
      subplot := 2
      #levels --> range(minimum(wy.data),stop=maximum(wy.data),length=16)
      levels --> range(minimum(q.v),stop=maximum(q.v),length=16)
      #transpose(wy.data)
      transpose(q.v)
    end
end

@recipe function plot(b::Body)
    x = [b.x; b.x[1]]
    y = [b.y; b.y[1]]
    linecolor --> mygreen
    fillrange --> 0
    fillcolor --> mygreen
    ratio := 1
    legend := :none
    grid := false
    x := x
    y := y
    ()
end

function RecipesBase.RecipesBase.apply_recipe(plotattributes::Dict{Symbol, Any}, bl::BodyList)
    series_list = RecipesBase.RecipeData[]
    for b in bl
        append!(series_list, RecipesBase.RecipesBase.apply_recipe(copy(plotattributes), b) )
    end
    series_list
end


text:
module DB
import Mongoc
import BSON
include("Initial.jl")
using .Initial
include("Node.jl")
using .Node

dbName = "Physics-Simulation"
client = nothing
database = nothing
collection = nothing

function insertToDB(data, id=nothing)
    document = juliaStructToBSON(data, id)

    push!(collection, document)
end

function connectToDB()
    global client = Mongoc.Client("localhost", 27017)
    global database = client[dbName]
    global collection = database[collectionName]
end

function juliaStructToBSON(data, id)
    json = ""

    for obj in data
        mass, radius, position, velocity = obj.mass, obj.radius, obj.position, obj.velocity

        json *= """{
            "mass":  $mass,
            "raduis": $radius,
            "position": $position,
            "velocity": $velocity
        },"""
    end

    if id ≠ nothing
        return Mongoc.BSON("""{
            "node": [""" * json[1:end - 1] * """],
            "_id": $id
        }""")
    else
        return Mongoc.BSON("{\"node\": [" * json[1:end - 1] * "]}")
    end
end

export connectToDB, insertToDB
end

max_stars_repo_name: PseudoDevelopers/Physics-Simulations

max_issues_repo_name: PseudoDevelopers/Physics-Simulations

max_forks_repo_name: PseudoDevelopers/Physics-Simulations


text:
# ---
# title: 113. Path Sum II
# id: problem113
# author: zhwang
# date: 2022-01-23
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/path-sum-ii/description/>
# hidden: true
# ---
#
# Given a binary tree and a sum, find all root-to-leaf paths where each path's
# sum equals the given sum.
#
# **Note:**  A leaf is a node with no children.
#
# **Example:**
#
# Given the below binary tree and `sum = 22`,
#
#
#
#           **5**
#          **/ \**
#         **4   8**
#        **/**   / **\**
#       **11**  13  **4**
#      /  **\**    **/** \
#     7    **2**  **5**   1
#
#
# Return:
#
#
#
#     [
#        [5,4,11,2],
#        [5,8,4,5]
#     ]
#
#
#
## @lc code=start
using LeetCode

function path_sum_113(node::Union{Nothing,TreeNode}, target::Int)::Vector{Vector{Int}}
    isnothing(node) && return Vector{Vector{Int}}[]
    isnothing(node.left) &&
        isnothing(node.right) &&
        return node.val != target ? Vector{Vector{Int}}[] : [[target]]
    paths = path_sum_113(node.left, target - node.val)
    paths = vcat(paths, path_sum_113(node.right, target - node.val))
    for path in paths
        pushfirst!(path, node.val)
    end
    return paths
end

## @lc code=end


text:
# ## Standard line colors

#md # ![FILE_NAME.png](images/FILE_NAME.png)

## ref https://jp.mathworks.com/matlabcentral/fileexchange/35229-matlab-plot-gallery-standard-line-colors?focused=6792895&tab=example

using SpecialFunctions

x = 0:0.2:10;
y0 = besselj.(0,x);
y1 = besselj.(1,x);
y2 = besselj.(2,x);
y3 = besselj.(3,x);
y4 = besselj.(4,x);
y5 = besselj.(5,x);
y6 = besselj.(6,x);

plot(x, y0, c=:red)
plot!(x, y1, c=:green)
plot!(x, y2, c=:blue)
plot!(x, y3, c=:cyan)
plot!(x, y4, c=:magenta)
plot!(x, y5, c=:yellow)
plot!(x, y6, c=:black)


text:
using ModularIndices
using Test, OffsetArrays, StaticArrays

@testset "0-d" begin
    for A in (fill(1), (x=Array{Int,0}(undef);x[]=1;x))
        @test A[Mod(5)] == 1
        @test A[Mod(-1)] == 1
        @test A[Mod(2:3)] == [1, 1]
        @test A[Mod(2, 30)] == [1, 1]
        @test A[Mod([5])] == [1]
        @test A[Mod(SVector(5))] == [1]
        A[Mod(2)] = 2
        @test A[] == 2
    end
end

@testset "1-d" begin
    for A in (1:10, collect(1:10), rand(10), @SVector(rand(10)), @MVector(rand(10)))
        @test A[Mod(11)] == A[1]
        @test A[Mod(12)] == A[2]
        @test A[Mod(1)] == A[1]
        @test A[Mod(1,2)] == A[[1,2]]
        @test A[Mod(11:15)] == A[1:5]
        @test A[Mod(SVector(11,12))] == A[[1,2]]
    end

    A = collect(1:5)
    A[Mod(7)] = 0
    @test A == [1,0,3,4,5]
end

@testset "2-d" begin
    for A in (rand(5,5), @SMatrix(rand(5,5)), view(rand(10,10), 1:5, 1:5))
        @test A[Mod(11), 2] == A[1, 2]
        @test A[2, Mod(11)] == A[2, 1]

        @test A[Mod(11:12), 2] == A[1:2, 2]
        @test A[Mod(11,12), 1] == A[1:2, 1]
    end

    A = OffsetArray(reshape(-10:13,3,:), -1, 0)
    # axes(A) == (Base.IdentityUnitRange(0:2), Base.IdentityUnitRange(1:8))
    @test A[Mod(3), 1] == A[0, 1]
    @test A[0, Mod(9)] == A[0, 1]
end

@testset "Invalid constructions" begin
    @test_throws MethodError Mod(1.5)
    @test_throws MethodError Mod([1.5, 2.5])
    @test_throws MethodError Mod(1.5, 2.5)
    @test_throws MethodError Mod(1:.1:2)
    @test_throws MethodError Mod((1,2,3))
end


text:
module Pyramids

using LinearAlgebra, FFTW, Images, Interpolations, Colors
using DSP # I need this only for conv2

export ImagePyramid, PyramidType, ComplexSteerablePyramid, LaplacianPyramid, GaussianPyramid
export subband, toimage, update_subband, update_subband!, test

"abstract supertype for the variety of pyramids that can be constructed."
abstract type PyramidType end
abstract type SimplePyramid <: PyramidType end

"""Type that indicates a complex steerable pyramid. [1]

[1] http://www.cns.nyu.edu/~eero/steerpyr/"""
struct ComplexSteerablePyramid <: PyramidType
end

"""Type that indicates a Laplacian pyramid. [1]

[1] persci.mit.edu/pub_pdfs/pyramid83.pdf"""
struct LaplacianPyramid <: SimplePyramid
end

"""Type that indicates a Gaussian pyramid. [1]

[1] http://persci.mit.edu/pub_pdfs/RCA84.pdf"""
struct GaussianPyramid <: SimplePyramid
end

"""Type that represents a concrete pyramidal representation of a given input image. Each type of pyramid has its own parameters. The basic construction method is, for example

```pyramid = ImagePyramid(im, ComplexSteerablePyramid(), scale=0.5^0.25)```

See the code for more information on optional arguments."""
struct ImagePyramid
    scale::Real
    num_levels::Int
    num_orientations::Int
    t::PyramidType
    pyramid_bands::Dict

    function ImagePyramid(pyr::ImagePyramid)
        return deepcopy(pyr)
    end

    function ImagePyramid(im::Array, t::ComplexSteerablePyramid; scale=0.5, min_size=15, num_orientations=8, max_levels=23, twidth=1)

        t = ComplexSteerablePyramid()
        scale = scale
        num_orientations = num_orientations

        im_dims = size(im)
        h = im_dims[1]
        w = im_dims[2]

        num_levels = min(ceil.(log2(minimum([h w]))/log2(1/scale) - (log2(min_size)/log2(1/scale))),max_levels);

        pyramid_bands, mtx, harmonics = build_complex_steerable_pyramid(im,
            num_levels, num_levels, order=num_orientations-1,
            twidth=twidth, scale=scale)

        return new(scale,num_levels,num_orientations,t,pyramid_bands)
    end

    function ImagePyramid(im::Array, t::GaussianPyramid; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
        this = new()

        pyramid_bands, num_levels = generate_gaussian_pyramid(im,
                min_size=min_size, max_levels=max_levels, filter=filter)

        scale = 0.5
        pyramid_bands = pyramid_bands
        num_orientations = 1
        num_levels = num_levels

        return new(scale,num_levels,num_orientations,t,pyramid_bands)
    end

    function ImagePyramid(im::Array, t::LaplacianPyramid; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])

        pyramid_bands, num_levels = generate_laplacian_pyramid(im,
                min_size=min_size, max_levels=max_levels, filter=filter)

        scale = 0.5
        pyramid_bands = pyramid_bands
        num_orientations = 1
        num_levels = num_levels

        return new(scale,num_levels,num_orientations,t,pyramid_bands)
    end
    # I am not sure this is needed
    # function ImagePyramid(pyramid_bands::Dict, scale, t, num_levels, num_orientations)
    #     return new(scale,num_levels,num_orientations,t,pyramid_bands)
    # end
end

##############################
# Functions

"""
    subband(pyramid::ImagePyramid, level; orientation=())

Returns the sub-band of the `pyramid` at `level` and at `orientation.` If orientation is not provided, `subband` assumes that the pyramid is not oriented.

Level 0 is always the high frequency residual."""
function subband(pyramid::ImagePyramid, level; orientation = ())
    if isempty(orientation)
        return pyramid.pyramid_bands[level]
    else
        return pyramid.pyramid_bands[level][orientation]
    end
end

function update_subband!(pyramid::ImagePyramid, level, new_subband; orientation = ())
    if isempty(orientation)
        pyramid.pyramid_bands[level] = copy(new_subband)
    else
        pyramid.pyramid_bands[level][orientation] = copy(new_subband)
    end

    return pyramid
end

function update_subband(pyramid::ImagePyramid, level, new_subband; orientation = ())
    newpyramid = ImagePyramid(pyramid)
    return update_subband!(newpyramid, level, new_subband, orientation=orientation)
end

"Converts a pyramid to a 2-D array (not of type `Image`!)"
function toimage(pyramid::ImagePyramid)
    if typeof(pyramid.t) <: ComplexSteerablePyramid
        im = reconstruct_complex_steerable_pyramid(pyramid, scale=pyramid.scale)
    elseif typeof(pyramid.t) <: LaplacianPyramid
        im = reconstruct_laplacian_pyramid(pyramid)
    elseif typeof(pyramid.t) <: GaussianPyramid
        im = subband(pyramid, 0)
    else
        error("Unsupported pyramid type $(typeof(pyramid.T))")
    end

    return im
end

##############################
# Private functions for building Gaussian and Laplacian pyramids.

function convolve_reflect(im, filter)
    filter_len = length(filter)
    filter_offset::Int = (filter_len-1)/2

    padded_im = zeros((size(im,1) + filter_len*2), (size(im,1) + filter_len*2))

    padded_im[(filter_len+1):(end-filter_len), (filter_len+1):(end-filter_len)] = im

    padded_im[(filter_len+1):(end-filter_len), 1:filter_len] =
        reverse(padded_im[(filter_len+1):(end-filter_len),
                (filter_len+2):(2*filter_len+1)], dims=2)
    padded_im[(filter_len+1):(end-filter_len), (end-filter_len+1):end] =
            reverse(padded_im[(filter_len+1):(end-filter_len),
                (end-2*filter_len):(end-filter_len-1)], dims=2)
    padded_im[1:filter_len, (filter_len+1):(end-filter_len),] =
        reverse(padded_im[(filter_len+2):(1+2*filter_len),
                (filter_len+1):(end-filter_len)], dims=1)
    padded_im[(end-filter_len+1):end, (filter_len+1):(end-filter_len)] =
            reverse(padded_im[(end-2*filter_len):(end-filter_len-1),
                (filter_len+1):(end-filter_len)], dims=1)

    new_im = DSP.conv2(filter, filter, padded_im)

    new_im = new_im[(1+filter_len+filter_offset):(end-filter_len-filter_offset), (1+filter_len+filter_offset):(end-filter_len-filter_offset)]
    return new_im
end

function reduce(im; filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    new_im = convolve_reflect(im, filter)
    return new_im[1:2:end, 1:2:end]
end

function expand(im; filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    new_im = zeros((size(im, 1)*2, size(im, 2)*2))
    new_im[1:2:end, 1:2:end] = im

    new_im = convolve_reflect(new_im, filter)
    return new_im
end

function generate_gaussian_pyramid(im; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    im = convert(Array{Float64}, copy(im))

    pyramid_bands = Dict{Integer, Array}()

    im_dims = collect(size(im))
    num_levels = min(max_levels, ceil.(Int, log2(minimum(im_dims)) - log2(min_size)))

    for i = 1:num_levels
        pyramid_bands[i-1] = im
        im = reduce(im, filter=filter)
    end

    return (pyramid_bands, num_levels)
end

function generate_laplacian_pyramid(im; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    im = convert(Array{Float64}, copy(im))

    pyramid_bands = Dict{Integer, Array}()

    im_dims = collect(size(im))

    num_levels = min(max_levels, ceil.(Int, log2(minimum(im_dims)) - log2(min_size)))

    filter_offset::Int = (length(filter)-1)/2

    for i = 1:num_levels
        reduced_im = reduce(im, filter=filter)
        next_im = expand(reduced_im)
        diff_im = im - next_im

        pyramid_bands[i-1] = diff_im

        im = reduced_im
    end

    pyramid_bands[num_levels] = im

    return (pyramid_bands, num_levels)
end

function reconstruct_laplacian_pyramid(pyramid::ImagePyramid)
    output_im::Array{Float64} = subband(pyramid, pyramid.num_levels)

    for i = (pyramid.num_levels-1):-1:0
        output_im = expand(output_im)
        output_im += subband(pyramid, i)
    end

    return output_im
end

##############################
# Private functions for building complex steerable pyramids.

function construct_steering_matrix(harmonics, angles; even = true)
    numh = 2*length(harmonics) - any(harmonics .== 0)

    imtx = zeros(length(angles), numh)
    col = 1
    for h in harmonics
        args = h * angles

        if h == 0
            imtx[:,col] = ones(angles)
            col += 1
        elseif !even
            imtx[:,col] = sin.(args)
            imtx[:,col+1] = -cos.(args)
            col += 2
        else
            imtx[:,col] = cos.(args)
            imtx[:,col+1] = sin.(args)
            col += 2
        end
    end

    r = rank(imtx)

    if (r != numh) && (r != length(angles))
        warning("Matrix is not full rank")
    end

    return pinv(imtx)
end

function raisedcosine(width=1, position=0, values=[0,1])
    sz = 256 # arbitrary
    X = collect( pi * (-sz-1:1) / (2 * sz) )
    Y = @. values[1] + ( (values[2]-values[1]) * cos(X)^2 )

    Y[1] = Y[2]
    Y[sz+3] = Y[sz+2]
    X = @.  position + (2*width/pi) * (X + pi/4)

    return (X, Y)
end

function build_complex_steerable_pyramid(im,
            height, nScales; order=3, twidth=1, scale=0.5)
    pyramid_bands = Dict{Integer, Union{Array, Dict{Integer, Array}}}()
    num_orientations = order + 1
    # generate steering matrix and harmonics info
    if mod.(num_orientations, 2) == 0
        # on 0.6 it was   ((0:((num_orientations/2)-1))*2+1)'
        harmonics = (1:2:(num_orientations/2-1)*2+1)
    else
        #on 0.6 it was ((0:((num_orientations-1)/2))*2)'
        harmonics = 0:2:(num_orientations-1)
    end

    steeringmatrix = construct_steering_matrix(harmonics,
                pi*(0:(num_orientations-1))/num_orientations, even=true)

    imdft = fftshift(fft(im))
    im_dims = size(im)
    ctr = ceil.(Int, (im_dims .+ 0.5) ./ 2. )

    angle = broadcast(atan,
            (((1:im_dims[1]) .- ctr[1]) ./ (im_dims[1]/2)),
            (((1:im_dims[2]) .- ctr[2]) ./ (im_dims[2]/2))') #SLOOW
    log_rad = broadcast((x,y) -> log2(sqrt(x^2 + y^2)),
            (((1:im_dims[1]) .- ctr[1]) ./ (im_dims[1]/2)),
            (((1:im_dims[2]) .- ctr[2]) ./ (im_dims[2]/2))') #SLOOW
    log_rad[ctr[1], ctr[2]] = log_rad[ctr[1], ctr[2]-1]

    Xrcos, Yrcos = raisedcosine(twidth, (-twidth/2), [0 1])
    Yrcos = sqrt.(Yrcos)
    YIrcos = @. sqrt(1 - Yrcos^2)

    # generate high frequency residual
    Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()) )
    Yrcosinterpolant_ext = extrapolate(Yrcosinterpolant,Flat())
    hi0mask = map(Yrcosinterpolant_ext,log_rad)
    hi0dft =  imdft .* hi0mask;

    pyramid_bands[0] = ifft(ifftshift(hi0dft));

    YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
    YIrcosinterpolant_ext = extrapolate(YIrcosinterpolant,Flat())
    lo0mask = map(YIrcosinterpolant_ext,log_rad)
    lo0dft = imdft .* lo0mask

    for ht = height:-1:1
        Xrcos = Xrcos .- log2(1/scale)
        order = num_orientations - 1

        cnst = (2^(2*order))*(factorial(order)^2)/(num_orientations*factorial(2*order))

        Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()))
        Yrcosinterpolant_ext = extrapolate(Yrcosinterpolant,Flat())
        himask = map(Yrcosinterpolant_ext,log_rad)

        # loop through each orientation band
        pyramid_level = Dict{Integer, Array}()

        for b in 1:num_orientations

            banddft = map(zip(lo0dft,himask,angle) ) do (lo0dft_i,himask_i,angle_i)
                ang = angle_i-pi*(b-1)/num_orientations
                a = (abs(mod(pi+ang, 2*pi) - pi) < pi/2) *
                            (2*sqrt(cnst) * (cos(ang)^order))
                (complex(0,-1)^(num_orientations-1) *
                             lo0dft_i * himask_i) * a
            end
            pyramid_level[b] = ifft(ifftshift(banddft))
        end

        pyramid_bands[height-ht+1] = pyramid_level

        dims = size(lo0dft)
        ctr = ceil.(Int, (dims .+ 0.5)./2)

        lodims = round.(Int, im_dims[1:2].* (scale^(nScales-ht+1)) )

        loctr = ceil.(Int, (lodims.+0.5)./2)
        lostart = @. ctr - loctr+1
        loend = @. lostart + lodims - 1

        log_rad = log_rad[lostart[1]:loend[1], lostart[2]:loend[2]]
        angle = angle[lostart[1]:loend[1], lostart[2]:loend[2]]
        lodft = lo0dft[lostart[1]:loend[1], lostart[2]:loend[2]]
        YIrcos = @. abs(sqrt(1 - Yrcos^2))

        YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
        YIrcosinterpolant_ext = extrapolate(YIrcosinterpolant,Flat())
        lomask = map(YIrcosinterpolant_ext,log_rad)

        lo0dft = lomask .* lodft
    end

    pyramid_bands[height+1] = real(ifft(ifftshift(lo0dft)))

    return (pyramid_bands, steeringmatrix, harmonics)
end

function make_angle_grid(sz, phase=0, origin=-1)
    if length(sz) == 1
        sz = [sz, sz]
    end

    if origin == -1
        origin = (sz + 1)/2
    end

    xramp = ones(round(Int, sz[1]), 1) * collect((1:sz[2]) .- origin[2])'
    yramp = collect((1:sz[1]) .- origin[1]) * ones(1, round.(Int, sz[2]))

    res = atan.(yramp, xramp)

    res = mod.(res .+ (pi-phase), 2*pi) .- pi

    return res
end

function reconstruct_steerable_pyramid(pyr::ImagePyramid;
        levs="all", bands="all", twidth=1, scale=0.5)
    dims = collect(size(subband(pyr, 0)))
    im_dft = zeros(Complex{Float64}, size(subband(pyr, 0)))

    ctr = ceil.(Int, (dims .+ 0.5) ./ 2)

    angle = broadcast(atan,
            (((1:dims[1]) .- ctr[1]) ./ (dims[1]/2)),
            (((1:dims[2]) .- ctr[2]) ./ (dims[2]/2))') #SLOOW
    log_rad = broadcast((x,y) -> log2(sqrt.(x.^2 + y.^2)),
        (((1:dims[1]) .- ctr[1]) ./ (dims[1]/2)),
        (((1:dims[2]) .- ctr[2]) ./ (dims[2]/2))') #SLOOW
    log_rad[ctr[1], ctr[2]] = log_rad[ctr[1], ctr[2]-1]
    log_rad0 = log_rad
    angle0 = angle

    Xrcos, Yrcos = raisedcosine(twidth, (-twidth/2), [0 1])
    Yrcos = sqrt.(Yrcos)
    YIrcos = @. sqrt(1 - Yrcos^2)

    # Start with low frequency residual
    low_dft = fftshift(fft(subband(pyr, pyr.num_levels+1)))

    lodims = collect(size(low_dft))
    loctr = ceil.(Int, (lodims .+ 0.5)./2)
    lostart = @. ctr - loctr + 1
    loend = @. lostart + lodims - 1

    log_rad = log_rad0[lostart[1]:loend[1], lostart[2]:loend[2]]
    angle = angle0[lostart[1]:loend[1], lostart[2]:loend[2]]

    Xrcos = Xrcos .- (log2(1/scale)*pyr.num_levels)
    YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
    YIrcosinterpolant_ext = extrapolate(YIrcosinterpolant,Flat())
    lomask = map(YIrcosinterpolant_ext,log_rad)

    im_dft[lostart[1]:loend[1], lostart[2]:loend[2]] += low_dft .* lomask

    # Accumulate mid-bamds
    for level = pyr.num_levels:-1:1
        lodims = collect(size(subband(pyr, level, orientation=1)))
        loctr = ceil.(Int, (lodims .+ 0.5) ./ 2)
        lostart = @. ctr - loctr+1
        loend = @. lostart + lodims - 1

        log_rad = copy(log_rad0[lostart[1]:loend[1], lostart[2]:loend[2]])
        angle = copy(angle0[lostart[1]:loend[1], lostart[2]:loend[2]])

        Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()))
        Yrcosinterpolant_ext = extrapolate(Yrcosinterpolant,Flat())
        himask = map(Yrcosinterpolant_ext,log_rad)

        Xrcos = Xrcos .+ log2(1/scale)
        order = pyr.num_orientations - 1

        cnst = ((complex(0,1))^(pyr.num_orientations-1)) *
            sqrt((2^(2*order))*(factorial(order)^2)/
                    (pyr.num_orientations*factorial(2*order)))

        for orientation = 1:pyr.num_orientations
            band_dft = fftshift(fft(subband(pyr, level, orientation=orientation)))
            ang = angle .- (pi*(orientation-1)/pyr.num_orientations)
            angle_mask = 2 .* (abs.(mod.(pi .+ ang, 2*pi) .- pi) .< pi/2) .*
                (cnst .* (cos.(ang).^order))

            im_dft[lostart[1]:loend[1], lostart[2]:loend[2]] += band_dft .*
             himask .* angle_mask
        end

        YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
        YIrcosinterpolant_ext = extrapolate(YIrcosinterpolant,Flat())
        lomask = map(YIrcosinterpolant_ext,log_rad)

        # Everything must be scaled by the low-frequency mask
        im_dft[lostart[1]:loend[1], lostart[2]:loend[2]] .*= lomask
    end

    # Add high frequency residual
    Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()) )
    Yrcosinterpolant_ext = extrapolate(Yrcosinterpolant,Flat())
    hi0mask = map(Yrcosinterpolant_ext,log_rad)
    im_dft += fftshift(fft(subband(pyr, 0))) .* hi0mask;

    return real(ifft(ifftshift(im_dft)))
end

# make the complex steerable pyramid real
function convert_complex_steerable_pyramid_to_real(pyramid::ImagePyramid;
            levs="all", bands="all", twidth=1, scale=0.5)
    pyramid = ImagePyramid(pyramid)
    num_levels = pyramid.num_levels
    num_orientations = pyramid.num_orientations

    for nsc in 1:num_levels
        dims = collect(size(subband(pyramid, nsc, orientation=1)))
        ctr = ceil.(Int, (dims .+ 0.5)./ 2)
        ang = make_angle_grid(dims, 0, ctr)
        ang[ctr[1], ctr[2]] = -pi/2

        for nor = 1:num_orientations
            ch = subband(pyramid, nsc, orientation=nor)

            ang0 = pi*(nor-1)/num_orientations
            xang = mod.(ang .- (ang0+pi), 2*pi) .- pi

            # this creates an angular mask
            amask = 2.0 .* (abs.(xang) .< pi/2) .+ (abs.(xang) .== pi/2)
            amask[ctr[1], ctr[2]] = 1.0
            amask[:,1] .= 1.0
            amask[1,:] .= 1.0

            # and masks the fft by it
            amask = fftshift(amask)
            ch = ifft(amask.*fft(ch))
            ch = 0.5 .* real(ch)

            # then creates a new pyramid
            update_subband!(pyramid, nsc, ch, orientation=nor)
        end
    end

    # and returns it as a real
    return pyramid
end

function reconstruct_complex_steerable_pyramid(pyramid::ImagePyramid;
            levs="all", bands="all", twidth=1, scale=0.5)
    real_pyramid = convert_complex_steerable_pyramid_to_real(pyramid,
        levs=levs, bands=bands, twidth=twidth, scale=scale)
    return reconstruct_steerable_pyramid(real_pyramid, levs=levs,
        bands=bands, twidth=twidth, scale=scale)
end

end


text:
#unit-multiplication.jl

#testing chunk multiplication
a = 0xFFFFFFFFFFFFFFFF
b = 0xFFFFFFFFFFFFFFFF

#unit test mult_exact.

################################################################################
# found errors
x = Unum{4,5}(0x0000000000000003, 0x5555555600000000, 0x0000, 0x0001, 0x001F)
y = Unum{4,5}(3)
#should be approximately 16, but usually gets 8.
@test x * y > Unum{4,5}(16)

x = Unum{4,7}(o64, Unums.ArrayNum{7}([0x3fff_ffff_ffff_fffe, 0xc000_0000_0000_0000]), z16, z16, 0x004F)
y = Unums.ArrayNum{7}([o64, z64])
Unums.frac_mul!(x, y, Val{2}, Val{true})
@test x.fraction == Unums.ArrayNum{7}([0x3fff_ffff_ffff_ffff, 0xffff_ffff_ffff_fffe])

#bizzare ubit setting.
p = Unum{4,7}(0x0000000000000001, UInt64[0x8000000000000000,0x0000000000000000], 0x0000, 0x0001, 0x0001)
q = Unum{4,7}(0x0000000000000007, UInt64[0x4000000000000000,0x0000000000000000], 0x0000, 0x0002, 0x0002)
@test Unums.is_exact(p * q)

#wierd multiplication results.

x = Ubound(neg_one(Unum{3,4}), pos_one(Unum{3,4}))
@test x * x == x
@test is_zero((x ^ 2).lower)
@test x ^ 3 == x
@test is_zero((x ^ 4).lower)

#making sure our heuristic multiplication works.

#single ulp.
UT = Unum{3,4}

right_exact = UT(2)
right_ulp   = Unums.inner_ulp!(UT(2))

@test Unums.outer_exact!(right_exact * right_ulp) == UT(4)

#single ulp that behaves differently.

lhs = Unum{3,5}(0x0000000000000001, 0x5555555500000000, 0x0001, 0x0002, 0x001F)
rhs = Unum{3,5}(0x000000000000000C, 0x3800000000000000, 0x0000, 0x0003, 0x001F)

@test lhs * rhs ≹ Unum{3,5}(13)

#double ulp.
x = UT(2)
Unums.inner_ulp!(x)
@test Unums.outer_exact!(x * x) == UT(4)

#some strange exact multiplications
#identified during bounds testing
x = Unums.inner_exact!(mmr(Unum{3,5}))
y = Unum{3,5}(0x0000000000000001, 0xFFFFFFFF00000000, 0x0001, 0x0001, 0x001F)
@test Unums.mul_exact(x, y, z16) == mmr(Unum{3,5})

#do another heuristic multiplication.  lub((1 - 1.5) ^ 2) == 2.25
x = Unum{3,5}(0x0000000000000001, 0x0000000000000000, 0x0001, 0x0001, 0x0000)
@test Unums.lub(x * x) == Unum{3,5}(0x0000000000000001, 0x2000000000000000, 0x0000, 0x0000, 0x001F)

#fixed an issue with ubound parity decisions.
x = Unum{3,5}(0x0000000000000006, 0x9CE4A90000000000, 0x0003, 0x0002, 0x001C)
y = Ubound(Unum{3,5}(0x0000000000000004, 0xEC58764D00000000, 0x0003, 0x0004, 0x001F), Unum{3,5}(0x0000000000000004, 0xEC58764B00000000, 0x0003, 0x0004, 0x001F))
@test x * y == Ubound(Unum{3,5}(0x0000000000000008, 0x8D0B111800000000, 0x0001, 0x0004, 0x001F), Unum{3,5}(0x0000000000000008, 0x8D0B112200000000, 0x0001, 0x0004, 0x001F))

#double ubound parity error.
x = Ubound(Unum{3,5}(0x0000000000000004, 0xB06061DA00000000, 0x0003, 0x0003, 0x001F), Unum{3,5}(0x0000000000000004, 0xB06061C700000000, 0x0003, 0x0003, 0x001F))
y = Ubound(Unum{3,5}(0x000000000000000C, 0x0A97BA7000000000, 0x0001, 0x0003, 0x001F), Unum{3,5}(0x000000000000000C, 0x0A97BAE700000000, 0x0001, 0x0003, 0x001F))

@test x * y == Ubound(Unum{3,5}(0x0000000000000003, 0xC2446A3500000000, 0x0003, 0x0001, 0x001F), Unum{3,5}(0x0000000000000003, 0xC244695500000000, 0x0003, 0x0001, 0x001F))

#inexact unum multiplication forgot to reset the fsize to maximum for inner target.

#=
x = Unum{3,6}(0x000000000000000E, 0xC7FFFFFE00000000, 0x0003, 0x0003, 0x001F)
y = Unum{3,6}(0x0000000000000018, 0xA000000000000000, 0x0002, 0x0004, 0x0002)
@test x * y == one(Unum{3,6})

x = Unum{3,6}(0x0000000000000018, 0xA000000000000000, 0x0002, 0x0004, 0x0002)
y = Ubound(Unum{3,6}(0x000000000000000E, 0xC800000300000000, 0x0003, 0x0003, 0x001F), Unum{3,6}(0x000000000000000E, 0xC7FFFFFE00000000, 0x0003, 0x0003, 0x001F))

@test x * y == Ubound(Unum{3,6}(0x0000000000000030, 0x727FFFFE60000000, 0x0001, 0x0005, 0x003F), Unum{3,6}(0x0000000000000030, 0x727FFFFF2FFFFFFF, 0x0001, 0x0005, 0x003F))

uft1 = unum_easy(Unum{4,6}, zero(UInt16), frac1, 1)
uft2 = unum_easy(Unum{4,6}, zero(UInt16), frac2, 1)
xft1 = convert(Float64, uft1)
xft2 = convert(Float64, uft2)
ifm3 = (reinterpret(UInt64, xft1 * xft2) << 12) & f64_frac_mask
ufm3 = (Unums.__mult_exact(uft1, uft2)).fraction & f64_frac_mask
@test ifm3 == ufm3
=#
#result:  A coding error caused this carry function in the multiplication
#operation to fail to trigger in the last 32-bit segment.  This has been fixed.

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Workflows
using Test

@testset "Workflows.jl" begin
    # Write your tests here.
end


text:
# utilities for Mesh2D
function opposing_face(mesh::Mesh2D, elem, local_face)
    @assert 0 < elem ≤ mesh.nelems
    @assert 0 < local_face ≤ 4

    neighbors = mesh.face_neighbors
    face = mesh.elem_faces[elem, local_face]

    if neighbors[face, 1] == elem && neighbors[face, 2] == local_face
        neigh_elem, neigh_lface = neighbors[face, 3], neighbors[face, 4]
    elseif neighbors[face, 3] == elem && neighbors[face, 4] == local_face
        neigh_elem, neigh_lface = neighbors[face, 1], neighbors[face, 2]
    end
    if neigh_elem == 0
        neigh_lface = local_face
    end
    reversed = neighbors[face, 5] == -1
    return neigh_elem, neigh_lface, reversed
end

function interior_faces(mesh::Mesh2D)
    tags = mesh.boundary_tags
    offset = mesh.face_boundary_offset

    if tags[1] ≠ 0
        return nothing
    else
        return mesh.face_boundary[1:offset[2]-1]
    end
end

function boundary_faces(mesh::Mesh2D, btag::Int)
    tags = mesh.boundary_tags
    offset = mesh.face_boundary_offset
    loc = findfirst(tags .== btag)
    if loc == nothing
        return nothing
    else
        return mesh.face_boundary[offset[loc]:offset[loc+1]-1]
    end
end

function boundary_face_info(mesh::Mesh2D, face::Int)
    face_neighbors = mesh.face_neighbors
    if face_neighbors[face, 1] == 0
        return tuple(face_neighbors[face, 3:4]...)
    else
        return tuple(face_neighbors[face, 1:2]...)
    end
end

function boundary_tag(mesh::Mesh2D, tag_name::Symbol)
    tag_names = mesh.boundary_tag_names
    tags = mesh.boundary_tags
    loc = findfirst(mesh.boundary_tag_names .== tag_name)
    if loc == nothing
        return nothing
    else
        return tags[loc]
    end
end

function unique_vertex_connectivity(mesh::Mesh2D, uvertno::Int)
    offset = mesh.uverts_offset
    conn = mesh.uverts_conn
    st, en = offset[uvertno], offset[uvertno+1]
    nconn = div(en - st, 2)
    return [tuple(conn[(st+(i-1)*2):(st+(i-1)*2+1)]...) for i = 1:nconn]
end

function vertex_coordinates(mesh::Mesh2D, elem)
    @assert 0 < elem ≤ mesh.nelems
    verts = mesh.elem_verts[elem, :]
    coords = mesh.coordinates
    return [tuple(coords[vert, :]...) for vert in verts]
end

# utilities for Mesh1D
# Note: for 1D meshes, faces are identical to vertices
function opposing_face(mesh::Mesh1D, elem, local_face)
    @assert 0 < elem ≤ mesh.nelems
    @assert 0 < local_face ≤ 2

    neighbors = mesh.vertex_neighbors
    face = mesh.elem_verts[elem, local_face]

    if neighbors[face, 1] == elem && neighbors[face, 2] == local_face
        neigh_elem, neigh_lface = neighbors[face, 3], neighbors[face, 4]
    elseif neighbors[face, 3] == elem && neighbors[face, 4] == local_face
        neigh_elem, neigh_lface = neighbors[face, 1], neighbors[face, 2]
    end
    if neigh_elem == 0
        neigh_lface = local_face
    end
    return neigh_elem, neigh_lface
end

function interior_faces(mesh::Mesh1D)
    tags = mesh.boundary_tags
    offset = mesh.vertex_boundary_offset

    if tags[1] ≠ 0
        return nothing
    else
        return mesh.vertex_boundary[1:offset[2]-1]
    end
end

function boundary_faces(mesh::Mesh1D, btag::Int)
    tags = mesh.boundary_tags
    offset = mesh.vertex_boundary_offset
    loc = findfirst(tags .== btag)
    if loc == nothing
        return nothing
    else
        return mesh.vertex_boundary[offset[loc]:offset[loc+1]-1]
    end
end

function boundary_face_info(mesh::Mesh1D, face::Int)
    face_neighbors = mesh.vertex_neighbors
    if face_neighbors[face, 1] == 0
        return tuple(face_neighbors[face, 3:4]...)
    else
        return tuple(face_neighbors[face, 1:2]...)
    end
end

function boundary_tag(mesh::Mesh1D, tag_name::Symbol)
    tag_names = mesh.boundary_tag_names
    tags = mesh.boundary_tags
    loc = findfirst(mesh.boundary_tag_names .== tag_name)
    if loc == nothing
        return nothing
    else
        return tags[loc]
    end
end

function unique_vertex_connectivity(mesh::Mesh1D, uvertno::Int)
    neighbors = mesh.vertex_neighbors
    vt = mesh.unique_verts[uvertno]
    if neighbors[vt, 1] == 0 && neighbors[vt, 3] ≠ 0
        return [tuple(neighbors[vt, 3:4]...)]
    elseif neighbors[vt, 1] ≠ 0 && neighbors[vt, 3] == 0
        return [tuple(neighbors[vt, 1:2]...)]
    else
        return [tuple(neighbors[vt, 1:2]...), tuple(neighbors[vt, 3:4]...)]
    end
end

function vertex_coordinates(mesh::Mesh1D, elem)
    @assert 0 < elem ≤ mesh.nelems
    verts = mesh.elem_verts[elem, :]
    coords = mesh.coordinates
    return [tuple(coords[vert, :]...) for vert in verts]
end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
module IOBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat

import Compat: UTF8String, view

const SUITE = BenchmarkGroup()

#################
# read (#12364) #
#################

function perf_read!(io)
    seekstart(io)
    x = 0
    while !(eof(io))
        x += read(io, UInt8)
    end
    return x
end

g = addgroup!(SUITE, "read", ["buffer", "stream", "string"])

testbuf = IOBuffer(randstring(RandUtils.SEED, 10^4))

g["read"]       = @benchmarkable perf_read!($testbuf)
g["readstring"] = @benchmarkable read($testbuf, String)

#################################
# serialization (#18633, #7893) #
#################################

g = addgroup!(SUITE, "serialization", ["buffer", "string"])

function serialized_buf(x)
    io = IOBuffer()
    serialize(io, x)
    return io
end

teststrings = [randstring(RandUtils.SEED, 32) for i=1:10^3]
teststrings_buf = serialized_buf(teststrings)

g["serialize", "Vector{String}"] = @benchmarkable serialize(io, $teststrings) setup=(io=IOBuffer())
g["deserialize", "Vector{String}"] = @benchmarkable (seek($teststrings_buf, 0); deserialize($teststrings_buf))

testdata = rand(RandUtils.SEED,1000,1000)
testdata_buf = serialized_buf(testdata)

g["serialize", "Matrix{Float64}"] = @benchmarkable serialize(io, $testdata) setup=(io=IOBuffer())
g["deserialize", "Matrix{Float64}"] = @benchmarkable (seek($testdata_buf, 0); deserialize($testdata_buf))

end # module


text:
# Structs
struct telegraph_points
     i₀::Vector{Int64}
     i₁::Vector{Int64}
end

struct telegraph_time_dist
     t₀::Vector{Float64}
     t₁::Vector{Float64}
end


text:
using StochasticDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear
Random.seed!(100)
prob = prob_sde_linear

## Solve and plot
println("Solve and Plot")
sol = solve(prob,EM(),dt=1//2^(4))
sol = solve(prob,RKMil(),dt=1//2^(4))
sol = solve(prob,RKMil_General(),dt=1//2^(4))
sol = solve(prob,SRI(),dt=1//2^(4))
sol = solve(prob,SRIW1(),dt=1//2^(4))
trajectories = 100
## Convergence Testing
println("Convergence Test on Linear")
dts = (1//2) .^(9:-1:4) #14->7 good plot with higher num Monte

sim = test_convergence(dts,prob,EM(),trajectories=trajectories)

sim2 = test_convergence(dts,prob,RKMil(),trajectories=trajectories)

sim21 = test_convergence(dts,prob,RKMil_General(),trajectories=trajectories)

sim3 = test_convergence(dts,prob,SRI(),trajectories=trajectories)

#TEST_PLOT && plot(plot(sim),plot(sim2),plot(sim3),layout=@layout([a b c]),size=(1200,600))

@test abs(sim.𝒪est[:l2]-.5) + abs(sim2.𝒪est[:l∞]-1) + abs(sim3.𝒪est[:final]-1.5)<.441  #High tolerance since low dts for testing!

@test abs(sim.𝒪est[:l2]-.5) + abs(sim21.𝒪est[:l∞]-1) + abs(sim3.𝒪est[:final]-1.5)<.441  #High tolerance since low dts for testing!

# test reinit
integrator = init(prob,EM(),dt=1//2^(4))
solve!(integrator)
reinit!(integrator)
solve!(integrator)

# test reinit
prob2 = SDEProblem((u,p,t)->prob.f(u,p,t),prob.g,prob.u0,prob.tspan)
integrator = init(prob2,EM(),dt=1//2^(4), tstops = [1//2], saveat = [1//3])
solve!(integrator)
reinit!(integrator)
solve!(integrator)


text:
# we need to import these so that we can extend them
import ParserCombinator: execute, success, failure
# defines hash() and == for us
using AutoHashEquals

# a simple example of a matcher that calls a sub-matcher, gets the result, and
# capitalizes the first letter found.  so it expects to receive a string.

# much of the work is already done for us by the Delegate related code in
# matchers.jl.  to use that, we need to follow the conventions of having a
# matcher with a matcher field and a state with a state field (they can have
# more fields - those are the minimal requirements).

# every matcher has a name field, which by default is the name of the matcher
# (using @with_names changes this to the variable name - see calc.jl)
@auto_hash_equals mutable struct Case<:Delegate
    name::Symbol
    matcher::Matcher
    Case(matcher) = new(:Case, matcher)
end

struct CaseState<:DelegateState
    state::State
end

# the Delegate code handles the initial call (with CLEAN), and failure.  all
# we need to do is add handling for success

function success(k::Config, m::Case, s, t, i, r::Value)
    # we don't care about the old state for this matcher (s), but we need to
    # save the child state (t), so that it can be used in backtracking.
    new_s = CaseState(t)
    # get the string contents from the child matcher
    # (nicer code would check this was a list containing a single string)
    contents::AbstractString = r[1]
    new_contents = uppercase(contents[1:1]) * contents[2:end]
    # and build the response from this matcher (see types.jl)
    Success(new_s, i, Any[new_contents])
end

@testset "case" begin

# now let's test that
@test parse_one("foo", Case(p".*")) == ["Foo"]

# to see what's happening in more detail, add debug logging:
@test parse_dbg("foo", Trace(Case(p".*"))) == ["Foo"]
# which gives
#  1:foo        00 Trace->Case
#  1:foo        01  Case->Pattern
#  4:           01  Case<-["foo"]
#  4:           00 Trace<-["Foo"]
# and we can see that Case calls Pattern, as expected; that the Pattern
# returns foo; and that Case transforms that to Foo.

end


text:
export vandpoly #export function vandpoly, make it public
Reduce.Preload() ## required otherwise Reduce has a fit -?

## ##The Vandermonde matrix method #https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/05Interpolation/vandermonde/#:~:text=The%20Vandermonde%20method%20is%20the,Gaussian%20elimination%20or%20PLU%20decomposition.

##Create macro genvandpoly with:
# INPUTS:
#  index - coefficient index we want to return
#  leny - number of points available for leny-1 degree polynomial interpolation
#  x - array with x coordinates of leny points
#  y - array with y coordinates of leny points
# OUTPUTS:
#  coefs[index] - the [$index]'s coefficient

macro genvandpoly(index,leny,x,y) ## generate poly functions using the Vandermonde matrix method
    leny = :($leny) #length of y #expression with number leny
    symX = [] #symmetric X matrix declaration
    symY = [] #symmetric Y matrix declaration
    A = Array{Expr}(undef, leny, leny) #Array of expressions (ie. Vandermonde matrix)
##we push an expression xᵢ,yᵢ into symX,symY
    for i = 1:leny #for all integers i from 1 to leny
        push!(symX, Meta.parse("x$i")) #meta.parse - parses a string to an expression
        push!(symY, Meta.parse("y$i")) #meta.parse - parses a string to an expression
    end
## ## Assemble A matrix #we go over all elements in A ∈ℜⁿˣⁿ matrix #we go over all elements in A and construct them as expressions
    for i = 1:leny  #i is the row number (ie. equation or point)
        for j = 1:leny #j is the column number (ie. power of x)
            pwr = leny-j #power is the total number of columns - index of current column
            var = symX[i] #variable that we need to raise to the power, took from symX
            A[i,j] = :($var^$pwr) #each element of A is an expression var^pwr, with var and pwr being numbers
        end
    end
## ## invert Vandermonde matrix "A" to find coefficients of polynomial
    coefs = Algebra.:*(Algebra.inv(A),symY) #coefs stores polynomial coefficients and is coefs=A\symY
    for i = 1:leny #for all integers i in [1,leny] (ie. all elements in coefs)
        str = repr(coefs[i]) #str is a string representation of element i from coefs
        for j = 1:leny
            str = replace(str,"y"*string(j)=>"$y"*"["*string(j)*"]") #replaces into str, the string yⱼ with the numerical value of y[j]
            str = replace(str,"x"*string(j)=>"$x"*"["*string(j)*"]") #replaces into str, the string xⱼ with the numerical value of x[j]
        end
        coefs[i] = Meta.parse(chop(str, head=2, tail=1)) #remove first 2 and last 1 characters from string str and parse it as an expression
    end
##return the coefficient with index number "index"
    return coefs[:($index)]
end

## Function vandpoly
#  INPUTS:
# x - vector with x coordinates from n points
# y - vector with y coordiantes from n points
#  OUTPUTS:
# Polynomial([cn,cn-1,...,c2,c1]) - Polynomial expression of the form c₁xⁿ⁻¹+c₂xⁿ⁻²+...+cₙ₋₁x+cₙ

function vandpoly(x,y) #used to return a polynomial of degree n-1 out of n datapoints in ℜ²
    xLength = length(x) #xLength is the length of x vector

    if xLength == 2 #if we have only 2 points we fit y=c₁x¹+c₂
        c1 = @genvandpoly(1, 2,:($x),:($y)) #c1 is the coefficient of x¹
        c2 = @genvandpoly(2, 2,:($x),:($y)) #c2 is the coefficient of x⁰
        return Polynomial([c2,c1]) #Return the polynomial c₁x¹+c₂

    elseif xLength == 3 #if we have only 3 points we fit y=c₁x²+c₂x¹+c₃
        c1 = @genvandpoly(1, 3,:($x),:($y)) #c1 is the coefficient of x²
        c2 = @genvandpoly(2, 3,:($x),:($y)) #c2 is the coefficient of x¹
        c3 = @genvandpoly(3, 3,:($x),:($y)) #c3 is the coefficient of x⁰
        return Polynomial([c3,c2,c1]) #Return the polynomial c₁x²+c₂x¹+c₃

    elseif xLength == 4 #if we have only 4 points we fit y=c₁x³+c₂x²+c₃x¹+c₄
        c1 = @genvandpoly(1, 4,:($x),:($y)) #c1 is the coefficient of x³
        c2 = @genvandpoly(2, 4,:($x),:($y)) #c2 is the coefficient of x²
        c3 = @genvandpoly(3, 4,:($x),:($y)) #c3 is the coefficient of x¹
        c4 = @genvandpoly(4, 4,:($x),:($y)) #c4 is the coefficient of x⁰
        return Polynomial([c4,c3,c2,c1]) #Return the polynomial c₁x³+c₂x²+c₃x¹+c₄

    elseif xLength == 5 #if we have only 5 points we fit y=c₁x⁴+c₂x³+c₃x²+c₄x¹+c₅
        c1 = @genvandpoly(1, 5,:($x),:($y)) #c1 is the coefficient of x⁴
        c2 = @genvandpoly(2, 5,:($x),:($y)) #c2 is the coefficient of x³
        c3 = @genvandpoly(3, 5,:($x),:($y)) #c3 is the coefficient of x²
        c4 = @genvandpoly(4, 5,:($x),:($y)) #c4 is the coefficient of x¹
        c5 = @genvandpoly(5, 5,:($x),:($y)) #c5 is the coefficient of x⁰
        return Polynomial([c5,c4,c3,c2,c1]) #Return the polynomial c₁x⁴+c₂x³+c₃x²+c₄x¹+c₅

    else #if we have more than 5 points we simply perform C\y
        Polynomial(inv([xval^j for xval in eachindex(x), j = 0:xlength - 1]) * y) # create poly coefficient by inverting array of linear equations
        println("Warning, function behaves slow when exceeding more than 5 datapoints")
    end
end


text:
function compute_normal(vertex,face)
    """
        compute_normal - compute the normal of a triangulation

          [normal,normalf] = compute_normal(vertex,face);

          normal(i,:) is the normal at vertex i.
          normalf(j,:) is the normal at face j.

          Copyright (c) 2004 Gabriel Peyre
    """
    epsilon = eps(Float64)
    nface = size(face,2)
    nvert = size(vertex,2)
    normal = zeros(3,nvert)

    # unit normals to the faces
    normalf = crossp(vertex[:,face[2,:]] - vertex[:,face[1,:]],
                     vertex[:,face[3,:]] - vertex[:,face[1,:]])
    d = sqrt(sum(normalf.^2,1))
    d[d .< epsilon] = 1
    normalf = normalf./repeat(d,outer=(3,1))

# unit normal to the vertex
    for i in 1:nface
        f = face[:,i]
        for j in 1:3
            normal[:,f[j]] = normal[:,f[j]] + normalf[:,i]
        end
    end

    # normalize
    d = sqrt(sum(normal.^2,1))
    d[d .< epsilon] = 1
    normal = normal./repeat(d,outer=(3,1))

    # enforce that the normal are outward
    v = vertex - repeat(mean(vertex,1),outer=(3,1))
    s = sum(v.*normal,2)
    if sum(s .> 0) < sum(s .< 0)
        # flip
        normal = -normal
        normalf = -normalf
    end
    return normal
end

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

function crossp(x,y)

    # x and y are (m,3) dimensional
    z = copy(x)
    z[1,:] = x[2,:].*y[3,:] - x[3,:].*y[2,:]
    z[2,:] = x[3,:].*y[1,:] - x[1,:].*y[3,:]
    z[3,:] = x[1,:].*y[2,:] - x[2,:].*y[1,:]

    return z
end


text:
export CoefficientHomotopy, CoefficientHomotopyCache, gamma

@enum CoeffHomotopyActiveCoeffs begin
    CH_COEFFS_EVAL
    CH_COEFFS_DT
    CH_COEFFS_UNKNOWN
end

"""
    CoefficientHomotopy(support, start_coeffs, target_coeffs)

Construct the homotopy ``H(x, t) = ∑_{a ∈ Aᵢ} (c_a t + (1-t)d_a) x^a`` where ``c_a`` are
the start coefficients and ``d_a`` the target coefficients.
"""
struct CoefficientHomotopy{S<:AbstractSystem,T1,T2} <: AbstractHomotopy
    system::S
    start_coeffs::Vector{Vector{T1}}
    target_coeffs::Vector{Vector{T2}}
end

function CoefficientHomotopy(
    exponents_polys::Union{<:Vector{<:Matrix},<:MPPolys},
    start::Vector{<:Vector},
    target::Vector{<:Vector},
)
    @assert all(length.(start) .== length.(target) .== _nterms.(exponents_polys))

    system = SPSystem(exponents_polys, start)
    # StaticPolynomials possibly changes the order of the terms.
    # Therefore the coefficient vectors are maybe no more correct.
    # We correct this by applying the permutation applied to the columns
    # of the support
    start_coeffs = deepcopy(start)
    target_coeffs = deepcopy(target)
    for (i, f) in enumerate(system.system.polys)
        p = SP.permutation(f)
        permute!(start_coeffs[i], p)
        permute!(target_coeffs[i], p)
    end

    CoefficientHomotopy(system, start_coeffs, target_coeffs)
end

_nterms(A::Matrix) = size(A, 2)
_nterms(f::MPPoly) = MP.nterms(f)

(H::CoefficientHomotopy)(x, t, c = cache(H, x, t)) = evaluate(H, x, t, c)

"""
    CoefficientHomotopyCache

An simple cache for `StartTargetHomotopy`s consisting of the caches for the
start and target system as well as a `Vector` and a `Matrix`.
"""
struct CoefficientHomotopyCache{C<:AbstractSystemCache,T} <: AbstractHomotopyCache
    system::C
    coeffs::Vector{Vector{T}}
    diffs::Vector{Vector{T}}
    t::Base.RefValue{ComplexF64}
    active_coeffs::Base.RefValue{CoeffHomotopyActiveCoeffs}
end

function cache(H::CoefficientHomotopy{S,T1,T2}, x, t) where {S,T1,T2}
    system = cache(H.system, x, t)
    U = promote_type(T1, T2, typeof(t))
    coeffs = map(c -> convert.(U, c), H.start_coeffs)
    diffs = map(H.start_coeffs, H.target_coeffs) do start, target
        convert(Vector{U}, start - target)
    end
    active_coeffs = Ref(CH_COEFFS_UNKNOWN)
    CoefficientHomotopyCache(system, coeffs, diffs, Ref(complex(NaN)), active_coeffs)
end

Base.size(H::CoefficientHomotopy) = size(H.system)

function update_coeffs!(cache::CoefficientHomotopyCache, H::CoefficientHomotopy, t)
    if t == cache.t[]
        if cache.active_coeffs[] != CH_COEFFS_EVAL
            set_coefficients!(H.system, cache.coeffs)
            cache.active_coeffs[] = CH_COEFFS_EVAL
        end
        return nothing
    end

    if isreal(t)
        rt = real(t)
        rt1 = 1.0 - t
        @inbounds for k = 1:length(cache.coeffs)
            c = cache.coeffs[k]
            start = H.start_coeffs[k]
            target = H.target_coeffs[k]
            for i in eachindex(c)
                c[i] = rt * start[i] + rt1 * target[i]
            end
        end
    else
        t1 = 1.0 - t
        @inbounds for k = 1:length(cache.coeffs)
            c = cache.coeffs[k]
            start = H.start_coeffs[k]
            target = H.target_coeffs[k]
            for i in eachindex(c)
                c[i] = t * start[i] + t1 * target[i]
            end
        end
    end
    set_coefficients!(H.system, cache.coeffs)
    cache.active_coeffs[] = CH_COEFFS_EVAL
    cache.t[] = t
    nothing
end

function update_coeffs_dt!(cache::CoefficientHomotopyCache, H::CoefficientHomotopy, t)
    if cache.active_coeffs[] != CH_COEFFS_DT
        set_coefficients!(H.system, cache.diffs)
        cache.active_coeffs[] = CH_COEFFS_DT
    end
    nothing
end

function evaluate!(u, H::CoefficientHomotopy, x, t, c::CoefficientHomotopyCache)
    update_coeffs!(c, H, t)
    @inbounds evaluate!(u, H.system, x, c.system)
    u
end

function dt!(u, H::CoefficientHomotopy, x, t, c::CoefficientHomotopyCache)
    update_coeffs_dt!(c, H, t)
    @inbounds evaluate!(u, H.system, x, c.system)
    u
end

function jacobian!(U, H::CoefficientHomotopy, x, t, c::CoefficientHomotopyCache)
    update_coeffs!(c, H, t)
    @inbounds jacobian!(U, H.system, x, c.system)
end

function evaluate_and_jacobian!(
    u,
    U,
    H::CoefficientHomotopy,
    x,
    t,
    c::CoefficientHomotopyCache,
)
    update_coeffs!(c, H, t)
    @inbounds evaluate_and_jacobian!(u, U, H.system, x, c.system)
    nothing
end

function evaluate(H::CoefficientHomotopy, x, t, c::CoefficientHomotopyCache)
    update_coeffs!(c, H, t)
    evaluate(H.system, x, c.system)
end

function dt(H::CoefficientHomotopy, x, t, c::CoefficientHomotopyCache)
    update_coeffs_dt!(c, H, t)
    evaluate(H.system, x, c.system)
end

function jacobian(H::CoefficientHomotopy, x, t, c::CoefficientHomotopyCache)
    update_coeffs!(c, H, t)
    jacobian(H.system, x, c.system)
end


text:
struct Assets <: Rewriter
    # Mapping of asset name (virtual path) to physical path
    assets::OrderedDict{AbstractPath, AbstractPath}
    isdirty::BitVector
end

function Assets(srcdir::AbstractPath; kwargs...)
    @assert isdir(srcdir)
    return Assets(srcdir, [relative(p, srcdir) for p in walkpath(srcdir)]; kwargs...)
end

function Assets(srcdir::AbstractPath, paths::AbstractVector{<:AbstractPath}; dstdir = p"assets/")
    assets = Dict(joinpath(dstdir, p) => absolute(joinpath(srcdir, p)) for p in paths)
    return Assets(assets)
end

function Assets(assets::Dict)
    isdirty = trues(length(assets))
    return Assets(assets, isdirty)
end

function postbuild(assets::Assets, project, builder::FileBuilder)
    for (i, (relp, srcpath)) in enumerate(assets.assets)
        dstpath = joinpath(absolute(builder.dir), relp)
        if !isfile(dstpath)
            mkpath(joinpath(builder.dir, parent(dstpath)))
            cp(srcpath, dstpath, force = true)
            assets.isdirty[i] = false
        end
    end
end

function filehandlers(assets::Assets, project, builder::FileBuilder)
    handlers = Dict()
    for (relp, fp) in assets.assets
        handlers[fp] = () -> cp(fp, joinpath(builder.dir, relp); force = true)
    end
    return handlers
end

function reset!(assets::Assets)
    fill!(assets.isdirty, true)
end


text:
"""
    ThreadedScans.scan!(op, xs::AbstractVector; [basesize]) -> xs

Compute inclusive scan of a vector `xs` with an associative binary function `op`
in parallel.
"""
function ThreadedScans.scan!(op, xs; basesize::Union{Nothing,Integer} = nothing)
    if basesize === nothing
        ntasks = Threads.nthreads()
    else
        ntasks = cld(length(xs), basesize)
    end
    return ThreadedScans.partitioned_hillis_steele!(op, xs; ntasks = ntasks)
end


text:
function distribute(graph::OptiGraph,to_workers::Vector{Int64};remote_name = :graph)
    error("distribute() has not yet been updated for the latest version of PipsSolver.jl")
end

#Distribute a modelgraph among workers.  Each worker should have the same master model.  Each worker will be allocated some of the nodes in the original modelgraph
# function distribute(graph::OptiGraph,to_workers::Vector{Int64};remote_name = :graph)
#     #NOTE: Linkconstraints keep their indices in new graphs, NOTE: Link constraint row index needs to match on each worker
#     #NOTE: Does not yet support subgraphs.  Aggregate first
#     #Create remote channel to store the nodes we want to send
#
#     #IDEA: Create a channel from the master process to each worker?
#     # channel_nodes = RemoteChannel(1)    #we will allocate and send nodes to workers
#     # channel_indices = RemoteChannel(1)
#     channel = RemoteChannel(1)
#     to_workers = sort(to_workers)
#
#     n_nodes = num_nodes(graph)
#     n_workers = length(to_workers)
#     nodes_per_worker = Int64(floor(n_nodes/n_workers))
#     nodes = all_nodes(mg)
#     node_indices = [getindex(graph,node) for node in nodes]
#
#
#     linkeqconstraints = OrderedDict()
#     linkineqconstraints = OrderedDict()
#     for edge in getedges(mg)
#         for (idx,link) in edge.linkeqconstraints
#             linkeqconstraints[idx] = link
#         end
#         for (idx,link) in edge.linkineqconstraints
#             linkineqconstraints[idx] = link
#         end
#     end
#
#     n_linkeq_cons =  length(linkeqconstraints)  #length(mg.linkeqconstraints)
#     n_linkineq_cons = length(linkineqconstraints) #length(mg.linkineqconstraints)
#
#     ineqlink_lb = zeros(n_linkineq_cons)
#     ineqlink_ub = zeros(n_linkineq_cons)
#
#     for (idx,link) in linkineqconstraints #mg.linkineqconstraints
#         if isa(link.set,MOI.LessThan)
#             ineqlink_lb[idx] = -Inf
#             ineqlink_ub[idx] = link.set.upper
#         elseif isa(link.set,MOI.GreaterThan)
#             ineqlink_lb[idx] = link.set.lower
#             ineqlink_ub[idx] = Inf
#         elseif isa(link.set,MOI.Interval)
#             ineqlink_lb[idx] = link.set.lower
#             ineqlink_ub[idx] = link.set.upper
#         end
#     end
#
#     #Allocate modelnodes onto provided workers
#     allocations = []
#     node_indices = []
#     j = 1
#     while  j <= n_nodes
#         if j + nodes_per_worker > n_nodes
#             push!(allocations,nodes[j:end])
#             push!(node_indices,[getindex(mg,node) for node in nodes[j:end]])
#         else
#             push!(allocations,nodes[j:j+nodes_per_worker - 1])
#             push!(node_indices,[getindex(mg,node) for node in nodes[j:j+nodes_per_worker - 1]])
#         end
#         j += nodes_per_worker
#     end
#
#     println("Distributing graph among workers: $to_workers")
#     remote_references = []
#     #Fill channel with sets of nodes to send
#     #TODO: Make this run in parallel
#     @sync begin
#         for (i,worker) in enumerate(to_workers)
#             @spawnat(1, put!(channel_nodes, allocations[i]))
#             @spawnat(1, put!(channel_indices, node_indices[i]))
#             ref1 = @spawnat worker begin
#                 Core.eval(Main, Expr(:(=), :nodes, take!(channel_nodes)))
#                 Core.eval(Main, Expr(:(=), :node_indices, take!(channel_indices)))
#             end
#             wait(ref1)
#             ref2 = @spawnat worker Core.eval(Main, Expr(:(=), remote_name, PipsSolver._create_worker_optigraph(getfield(Main,:nodes),getfield(Main,:node_indices),
#             n_nodes,n_linkeq_cons,n_linkineq_cons,ineqlink_lb,ineqlink_ub)))
#             push!(remote_references,ref2)
#         end
#         #return remote_references
#     end
#     return remote_references
# end
#
# function _create_worker_optigraph(optinodes::Vector{OptiNode},node_indices::Vector{Int64},n_nodes::Int64,n_linkeq_cons::Int64,n_linkineq_cons::Int64,
#     link_ineq_lower::Vector,link_ineq_upper::Vector)
#
#     graph = OptiGraph()
#     graph.node_idx_map = Dict{OptiNode,Int64}()
#
#     #Add nodes to worker's graph.  Each worker should have the same number of nodes, but some will be empty.
#     for i = 1:n_nodes
#         add_node!(graph)
#     end
#
#     #Populate models for given nodes
#     for (i,node) in enumerate(modelnodes)
#         index = node_indices[i]  #need node index in highest level
#         new_node = getnode(graph,index)
#         set_model(new_node,getmodel(node))
#     end
#     # We need the graph to have the partial constraints over graph nodes
#     # Add link constraints
#     # graph.linkeqconstraints = _add_linkeq_terms(modelnodes)
#     # graph.linkineqconstraints = _add_linkineq_terms(modelnodes)
#
#     #Setup new graph linkconstraints
#     linkeqconstraints = _add_linkeq_terms(modelnodes)
#     linkineqconstraints = _add_linkineq_terms(modelnodes)
#
#     #Need to match both equality and inequality
#     #Add linkconstraints, then fix indices
#     for (idx,link) in linkeqconstraints
#         cref = Plasmo.add_link_equality_constraint(graph,JuMP.ScalarConstraint(link.func,link.set);eq_idx = idx)
#
#         # #need to change the constraint index to match the original linkconstraint
#         # linkedge = cref.linkedge
#         # old_idx = graph.linkeqconstraint_index
#         #
#         # #change index
#         # linkedge.linkeqconstraints[idx] = linkedge.linkeqconstraints[old_idx]
#         #
#         # #we don't always want to delete
#         # if old_idx != idx &&  linkedge.linkeqconstraints[old_idx]
#         #     delete!(linkedge.linkeqconstraints,old_idx)
#         # end
#
#     end
#
#
#
#     for (idx,link) in linkineqconstraints
#         cref = Plasmo.add_link_inequality_constraint(graph,JuMP.ScalarConstraint(link.func,link.set);ineq_idx = idx)
#         # linkedge = cref.linkedge
#         # old_idx = graph.linkineqconstraint_index
#         # linkedge.linkineqconstraints[idx] = linkedge.linkineqconstraints[old_idx] #old_idx is the wrong index
#         #
#         # if old_idx != idx
#         #     delete!(linkedge.linkineqconstraints,old_idx)
#         # end
#     end
#
#     #Tell the worker how many linkconstraints the graph actually has
#     graph.obj_dict[:n_linkeq_cons] = n_linkeq_cons
#     graph.obj_dict[:n_linkineq_cons] = n_linkineq_cons
#     graph.obj_dict[:linkineq_lower] = link_ineq_lower
#     graph.obj_dict[:linkineq_upper] = link_ineq_upper
#     return graph
# end
#
# function _add_linkeq_terms(modelnodes::Vector{OptiNode})
#     linkeqconstraints = OrderedDict()
#     for node in modelnodes
#         partial_links = node.partial_linkeqconstraints
#         for (idx,linkconstraint) in partial_links
#             if !(haskey(linkeqconstraints,idx))   #create link constraint
#                 new_func = linkconstraint.func
#                 set = linkconstraint.set
#                 linkcon = LinkConstraint(new_func,set)
#                 linkeqconstraints[idx] = linkcon
#             else #update linkconstraint
#                 newlinkcon = linkeqconstraints[idx]
#                 nodelinkcon = node.partial_linkeqconstraints[idx]
#                 newlinkcon = LinkConstraint(newlinkcon.func + nodelinkcon.func,newlinkcon.set)
#                 linkeqconstraints[idx] = newlinkcon
#             end
#         end
#     end
#     return linkeqconstraints
# end
#
# function _add_linkineq_terms(modelnodes::Vector{OptiNode})
#     linkineqconstraints = OrderedDict()
#     for node in modelnodes
#         partial_links = node.partial_linkineqconstraints
#         for (idx,linkconstraint) in partial_links
#             if !(haskey(linkineqconstraints,idx))   #create link constraint
#                 new_func = linkconstraint.func
#                 set = linkconstraint.set
#                 linkcon = LinkConstraint(new_func,set)
#                 linkineqconstraints[idx] = linkcon
#             else #update linkconstraint
#                 newlinkcon = linkineqconstraints[idx]
#                 nodelinkcon = node.partial_linkineqconstraints[idx]
#                 newlinkcon = LinkConstraint(newlinkcon.func + nodelinkcon.func,newlinkcon.set)
#                 linkineqconstraints[idx] = newlinkcon
#             end
#         end
#     end
#     return linkineqconstraints
# end