text:
"""
Memoize functions for fitness evaluation
"""

"""
    memoize_fitness(f_fitness::Function, n_outputs::Int)

Memoize fitness function.
Because `foo_i` is auto-differentiated with ForwardDiff, our cache needs to
work when `x` is a `Float64` and a `ForwardDiff.Dual`.
"""
function memoize_fitness(f_fitness::Function, nx::Int, n_outputs::Int)
    last_x, last_f = ones(nx), ones(n_outputs)
    last_dx, last_dfdx = ones(nx), ones(nx)
    function f_fitness_i(i::Int, x::T...) where {T<:Real}
        if T == Float64
            if x != last_x
                last_x, last_f = x, f_fitness(x...)
            end
            return last_f[i]::T
        else
            if x != last_dx
                last_dx, last_dfdx = x, f_fitness(x...)
            end
            return last_dfdx[i]::T
        end
    end
    return [(x...) -> f_fitness_i(i, x...) for i in 1:n_outputs]
end

"""
    memoize_fitness_gradient(f_fitness::Function, n_outputs::Int, nx::Int)

Create memoized gradient computation with FiniteDifferences.forward_fdm()
"""
function memoize_fitness_gradient(
    f_fitness::Function,
    nx::Int,
    nfitness::Int,
    order::Int=2
)
	diff_f = "forward"
    return memoize_fitness_gradient(f_fitness, nx, nfitness, diff_f, order)
end

"""
    memoize_fitness_gradient(f_fitness::Function, nfitness::Int, diff_f::Function, order::Int=2)

Create memoized gradient computation with method specified by `diff_f`
    - `diff_f = "forward"` use `FiniteDifferences.forward_fdm()`
    - `diff_f = "central"` use `FiniteDifferences.central_fdm()`
    - `diff_f = "backward"` use `FiniteDifferences.backward_fdm()`
"""
function memoize_fitness_gradient(
    f_fitness::Function,
    nx::Int,
    nfitness::Int,
    diff_f::String,
    order::Int=2
)
	# create memoized version
	f_fitness_nosplat(x) = f_fitness(x...)
    last_x, last_f = ones(nx), zeros(nx,nfitness)

    # Memoized function for gradient and jacobians
    if cmp(diff_f, "forward")==0

        function foo_fwd(i::Int, df::AbstractVector{T}, x...) where {T} #::T...)# where {T<:Real}
            if x != last_x
                # update
                last_x = x
            	jac2 = FiniteDifferences.jacobian(FiniteDifferences.forward_fdm(order, 1), f_fitness_nosplat, x)
                last_f = jac2[1]
            end
            df[:] = last_f[i,:]
            return
        end
        return [(df, x...) -> foo_fwd(i, df, x...) for i in 1:nfitness]

    elseif cmp(diff_f, "central")==0

        function foo_cent(i::Int, df::AbstractVector{T}, x...) where {T} #::T...)# where {T<:Real}
            if x != last_x
                # update
                last_x = x
                jac2 = FiniteDifferences.jacobian(FiniteDifferences.central_fdm(order, 1), f_fitness_nosplat, x)
                last_f = jac2[1]
            end
            df[:] = last_f[i,:]
            return
        end
        return [(df, x...) -> foo_cent(i, df, x...) for i in 1:nfitness]

    elseif cmp(diff_f, "backward")==0

        function foo_bck(i::Int, df::AbstractVector{T}, x...) where {T} #::T...)# where {T<:Real}
            if x != last_x
                # update
                last_x = x
                jac2 = FiniteDifferences.jacobian(FiniteDifferences.backward_fdm(order, 1), f_fitness_nosplat, x)
                last_f = jac2[1]
            end
            df[:] = last_f[i,:]
            return
        end
        return [(df, x...) -> foo_bck(i, df, x...) for i in 1:nfitness]

    end
end

ext: jl

max_stars_repo_name: Yuricst/NLPSaUT

max_issues_repo_name: Yuricst/NLPSaUT

max_forks_repo_name: Yuricst/NLPSaUT


text:
const Node = Any

struct HyperEdge{N}
    label::String
    source::Array{N}
    target::Array{N}
end

struct HyperGraph{N} <: AbstractArray{HyperEdge, 1}
    edges::Array{HyperEdge{N}}
end

# minimum methods for abstract array; I delegate everything to the edges array
Base.size(HG::HyperGraph) = size(HG.edges)

Base.IndexStyle(::Type{<: HyperGraph}) = IndexLinear()

Base.getindex(HG::HyperGraph, i::Int) = getindex(HG.edges, i)

const HRGrammarRules = Dict{String,HyperGraph{String}}


text:
module AutonomousControl

using NLOptControl
using VehicleModels

include("CaseModule.jl")
using .CaseModule

export
      initializeAutonomousControl,
      autonomousControl,
      updateAutoParams!

"""
n=initializeAutonomousControl(c);
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/1/2017, Last Modified: 7/05/2017 \n
--------------------------------------------------------------------------------------\n
"""
function initializeAutonomousControl(c)
 pa=Vpara(x_min=c.m.Xlims[1],x_max=c.m.Xlims[2],y_min=c.m.Ylims[1],y_max=c.m.Ylims[2],sr_min=-0.18,sr_max=0.18);
 @unpack_Vpara pa
 XF=[c.g.x_ref, c.g.y_ref, NaN, NaN, NaN, NaN, NaN, NaN];
 XL=[x_min, y_min, NaN, NaN, psi_min, sa_min, u_min, NaN];
 XU=[x_max, y_max, NaN, NaN, psi_max, sa_max, u_max, NaN];
 CL = [sr_min, jx_min]; CU = [sr_max, jx_max];
 n=define(numStates=8,numControls=2,X0=copy(c.m.X0),XF=XF,XL=XL,XU=XU,CL=CL,CU=CU)
 n.s.tf_max=10.0;
 n.params=[pa];   # vehicle parameters

 # set mpc parameters
 initializeMPC!(n;FixedTp=c.m.FixedTp,PredictX0=c.m.PredictX0,tp=c.m.tp,tex=copy(c.m.tex),max_iter=c.m.mpc_max_iter);
 n.mpc.X0=[copy(c.m.X0)];
 n.mpc.plantEquations=ThreeDOFv2;
 n.mpc.modelEquations=ThreeDOFv2;

 # define tolerances
 X0_tol=[0.05,0.05,0.05,0.05,0.01,0.001,0.05,0.05];
 XF_tol=[c.m.sigma,c.m.sigma,NaN,NaN,NaN,NaN,NaN,NaN];
 defineTolerances!(n;X0_tol=X0_tol,XF_tol=XF_tol);

         # 1  2  3  4  5    6   7   8
 names = [:x,:y,:v,:r,:psi,:sa,:ux,:ax];
 descriptions = ["X (m)","Y (m)","Lateral Velocity (m/s)", "Yaw Rate (rad/s)","Yaw Angle (rad)", "Steering Angle (rad)", "Longitudinal Velocity (m/s)", "Longitudinal Acceleration (m/s^2)"];
 states!(n,names,descriptions=descriptions)

          # 1   2
 names = [:sr,:jx];
 descriptions = ["Steering Rate (rad/s)","Longitudinal Jerk (m/s^3)"];
 controls!(n,names,descriptions=descriptions)

 # dynamic constraints and additional constraints
 dx,con,tire_expr=ThreeDOFv2_expr(n)
 dynamics!(n,dx)
 constraints!(n,con)

 # solver settings
 solver_time = (c.m.solver==:Ipopt) ? :max_cpu_time : :maxtime_real
 SS=((:name=>c.m.solver),(:mpc_defaults=>true),(solver_time=>100.)) # let the solver have more time for the initalization

 # configure problem
 configure!(n;(Nck=c.m.Nck),(:integrationScheme=>c.m.integrationScheme),(:finalTimeDV=>true),(:solverSettings=>SS))
 x=n.r.x[:,1];y=n.r.x[:,2];psi=n.r.x[:,5]; # pointers to JuMP variables

 #################################
 # obstacle aviodance constraints
 ################################
 Q = size(c.o.A)[1]; # number of obstacles
 @NLparameter(n.mdl, a[i=1:Q] == copy(c.o.A[i]));
 @NLparameter(n.mdl, b[i=1:Q] == copy(c.o.B[i]));
 @NLparameter(n.mdl, X_0[i=1:Q] == copy(c.o.X0[i]));
 @NLparameter(n.mdl, Y_0[i=1:Q] == copy(c.o.Y0[i]));
 @NLparameter(n.mdl, speed_x[i=1:Q] == copy(c.o.s_x[i]));
 @NLparameter(n.mdl, speed_y[i=1:Q] == copy(c.o.s_y[i]));
 obs_params=[a,b,X_0,Y_0,speed_x,speed_y,Q];

 # obstacle postion after the intial postion
 X_obs=@NLexpression(n.mdl, [j=1:Q,i=1:n.numStatePoints], X_0[j] + speed_x[j]*n.tV[i]);
 Y_obs=@NLexpression(n.mdl, [j=1:Q,i=1:n.numStatePoints], Y_0[j] + speed_y[j]*n.tV[i]);

 # constraint on position
 obs_con=@NLconstraint(n.mdl, [j=1:Q,i=1:n.numStatePoints-1], 1 <= ((x[(i+1)]-X_obs[j,i])^2)/((a[j]+c.m.sm)^2) + ((y[(i+1)]-Y_obs[j,i])^2)/((b[j]+c.m.sm)^2));
 newConstraint!(n,obs_con,:obs_con);

 ####################
 # LiDAR constraints
 ###################
 # ensure that the final x and y states are near the LiDAR boundary
 @NLparameter(n.mdl, LiDAR_param_1==(c.m.Lr + c.m.L_rd)^2);
 @NLparameter(n.mdl, LiDAR_param_2==(c.m.Lr - c.m.L_rd)^2);
 LiDAR_edge_high = @NLconstraint(n.mdl,[j=1], (x[end]-x[1])^2+(y[end]-y[1])^2  <= LiDAR_param_1);
 LiDAR_edge_low = @NLconstraint(n.mdl,[j=1], (x[end]-x[1])^2+(y[end]-y[1])^2  >= LiDAR_param_2);
 newConstraint!(n,LiDAR_edge_high,:LiDAR_edge_high);
 newConstraint!(n,LiDAR_edge_low,:LiDAR_edge_low);

# constrain all state points to be within LiDAR boundary
 LiDAR_range = @NLconstraint(n.mdl, [j=1:n.numStatePoints-1], (x[j+1]-x[1])^2+(y[j+1]-y[1])^2 <= (c.m.Lr + c.m.L_rd)^2 );
 newConstraint!(n,LiDAR_range,:LiDAR_range);

 if goalRange!(n,c)   # relax LiDAR boundary constraints
    setvalue(LiDAR_param_1, 1e6)
    setvalue(LiDAR_param_2,-1e6)
else                  # relax constraints on the final x and y position
    for st=1:2;for k in 1:2; setRHS(n.r.xf_con[st,k], 1e6); end; end
end
 LiDAR_params=[LiDAR_param_1,LiDAR_param_2]
 #####################
 # objective function
 ####################
 # parameters
 if goalRange!(n,c)
  println("\n goal in range")
  @NLparameter(n.mdl, w_goal_param == 0.0)
  @NLparameter(n.mdl, w_psi_param == 0.0)
 else
  @NLparameter(n.mdl, w_goal_param == c.w.goal)
  @NLparameter(n.mdl, w_psi_param == c.w.psi)
 end
 obj_params=[w_goal_param,w_psi_param]

 # penalize distance to goal
 goal_obj=@NLexpression(n.mdl,w_goal_param*((x[end] - c.g.x_ref)^2 + (y[end] - c.g.y_ref)^2)/((x[1] - c.g.x_ref)^2 + (y[1] - c.g.y_ref)^2 + EP))

 # penalize difference between final heading angle and angle relative to the goal NOTE currently this is broken becasue atan2() is not available
 #psi_frg=@NLexpression(n.mdl,asin(c.g.y_ref-y[end])/(acos(c.g.x_ref-x[end]) + EP) )
 #psi_obj=@NLexpression(n.mdl,w_psi_param*(asin(sin(psi[end] - psi_frg))/(acos(cos(psi[end] - psi_frg)) + EP) )^2 )
 psi_obj=0;
 # psi_obj=@NLexpression(n.mdl,w_psi_param*(asin(sin(psi[end] - asin(c.g.y_ref-y[end])/acos(c.g.x_ref-x[end])))/(acos(cos(psi[end] - asin(c.g.y_ref-y[end])/acos(c.g.x_ref-x[end]))) + EP) )^2 )

 # soft constraints on vertical tire load
 tire_obj=integrate!(n,tire_expr);

 # minimizing the integral over the entire prediction horizon of the line that passes through the goal
 #haf_obj=integrate!(n,:( $c.w.haf*( sin($c.g.psi_ref)*(x[j]-$c.g.x_ref) - cos($c.g.psi_ref)*(y[j]-$c.g.y_ref) )^2 ) )
 haf_obj=0;
 # penalize control effort
 ce_obj=integrate!(n,:($c.w.ce*($c.w.sa*(sa[j]^2)+$c.w.sr*(sr[j]^2)+$c.w.jx*(jx[j]^2))) )

 @NLobjective(n.mdl, Min, goal_obj + psi_obj + c.w.Fz*tire_obj + haf_obj + c.w.time*n.tf + ce_obj )

 #########################
 # intial optimization (s)
 ########################
 n.s.save=false; n.s.MPC=false; n.s.evalConstraints=false;
 if n.s.save
  warn("saving initial optimization results where functions where cashed!")
 end
 for k in 1:1
  optimize!(n);
  if n.r.status==:Optimal; break; end
 end

 # modifify the maximum solver time
 SS=Dict((:name=>c.m.solver),(:mpc_defaults=>true),(solver_time=>c.m.max_cpu_time))
 defineSolver!(n,SS)

         #  1      2          3          4
 n.params=[pa,obs_params,LiDAR_params,obj_params];

 n.s.save=true; # settings

 #n.s.evalConstraints=true

 return n
end

"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/20/2017, Last Modified: 7/05/2017 \n
--------------------------------------------------------------------------------------\n
"""
function updateAutoParams!(n,c)

  # obstacle information-> only show if it is in range at the start TODO
  goal_in_range = goalRange!(n,c)
  if goal_in_range # TODO make a flag that indicates this switch has been flipped
    println("goal is in range")

    # enforce final state constraints on x and y position
    for st=1:2;
      setRHS(n.r.xf_con[st,1], +(n.XF[st]+n.XF_tol[st]));
      setRHS(n.r.xf_con[st,2], -(n.XF[st]-n.XF_tol[st]));
    end

    # relax LiDAR constraints
    setvalue(n.params[3][1], 1e6)
    setvalue(n.params[3][2],-1e6)

    # remove terms in cost function
    setvalue(n.params[4][1],0.0)
    setvalue(n.params[4][2],0.0)
  end
  #NOTE assuming it is not going in and out of range

 return goal_in_range
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 7/04/2017, Last Modified: 7/04/2017 \n
--------------------------------------------------------------------------------------\n
"""
function goalRange!(n,c)
 return ( (n.mpc.X0[end][1] - c.g.x_ref)^2 + (n.mpc.X0[end][2] - c.g.y_ref)^2 )^0.5 < c.m.Lr
end

end # module


text:
    A = [1 2 ; 3 4]
    B = [10, 10]
    sum!(B, A)
    sum!(A, B)

    using ReverseDiffSource

        Nn = [20,10,1]
        function ann(l1, l2, x1)
            x2 = l1 * x1
            x2 = max(0., x2) #relu
            x3 = l2 * x2
            x3 = (max(0., x3[1]))^2 #relu
        end

        # ls = Array[ nn.ws...]
        ls = init(Nn)
        dann = rdiff(ann, (ls..., ones(20) ))
        out = dann(ls..., ones(20) )

    nn = S.NN((ls[1], ls[2]), (S.Relu(), S.Relu()))
    ts2 = S.cycle(nn, ones(20), zeros(1))

        out[2]
        ts2.dws[1]

        out[3]
        ts2.dws[2]

        ts2.ds[end]

        out[3] ./ ts2.dws[2]

    r1 = S.calc(nn, ones(20))[1] # 1.16217
    r2 = ann(ls..., ones(20))
    r1-r2

    ls[1][9,1] += 1e-10
    (ann(ls..., ones(20)) - r1) * 1e10  # -0.8600

    nn.ws[1][9,1] += 1e-10
    (S.calc(nn, ones(20)) - r2) * 1e10  # -0.8600

    reload("SimpleNeuralNets"); S=SimpleNeuralNets

    nn = S.NN([10,10,10,10,1], [ S.Tanh(), S.Tanh()] )
    ts = S.TrainState(nn, 10)

    size(nn)
    size(ts)

    ts2 = S.cycle(nn, ones(20,10), ones(5,10))
    S.cycle!(ts, nn, ones(10,10), ones(1,10))

    S.depth(nn)

    s0 = S.calc(nn, rand(20))
    copy!(s0.ss[1], ones(20))
    S.compatible(nn, s0)
    S.calc!(nn, s0)

    s0 = S.calc(nn, rand(20,10))
    copy!(s0.ss[1], ones(20))
    S.compatible(nn, s0)
    S.calc!(nn, s0)

    a = (1,2,3)
    typeof(a)
    a[1]
    a[1] = 2

    nn = S.NN((randn(1,1),), (S.Relu(),))
    S.calc(nn, [1.])

    nn = S.NN((randn(1,1),), (S.SoftRelu(),))
    S.calc(nn, [1.])

########### ANN3 #################################"

    import Bokeh
    # using ReverseDiffSource
    reload("SimpleNeuralNets"); s = SimpleNeuralNets

    f(x) = exp(-(x-1)^2/5)*cos(x)*sign(x)

    xs = collect(-5:0.001:5)
    X = Float64[ x for x in xs ]'
    Y = Float64[ f(x) for x in xs ]'

    nn = s.NN([1,10,10,10,1], [ s.Relu(), s.Tanh()] )

    ts = s.TrainingSet(X, Y)
    any(ts.weights .!= 1.0)

    λ0, μ = 1e-3, 0.9
    s.sgd(nn, ts, λ0, μ, freq=200, steps=1000)

    ts2 = s.TrainingSet(X, Y, SqEuclidean(), rand(size(X,2)))
    s.sgd(nn, ts2, λ0, μ, freq=200, steps=1000)

    s.colwise(a,b,c)
    length(a.weights)

    px = s.calc(nn, X)
    sum(s.colwise(ts2.metric, px, ts2.Y))

    λ0, μ = 1e-3, 0.9
    s.sgd(nn, λ0, μ, X, Y, freq=200, steps=10000)
    λ0, μ = 1e-5, 0.999
    s.sgd(nn, λ0, μ, X, Y, freq=200, steps=10000)

    Bokeh.hold(true)
    Bokeh.plot(f, -5:5, "k-")
    px = collect(-5:0.01:5)
    py = s.calc(nn, px')'
    Bokeh.plot(px, py, "b-")
    Bokeh.showplot("/tmp/plot.html")
    Bokeh.hold(false)

    nn = s.NN([1,10,10,10,1], [ s.SoftRelu(), s.Tanh()] )
    λ0, μ = 1e-3, 0.9
    s.sgd(nn, λ0, μ, X, Y, freq=200, steps=10000)
    λ0, μ = 1e-3, 0.
    s.sgd(nn, λ0, μ, X, Y, freq=1000, steps=10000)

    Bokeh.hold(true)
    Bokeh.plot(f, -5:5, "k-")
    px = collect(-4:0.1:4)
    py = s.calc(nn, px')'
    Bokeh.plot(px, py, "b-")
    Bokeh.showplot("/tmp/plot.html")
    Bokeh.hold(false)

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
@enum PeelSetting peelf peellen
include("./peelAndBoundSubroutines.jl")
include("./peelAndBound.jl")
include("./parallelPAB.jl")
include("./branchAndBound.jl")


text:
module NormalizeQuantiles

export normalizeQuantiles

export sampleRanks
export qnTiesMethods,tmMin,tmMax,tmOrder,tmReverse,tmRandom,tmAverage

using Distributed
using SharedArrays
using Random
using Statistics

@enum qnTiesMethods tmMin tmMax tmOrder tmReverse tmRandom tmAverage

function checkForNotANumber(x::Any)
    (!isa(x,Integer) && !isa(x,Real)) || isnan(x)
end

@doc "
### qnmatrix::Array{Float64} function normalizeQuantiles(matrix::AbstractArray)
Calculate the quantile normalized data for the input matrix

Parameter:
    matrix::AbstractArray
The input data as an array of values interpreted as matrix(rows,columns)

Return value:
    qnmatrix::Array{Float64}
The quantile normalized data as Array{Float64}

Example:

    using NormalizeQuantiles

    array = [ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ]
    qn = normalizeQuantiles(array)

    row = 2
    column = 2
    array=convert(Array{Any},array)
    array[row,column] = missing
    qn = normalizeQuantiles(array)
"
function normalizeQuantiles(matrix::AbstractArray)
    if ndims(matrix) > 2
        throw(ArgumentError("normalizeQuantiles expects an array of dimension 2"))
    end
    nrows=size(matrix,1)
    ncols=size(matrix,2)
    # preparing the result matrix
    qnmatrix=SharedArray{Float64}(nrows,ncols)
    if ncols>0 && nrows>0
        # foreach column: sort the values without NAs; put NAs (if any) back into sorted list
        NormalizeQuantiles.sortColumns!(matrix,qnmatrix,nrows,ncols)
        # foreach row: set all values to the mean of the row, except NAs
        NormalizeQuantiles.meanRows!(qnmatrix,nrows)
        # foreach column: equal values in original column should all be mean of normalized values
        # foreach column: reorder the values back to the original order
        NormalizeQuantiles.equalValuesInColumnAndOrderToOriginal!(matrix,qnmatrix,nrows,ncols)
    end
    convert(Array{Float64},qnmatrix)
end

function sortColumns!(matrix::AbstractArray,qnmatrix::SharedArray{Float64},nrows,ncols)
    @inbounds @sync @distributed for column = 1:ncols
        sortcol=matrix[:,column]
        goodIndices=[ !NormalizeQuantiles.checkForNotANumber(x) for x in sortcol ]
        missingIndices=(1:nrows)[.!goodIndices]
        sortcol[goodIndices]=Array{Float64}(sortcol[goodIndices])
        length(missingIndices)>0 ? sortcol[missingIndices].=NaN : nothing
        sort!(sortcol)
        for missingPos in missingIndices
            sortcol[(missingPos+1):end]=sortcol[missingPos:(end-1)]
            sortcol[missingPos]=NaN
        end
        qnmatrix[:,column]=sortcol
    end
end

function meanRows!(qnmatrix::SharedArray{Float64},nrows)
    @inbounds @sync @distributed for row = 1:nrows
        goodIndices=[ ! NormalizeQuantiles.checkForNotANumber(x) for x in qnmatrix[row,:] ]
        rowmean=mean(qnmatrix[row,goodIndices])
        qnmatrix[row,goodIndices].=rowmean
    end
end

function equalValuesInColumnAndOrderToOriginal!(matrix::AbstractArray,qnmatrix::SharedArray{Float64},nrows,ncols)
    qncol=Array{Float64}(undef,nrows,1)
    @inbounds @sync @distributed for column = 1:ncols
        goodIndices=[ !NormalizeQuantiles.checkForNotANumber(x) for x in vec(matrix[:,column]) ]
        sortp=sortperm([ Float64(x) for x in vec(matrix[goodIndices,column]) ])
        goodIndices2=[ !NormalizeQuantiles.checkForNotANumber(x) for x in vec(qnmatrix[:,column]) ]
        if length(sortp)>0
            NormalizeQuantiles.setMeanForEqualOrigValues(matrix[goodIndices,column][sortp],qnmatrix,column,goodIndices)
        end
        fill!(qncol,NaN)
        qncol[(1:nrows)[goodIndices][sortp]]=vec(qnmatrix[(1:nrows)[goodIndices2],column])
        qnmatrix[:,column]=qncol
    end
end

function setMeanForEqualOrigValues(sortedArrayNoNAs::AbstractArray,qnmatrix::SharedArray{Float64},column::Int,goodIndices::Array{Bool})
    nrows=length(sortedArrayNoNAs)
    foundIndices=zeros(Int,nrows)
    goodIndices2=(1:length(goodIndices))[goodIndices]
    count=1
    lastValue=sortedArrayNoNAs[1]
    for i in 2:nrows
        nextValue=sortedArrayNoNAs[i]
        if nextValue==lastValue
            if count==1
                foundIndices[count]=goodIndices2[i-1]
                count+=1
                foundIndices[count]=goodIndices2[i]
            else
                foundIndices[count]=goodIndices2[i]
            end
            count+=1
        else
            if count>1
                qnmatrix[foundIndices[1:(count-1)],column].=mean([ x for x in qnmatrix[foundIndices[1:(count-1)],column] ])
            end
            count=1
        end
        lastValue=nextValue
    end
    if count>1
        qnmatrix[foundIndices[1:(count-1)],column].=mean([ x for x in qnmatrix[foundIndices[1:(count-1)],column] ])
    end
end

@doc "
### (Array{Union{Missing,Int}},Dict{Int,Array{Int}}) sampleRanks(array::AbstractArray;tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)

Calculate ranks of the values of a given vector.

Parameters:

    array: the input array
    tiesMethod: the method how ties (equal values) are treated
       possible values: tmMin tmMax tmOrder tmReverse tmRandom tmAverage
       default is tmMin
    naIncreasesRank: if true than any NA increases the following ranks by 1
    resultMatrix: if true than return a dictionary of rank keys and array of indices values

Example:

    using NormalizeQuantiles

    a = [ 5.0 2.0 4.0 3.0 1.0 ]

    (r,m)=sampleRanks(a)

    (r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=false,resultMatrix=true)

r is the vector of ranks.

m is a dictionary with rank as keys and as value the indices of all values of this rank.

"
function sampleRanks(array::AbstractArray;tiesMethod::qnTiesMethods=tmMin,naIncreasesRank::Bool=false,resultMatrix::Bool=false)
    nrows=length(array)
    goodIndices=falses(nrows)
    naCounts=zeros(Int,nrows)
    naCount=0
    reducedIndex=1
    for arrayIndex in 1:length(array)
        if NormalizeQuantiles.checkForNotANumber(array[arrayIndex])
            naCount+=1
        else
            goodIndices[arrayIndex]=true
            naCounts[reducedIndex]=naCount
            reducedIndex+=1
        end
    end
    reducedArraySorted=[ Float64(x) for x in array[goodIndices] ]
    reducedArraySortedIndices=sortperm(reducedArraySorted)
    reducedArraySorted=reducedArraySorted[reducedArraySortedIndices]
    result=Array{Union{Missing,Int}}(missing,nrows)
    rankMatrix=Dict{Int,Array{Int}}()
    group = Array{Int,1}()
    firstFound = true
    lastFound = undef
    nextRank=1
    rankIncrement=0
    doIncreaseRank=false
    headingNA=false
    groupIndex=1
    for resultIndex in 1:nrows
        if goodIndices[resultIndex]
            if headingNA
                headingNA=false
                if doIncreaseRank
                    nextRank+=rankIncrement
                    rankIncrement=0
                    doIncreaseRank=false
                end
            end
            reducedArrayIndex=reducedArraySortedIndices[groupIndex]
            if !firstFound && lastFound != reducedArraySorted[groupIndex]
                nextRank = NormalizeQuantiles.setRank!(group,nextRank,tiesMethod,result,resultMatrix,rankMatrix)
                group = [naCounts[reducedArrayIndex]+reducedArrayIndex]
                if doIncreaseRank
                    nextRank+=rankIncrement
                    rankIncrement=0
                    doIncreaseRank=false
                end
            else
                push!(group,naCounts[reducedArrayIndex]+reducedArrayIndex)
            end
            lastFound = reducedArraySorted[groupIndex]
            firstFound = false
            groupIndex+=1
        else
            if naIncreasesRank
                doIncreaseRank=true
                rankIncrement+=1
            end
            if firstFound
                headingNA=true
            end
        end
    end
    NormalizeQuantiles.setRank!(group,nextRank,tiesMethod,result,resultMatrix,rankMatrix)
    (result,rankMatrix)
end

function setRank!(group::Array{Int},nextRank::Int,tiesMethod::qnTiesMethods,result::Array{Union{Missing,Int}},resultMatrix::Bool,rankMatrix::Dict{Int,Array{Int}})
    ranksCount=length(group)
    ranks=nextRank:(nextRank+ranksCount-1)
    if tiesMethod==tmMin
        minRank=minimum(ranks)
        result[group].=minRank
        if resultMatrix
            rankMatrix[minRank]=group
        end
        nextRank+=1
    elseif tiesMethod==tmMax
        maxRank=maximum(ranks)
        result[group].=maxRank
        if resultMatrix
            rankMatrix[maxRank]=group
        end
        nextRank=maxRank+1
    elseif tiesMethod==tmOrder
        maxRank=maximum(ranks)
        result[group]=ranks
        if resultMatrix
            for rankIndex in 1:length(ranks)
                rankMatrix[ranks[rankIndex]]=[group[rankIndex]]
            end
        end
        nextRank=maxRank+1
    elseif tiesMethod==tmReverse
        maxRank=maximum(ranks)
        ranks=reverse(ranks,dims=1)
        result[group]=ranks
        if resultMatrix
            for rankIndex in 1:length(ranks)
                rankMatrix[ranks[rankIndex]]=[group[rankIndex]]
            end
        end
        nextRank=maxRank+1
    elseif tiesMethod==tmRandom
        maxRank=maximum(ranks)
        ranks=ranks[randperm(ranksCount)]
        result[group]=ranks
        if resultMatrix
            for rankIndex in 1:length(ranks)
                rankMatrix[ranks[rankIndex]]=[group[rankIndex]]
            end
        end
        nextRank=maxRank+1
    elseif tiesMethod==tmAverage
        rankMean=round(Int,mean(ranks))
        result[group].=rankMean
        if resultMatrix
            rankMatrix[rankMean]=group
        end
        nextRank=rankMean+1
    end
    nextRank
end

end # module


text:
module Uncomprehensions

using MLStyle
using MLStyle.MatchCore: mangle
using MLStyle.Infras: @format

def_pattern(Uncomprehensions,
    predicate = (x -> @match x begin
        :[$_ for $_ in $_ if $_] ||
        :[$_ for $_ in $_]    => true
        _                     => false
    end),
    rewrite = (tag, case, mod) -> begin
        @match case begin
            :[$expr for $iter in $seq if $cond]  ||
            :[$expr for $iter in $seq] && Do(cond=nothing) => begin
                function (body)
                    iter_var = mangle(mod)

                    # the name of element that'll be pushed to the result vector
                    produced_elt_var = mangle(mod)

                    inner_test_var = mangle(mod)
                    seq_var = mangle(mod)

                    # pattern is `seq`, return expr is `body`
                    decons_seq = mk_pattern(seq_var, seq, mod)(body)
                    lnode = LineNumberNode(1, "Uncomprehensions.jl")
                    if cond === nothing
                        # pattern is `expr`, return expr is `iter`
                        decons_elt = mk_pattern(iter_var, expr, mod)(iter)
                        @format [
                            seq_var, iter_var, inner_test_var, produced_elt_var,
                            body, tag, decons_elt, decons_seq, push!, lnode
                        ] quote
                            if tag isa Vector
                                let seq_var = [], inner_test_var = true

                                    for iter_var in tag
                                        produced_elt_var = decons_elt
                                        if produced_elt_var === failed
                                            inner_test_var = false
                                            break
                                        end
                                        push!(seq_var, produced_elt_var)
                                    end

                                    if inner_test_var
                                        decons_seq
                                    else
                                        failed
                                    end
                                end
                            else
                                failed
                            end
                        end
                    else
                        mid_var = mangle(mod)
                        decons_elt = mk_pattern(iter_var, expr, mod)(
                            @format [cond, iter, mid_var] quote
                                mid_var = iter
                                mid_var === failed ? failed : (cond, mid_var)
                            end
                        )
                        @format [
                            seq_var, iter_var, inner_test_var, produced_elt_var,
                            body, tag, decons_elt, decons_seq, push!, lnode
                        ] quote
                            if tag isa Vector
                                let seq_var = [], inner_test_var = true

                                    for iter_var in tag
                                        produced_elt_var = decons_elt
                                        if produced_elt_var === failed
                                            inner_test_var = false
                                            break
                                        end
                                        if produced_elt_var[1]
                                            push!(seq_var, produced_elt_var[2])
                                        end
                                    end

                                    if inner_test_var
                                        decons_seq
                                    else
                                        failed
                                    end
                                end
                            else
                                failed
                            end
                        end
                    end
                end
            end
        end
    end
)
end


text:
using Distributions
using Test
using KernelGoodnessOfFit

@testset "Univeriate Normal" begin
    xs = randn(1, 100);
    q = Normal(1.0, 1.0)

    t1 = FSSDrand(xs, q)

    println(t1)

    @test pvalue(t1) ≥ 0.0
end


text:
@ms using LinearAlgebra

"""
Instance data structure definition
"""
mutable struct Instance

    name::String # name of the instance

    n::Int # number of airdromes
    I::Int # departure airdrome
    F::Int # arrival airdrome
    k::Int # minimum number of airdromes to be visited
    m::Int # total number of groups
    R::Array{Int} # association airdrome-group
    coordinates::Array{Tuple{Int, Int}} # geographical coordinates of the airdromes
    D::Dict{Tuple{Int, Int}, Int} # rounded distance between the airdromes
    Δ::Int # maximum distance that can be traveled without landing to refuel
    feasible_arcs::Set{Tuple{Int, Int}} # feasible arcs regarding the threshold distance Δ

    """
    Instance object builder
    """
    function Instance()

        this = new()

        # initialization of the containers which will be filled dynamically
        this.R = Int[]
        this.coordinates = Tuple{Int, Int}[]
        this.D = Dict{Tuple{Int, Int}, Int}()
        this.feasible_arcs = Set{Tuple{Int, Int}}()

        # read data from the input file
        read_input_file!(this, args["input"])

        # compute integer distances using from coordinates and update the feasible arcs set
        compute_arc_weights!(this)

        return this

    end

end

"""
Input file reader
"""
function read_input_file!(inst::Instance, path::String)

    # the name of the instance is the name of the file
    base, ext = splitext(basename(path))
    inst.name = base

    # collect lines
    lines = readlines(path)

    # read data

    inst.n = parse(Int, lines[1])
    inst.I = parse(Int, lines[2])
    inst.F = parse(Int, lines[3])
    inst.k = parse(Int, lines[4])
    inst.m = parse(Int, lines[5])

    for r in split(lines[7])
        push!(inst.R, parse(Int, r))
    end

    inst.Δ = parse(Int, lines[9])

    for i in 11:10+inst.n
        x_str, y_str = split(lines[i])
        push!(inst.coordinates, (parse(Int, x_str), parse(Int, y_str)))
    end

end

"""
Computes the weight of each arc and concludes on their feasibility
"""
function compute_arc_weights!(inst::Instance)

    for i in 1:inst.n
        for j in i+1:inst.n # we use the symmetry of the distance

            # compute distance using LinearAlgebra
            inst.D[(i, j)] = round(Int, norm(inst.coordinates[i].-inst.coordinates[j]))
            inst.D[(j, i)] = inst.D[(i, j)]

            # conclude on the feasibility of the arc
            if inst.D[(i, j)] <= inst.Δ
                push!(inst.feasible_arcs, (i,j), (j,i))
            end

        end
    end

end

max_stars_repo_name: antcrepin/flight-route-optimizer

max_issues_repo_name: antcrepin/flight-route-optimizer

max_forks_repo_name: antcrepin/flight-route-optimizer


text:
module RememberToUpdateRegistryCI

using GitHub
using Pkg

# Some of the code in this file is taken from:
# CompatHelper.jl (https://github.com/bcbi/CompatHelper.jl)

struct AlwaysAssertionError <: Exception
end

@inline function always_assert(cond::Bool)::Nothing
    cond || throw(AlwaysAssertionError())
    return nothing
end

function get_all_pull_requests(repo::GitHub.Repo,
                               state::String;
                               auth::GitHub.Authorization,
                               per_page::Integer = 100,
                               page_limit::Integer = 100)
    all_pull_requests = Vector{GitHub.PullRequest}(undef, 0)
    myparams = Dict("state" => state,
                    "per_page" => per_page,
                    "page" => 1)
    prs, page_data = GitHub.pull_requests(repo;
                                          auth=auth,
                                          params = myparams,
                                          page_limit = page_limit)
    append!(all_pull_requests, prs)
    while haskey(page_data, "next")
        prs, page_data = GitHub.pull_requests(repo;
                                              auth=auth,
                                              page_limit = page_limit,
                                              start_page = page_data["next"])
        append!(all_pull_requests, prs)
    end
    unique!(all_pull_requests)
    return all_pull_requests
end

_repos_are_the_same(::GitHub.Repo, ::Nothing) = false
_repos_are_the_same(::Nothing, ::GitHub.Repo) = false
_repos_are_the_same(::Nothing, ::Nothing) = false
function _repos_are_the_same(x::GitHub.Repo, y::GitHub.Repo)
    if x.name == y.name && x.full_name == y.full_name &&
                           x.owner == y.owner &&
                           x.id == y.id &&
                           x.url == y.url &&
                           x.html_url == y.html_url &&
                           x.fork == y.fork
       return true
    else
        return false
    end
end

function exclude_pull_requests_from_forks(repo::GitHub.Repo, pr_list::Vector{GitHub.PullRequest})
    non_forked_pull_requests = Vector{GitHub.PullRequest}(undef, 0)
    for pr in pr_list
        always_assert(_repos_are_the_same(repo, pr.base.repo))
        if _repos_are_the_same(repo, pr.head.repo)
            push!(non_forked_pull_requests, pr)
        end
    end
    return non_forked_pull_requests
end

function only_my_pull_requests(pr_list::Vector{GitHub.PullRequest}; my_username::String)
    _my_username_lowercase = lowercase(strip(my_username))
    n = length(pr_list)
    pr_is_mine = BitVector(undef, n)
    for i = 1:n
        pr_user_login = pr_list[i].user.login
        if lowercase(strip(pr_user_login)) == _my_username_lowercase
            pr_is_mine[i] = true
        else
            pr_is_mine[i] = false
        end
    end
    my_pr_list = pr_list[pr_is_mine]
    return my_pr_list
end

function create_new_pull_request(repo::GitHub.Repo;
                                 base_branch::String,
                                 head_branch::String,
                                 title::String,
                                 body::String,
                                 auth::GitHub.Authorization)
    params = Dict{String, String}()
    params["title"] = title
    params["head"] = head_branch
    params["base"] = base_branch
    params["body"] = body
    result = GitHub.create_pull_request(repo; params = params, auth = auth)
    return result
end

function git_commit(message)::Bool
    return try
        success(`git commit -m "$(message)"`)
    catch
        false
    end
end

function generate_username_mentions(usernames::AbstractVector)::String
    intermediate_result = ""
    for username in usernames
        _username = filter(x -> x != '@', strip(username))
        if length(_username) > 0
            intermediate_result = intermediate_result * "\ncc: @$(_username)"
        end
    end
    final_result = convert(String, strip(intermediate_result))
    return final_result
end

function set_git_identity(username, email)
    run(`git config user.name "$(username)"`)
    run(`git config user.email "$(email)"`)
    return nothing
end

function create_new_pull_request(repo::GitHub.Repo;
                                 base_branch::String,
                                 head_branch::String,
                                 title::String,
                                 body::String,
                                 auth::GitHub.Authorization)
    params = Dict{String, String}()
    params["title"] = title
    params["head"] = head_branch
    params["base"] = base_branch
    params["body"] = body
    result = GitHub.create_pull_request(repo; params = params, auth = auth)
    return result
end

function main(relative_path;
              registry,
              github_token = ENV["GITHUB_TOKEN"],
              master_branch = "master",
              pr_branch = "github_actions/remember_to_update_registryci-$(Base.VERSION.major).$(Base.VERSION.minor).$(Base.VERSION.patch)",
              pr_title = "Update RegistryCI.jl by updating the .ci/Manifest.toml file ($(Base.VERSION.major).$(Base.VERSION.minor).$(Base.VERSION.patch))",
              cc_usernames = String[],
              # my_username = ,
              my_email = "41898282+github-actions[bot]@users.noreply.github.com",
              old_julia_version::VersionNumber)
    original_project = Base.active_project()
    original_directory = pwd()

    OLDMAJOR = old_julia_version.major
    OLDMINOR = old_julia_version.minor
    OLDPATCH = old_julia_version.patch

    @info "" old_julia_version OLDMAJOR OLDMINOR OLDPATCH

    tmp_dir = mktempdir()
    atexit(() -> rm(tmp_dir; force = true, recursive = true))
    cd(tmp_dir)

    auth = GitHub.authenticate(github_token)
    my_username = try
        GitHub.whoami(; auth = auth).login
    catch
        "github-actions[bot]"
    end
    my_repo = GitHub.repo(registry; auth = auth)
    registry_url_with_auth = "https://x-access-token:$(github_token)@github.com/$(registry)"
    _all_open_prs = get_all_pull_requests(my_repo, "open"; auth = auth)
    _nonforked_prs = exclude_pull_requests_from_forks(my_repo, _all_open_prs)
    pr_list = only_my_pull_requests(_nonforked_prs; my_username = my_username)
    pr_titles = Vector{String}(undef, length(pr_list))
    for i = 1:length(pr_list)
        pr_titles[i] = convert(String, strip(pr_list[i].title))::String
    end

    username_mentions_text = generate_username_mentions(cc_usernames)

    run(`git clone $(registry_url_with_auth) REGISTRY`)
    cd("REGISTRY")
    run(`git checkout $(master_branch)`)
    run(`git checkout -B $(pr_branch)`)
    cd(relative_path)
    main_manifest = "Manifest.toml"
    old_manifest = "Manifest.$(OLDMAJOR).$(OLDMINOR).$(OLDPATCH).toml"
    main_manifest_filename = joinpath(pwd(), main_manifest)
    old_manifest_filename = joinpath(pwd(), old_manifest)
    _backup_filename = joinpath(pwd(), "Manifest.toml.original")
    is_old_julia_version = (Base.VERSION.major == OLDMAJOR) && (Base.VERSION.minor == OLDMINOR) && (Base.VERSION.patch == OLDPATCH)
    if is_old_julia_version
        mv(main_manifest, _backup_filename; force = true)
        rm(old_manifest; force = true)
    else
        rm(main_manifest; force = true)
    end
    Pkg.activate(pwd())
    Pkg.instantiate()
    Pkg.update()
    if is_old_julia_version
        mv(main_manifest, old_manifest; force = true)
        mv(_backup_filename, main_manifest; force = true)
    else
    end
    set_git_identity(my_username, my_email)
    try
        if is_old_julia_version
            run(`git add $(old_manifest)`)
        else
            run(`git add $(main_manifest)`)
        end
    catch
    end
    commit_was_success = git_commit("Update .ci/Manifest.toml")
    @info("commit_was_success: $(commit_was_success)")
    if commit_was_success
        run(`git push -f origin $(pr_branch)`)
        if pr_title in pr_titles
            @info("An open PR with the title already exists", pr_title)
        else
            new_pr_body = strip(string("This pull request updates ",
                                       "RegistryCI.jl by updating the ",
                                       "`.ci/Manifest.toml` file.\n\n",
                                       username_mentions_text))
            _new_pr_body = convert(String, strip(new_pr_body))
            create_new_pull_request(my_repo;
                                    base_branch = master_branch,
                                    head_branch = pr_branch,
                                    title = pr_title,
                                    body = _new_pr_body,
                                    auth = auth)
        end
    end

    cd(original_directory)
    rm(tmp_dir; force = true, recursive = true)
    Pkg.activate(original_project)
    return commit_was_success
end

end # end module RememberToUpdateRegistryCI

max_stars_repo_name: kdyrhage/General

max_issues_repo_name: kdyrhage/General

max_forks_repo_name: kdyrhage/General


text:
using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()

using Documenter, BigBed

makedocs(
    format = Documenter.HTML(
        edit_branch = "develop"
    ),
    sitename = "BigBed.jl",
    pages = [
        "Home" => "index.md",
    ],
    authors = "Kenta Sato, D. C. Jones, Ben J. Ward, The BioJulia Organisation and other contributors."
)
deploydocs(
    repo = "github.com/BioJulia/BigBed.jl.git",
    devbranch = "develop"
)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# types used to select a column
const ColSelector = Union{Symbol,Integer,AbstractString}

"""
    ColSpec

`ColSpec` is a union of types used to filter columns.
The `ColSpec` type together with the `ColSelector` union type and
the `choose` function form the ColSpec interface.

To implement the ColSpec interface, the following steps must be performed:

1. add colspec fied:
```julia
struct MyTransform{S<:ColSpec,#= other type params =#}
  colspec::S
  # other fileds
end
```
2. use `choose(colspec, names)` function in apply:
```julia
function apply(transform::MyTransform, table)
  names = Tables.columnnames(table)
  # selected column names
  snames = choose(transform.colspec, names)
  # code...
end
```

If you need to create constructors that accept
individual column selectors use the `ColSelector` type. Example:
```julia
function MyTransform(args::T...) where {T<:ColSelector}
  # code...
end
```
"""
const ColSpec = Union{Vector{T},NTuple{N,T},Regex,Colon} where {N,T<:ColSelector}

"""
    choose(colspec::ColSpec, names) -> Vector{Symbol}

Choose column `names` using `colspec`.

# Examples

```julia
julia> names = (:a, :b, :c, :d, :e, :f);

julia> choose(["a", "c", "e"], names)
3-element Vector{Symbol}:
 :a
 :c
 :e

julia> choose((1, 3, 5), names)
3-element Vector{Symbol}:
 :a
 :c
 :e

julia> choose(r"[ace]", names)
3-element Vector{Symbol}:
 :a
 :c
 :e
```
"""
function choose(colspec::Vector{Symbol}, names)
  # validate columns
  @assert !isempty(colspec) "Invalid column selection."
  @assert colspec ⊆ names "Invalid column selection."
  return colspec
end

choose(colspec::Vector{<:AbstractString}, names) =
  choose(Symbol.(colspec), names)

choose(colspec::Vector{<:Integer}, names::Vector) =
  choose(names[colspec], names)

choose(colspec::Vector{<:Integer}, names::Tuple) =
  choose(colspec, collect(names))

choose(colspec::NTuple{N,<:ColSelector}, names) where {N} =
  choose(collect(colspec), names)

function choose(colspec::Regex, names::Vector)
  fnames = filter(n -> occursin(colspec, String(n)), names)
  choose(fnames, names)
end

choose(colspec::Regex, names::Tuple) =
  choose(colspec, collect(names))

choose(::Colon, names::Vector) = names
choose(::Colon, names::Tuple) = collect(names)


text:
@inline function (-){G<:Grasp,R<:Real}(a::Rvl{G,R})
   Rvl{negate(G), R}(-a.hi, -a.lo)
end

function (-){G1<:Grasp,G2<:Grasp,R<:Real}(a::Rvl{G1,R}, b::Rvl{G2,R})
    aLoIsOpen, aHiIsOpen = boundaries(G1)
    bLoIsOpen, bHiIsOpen = boundaries(G2)

    abGrasp = boundaries( (aLoIsOpen|bHiIsOpen), (aHiIsOpen|bLoIsOpen) )

    lo = (-)(a.lo, b.hi, RoundDown)
    hi = (-)(a.hi, b.lo, RoundUp)

    Rvl{abGrasp,R}(lo, hi)
end

(-){G<:Grasp,R<:Real}(a::Rvl{G,R}, b::R) = (-)(a, Rvl{G,R}(b))
(-){G<:Grasp,R<:Real}(a::R, b::Rvl{G,R}) = (-)(Rvl{G,R}(a), b)

(-){G<:Grasp,I<:WorkInt}(a::Rvl{G,Float64}, b::I) = (-)(a, Rvl{G,Float64}(convert(Float64,b)))
(-){G<:Grasp,I<:WorkInt}(a::I, b::Rvl{G,Float64}) = (-)(Rvl{G,Float64}(convert(Float64,a)), b)
(-){G<:Grasp,I<:WorkInt}(a::Rvl{G,Float32}, b::I) = (-)(a, Rvl{G,Float32}(convert(Float32,b)))
(-){G<:Grasp,I<:WorkInt}(a::I, b::Rvl{G,Float32}) = (-)(Rvl{G,Float32}(convert(Float32,a)), b)


text:
module AppendedArraysTests

using Test
using Gridap.Arrays

a = collect(11:20)
b = collect(Float64,101:120)
c = lazy_append(a,b)
r = vcat(a,b)
test_array(c,r)

a = collect(Float64,11:20)
b = collect(Float64,101:120)
c = lazy_append(a,b)
r = vcat(a,b)
test_array(c,r)

@test sum(c) == sum(r)

d = apply(-,c)
r = -c
test_array(d,r)
@test isa(d,AppendedArray)

e = apply(-,c,d)
r = c-d
test_array(e,r)
@test isa(e,AppendedArray)

d = apply(Float64,-,c)
r = -c
test_array(d,r)
@test isa(d,AppendedArray)

e = apply(Float64,-,c,d)
r = c-d
test_array(e,r)
@test isa(e,AppendedArray)

end # module


text:
using MPI
using MPIArrays
using Test
using Random

MPI.Init()
const comm = MPI.COMM_WORLD
const rank = MPI.Comm_rank(comm)
const nb_procs = MPI.Comm_size(comm)

if rank == 0
@testset "Partitioning" begin
matpart = MPIArrays.ContinuousPartitioning([2,5,3], [2,3])
@test matpart[1,1] == (1:2,1:2)
@test matpart[3,2] == (8:10,3:5)
@test MPIArrays.partition_sizes(matpart) == ([2,5,3], [2,3])

vecpart = MPIArrays.ContinuousPartitioning([2,5,3])
@test vecpart[1] == (1:2,)
@test vecpart[2] == (3:7,)
@test vecpart[3] == (8:10,)
@test MPIArrays.partition_sizes(vecpart) == ([2,5,3],)

end

end

@testset "ArrayInterface" begin

function do_test(partitions::Vararg{Int,N}) where N
    sizes = (partitions .* 4) .+ (partitions .% (1:N))
    serial_A = reshape(1:prod(sizes), sizes...)
    A = MPIArray{Int}(comm, partitions, sizes...)
    forlocalpart!(lp -> lp .= serial_A[localindices(A, rank)...], A)
    @test all(A .== serial_A)
    sync(A)
    if rank == 1
        A .= serial_A
    end
    sync(A)
    @test all(A .== serial_A)
    return A
end

vec = do_test(nb_procs)
@test size(vec) == (nb_procs*4,)

if nb_procs % 2 == 0
    p1 = nb_procs ÷ 2
    p2 = nb_procs ÷ p1
else
    p1 = nb_procs
    p2 = 1
end
@test p1*p2 == nb_procs
mat = do_test(p1, p2)

ten = MPIArray{Float64}(comm, (1,1,nb_procs), 5,6,7)
sync(ten)
@test size(ten.partitioning) == (1,1,nb_procs)

forlocalpart!(lp -> lp .= rank+1, vec)
forlocalpart!(lp -> lp .= rank+1, mat)
MPI.Win_sync(mat.win)
MPI.Win_sync(vec.win)

sync(mat)

ab1 = MPIArrays.Block(mat, 3:size(mat,1), 3:5)
ab2 = MPIArrays.Block(mat, 1:4, 1:3)
ab3 = MPIArrays.Block(mat, 1:2, 2:5)
vb = MPIArrays.Block(vec, 3:length(vec))

if rank == 0
    println("vector:")
    display(vec)
    println()
    println("matrix:")
    display(mat)
    println()
end

sync(mat)

mat_redist = redistribute(mat, 1, nb_procs)

sync(mat_redist)

if rank == 0
    @test all(mat_redist .== mat)
end

sim1 = similar(mat, Float64)
@test sim1.partitioning == mat.partitioning
@test typeof(sim1[1]) == Float64
@test size(sim1) == size(mat)
sim2 = similar(mat, (10,))
@test size(sim2) == (10,)
sim3 = similar(mat, (3,4,5))
@test size(sim3) == (3,4,5)

sync(sim3)

# test manual distribution
distrib_vec = MPIArray{Int}(MPI.COMM_WORLD, fill(5, nb_procs), [10])
@test size(distrib_vec) == (5*nb_procs, 10)

end

@testset "ArrayBlock" begin

dims = MPIArray{Int}(comm, (nb_procs,), 2,)

if rank == 0
    dims[1] = rand(nb_procs:5*nb_procs+1)
    dims[2] = rand(nb_procs+1:5*nb_procs)
end

sync(dims)

nrows = dims[1]
ncols = dims[2]

if rank == 0
    println("Running ArrayBlock test with size $nrows×$ncols")
end

mat1 = MPIArray{Float64}(comm, (nb_procs,1), nrows, ncols)
forlocalpart!(rand!, mat1)
sync(mat1)

localblock = mat1[localindices(mat1)...]
localmat = allocate(localblock)
getblock!(localmat, localblock)
forlocalpart(mat1) do A
    @test(all(A .== localmat))
end

gb = GlobalBlock(localmat, localblock)

for (I,J) in zip(eachindex(gb), eachindex(localmat))
    @test gb[I] == localmat[J]
end

ghosted = GhostedBlock(mat1)
for I in eachindex(mat1)
    push!(ghosted, I)
end
sort!(ghosted)
sync(ghosted)
for I in eachindex(mat1)
    @test mat1[I] == getglobal(ghosted,I)
end

serial_mat2 = reshape(1:nrows*ncols,nrows,ncols)
mat2 = MPIArray{Int}(comm, (1,nb_procs), nrows, ncols)
forlocalpart!(mat2) do A
    A .= serial_mat2[localindices(mat2)...]
end

sync(mat2)

fullblock = mat2[1:nrows, 1:ncols]
fullmat = allocate(fullblock)
getblock!(fullmat, fullblock)
@test all(fullmat .== serial_mat2)

sync(mat2)

forlocalpart!(mat2) do A
    fill!(A,0)
end

sync(mat2)
putblock!(fullmat, fullblock)
sync(mat2)

forlocalpart!(mat2) do A
    @test all(A .== serial_mat2[localindices(mat2)...])
end

sync(mat2)
putblock!(fullmat, fullblock, +)
sync(mat2)

forlocalpart!(mat2) do A
    @test all(A .== serial_mat2[localindices(mat2)...] .* (nb_procs+1))
end

end

@testset "Filter" begin

localarray = fill(rank+1, 2*(rank+1))
vec_from_local = MPIArray(localarray, nb_procs)
sync(vec_from_local)
@test size(vec_from_local) == (sum(1:nb_procs)*2,)
@test forlocalpart(lp -> all(lp .== localarray), vec_from_local)

localarray = fill(rank+1, 5, 2*(rank+1))
mat_from_local = MPIArray(localarray, 1, nb_procs)
@test size(mat_from_local) == (5, sum(1:nb_procs)*2,)
@test forlocalpart(lp -> all(lp .== localarray), mat_from_local)

testrange = collect(1:100)
mpirange = MPIArray{Int}(length(testrange))
if rank == 0
    putblock!(testrange,mpirange[:])
end
sync(mpirange)

filtered_range = filter(isodd, mpirange)
@test length(filtered_range) == length(testrange) ÷ 2
@test all(filtered_range .== filter(isodd, testrange))

localstart_before_filter = localindices(mpirange)[1][1]

sync(mpirange)

filter!(x -> x <= 25, mpirange)

if localstart_before_filter > 25
    @test forlocalpart(lp -> length(lp) == 0, mpirange)
end

sync(mpirange)

redist_range = redistribute(mpirange)
@test forlocalpart(lp -> length(lp) >= 25 ÷ nb_procs, redist_range)

end

MPI.Finalize()


text:
include("mybicgstab.jl")

#one-sided FEAST for linear eigenvalue problems
function feast_linear(A,B,x0,nc,emid,ra,rb,eps,maxit;log=false, insideEps=1e-2, set_nin=false)
	 (n,m0)=size(x0)

	 integrand(z,x,lest,data,resvecs)=\(z*B-A,B*x)

	 function rrsolve(Q)
		  Aq=Q'*A*Q
		  Bq=Q'*B*Q

		  le,xe=eigen(Aq,Bq)
		  return (le,Q*xe)
	 end

	 Tf(l,x)=(B*x*diagm( 0 => l)-A*x)

	 return feast_core(Tf,integrand,rrsolve,x0,nc,emid,ra,rb,eps,maxit;log=log,insideEps=insideEps, set_nin=set_nin)

end

using SparseArrays
#two-sided FEAST for linear eigenvalue problems
function feastNS_linear(A,B,x0,y0,nc,emid,ra,rb,eps,maxit; log=false,insideEps=1e-2, solver=(\), factorizer=lu, FactoType=SuiteSparse.UMFPACK.UmfpackLU{promote_type(A,B),Int64})
	 (n,m0)=size(x0)
	 if eltype(A) <: Complex || eltype(B) <: Complex
			 #A32, B32 = convert(SparseMatrixCSC{ComplexF32,Int64}, A), convert(SparseMatrixCSC{ComplexF32,Int64}, B)
	else
		#A32, B32 = convert(SparseMatrixCSC{Float32,Int64}, A), convert(SparseMatrixCSC{Float32,Int64}, B)
	end
	A32, B32 = A, B

	integrand, hintegrand, facto, hfacto, Tf, hTf = let A32=A32, B32=B32, A=A, B=B
		 ((factor,x,lest,data,resvecs)->solver(factor,B32*x),
		 (factor,y,lest,hdata,hresvecs)->solver(factor,B32'*y),
		 z->factorizer(z*B32-A32),
		 z->factorizer(convert(SparseMatrixCSC, (z*B32-A32)')) ,
		 (l,x)->begin
			 println("X times EV")
			 tmp1 = x*spdiagm(0 => l)
			 println("B times (X times EV)")
			 B1 = B*tmp1
			 println("A times X")
			 A1 = A*x
			 return B1 - A1
		end,
		(l,y)->begin
			println("Y times EV*")
			tmp1 = y*spdiagm(0 => conj.(l))
			println("B* times (Y times EV*)")
			B1 = B'*tmp1
			println("A* times X")
			A1 = A'*y
			return B1 - A1
	  end)
	 end
	 function rrsolvens(Q,R)
		  Aq=R'*(A*Q)
		  Bq=R'*(B*Q)
		  (ye,le,xe)=nseig(Aq,Bq)
		  return (R*ye,le,Q*xe)
	 end

	 #hTf(l,y)=B'*(y*spdiagm(0 => l)')-A'*y

	 #use SVD to B-biorthogonalize subspaces
	 function biortho(Q,R)
		  Bq=R'*(B*Q)
		  (u,s,v)=svd(Bq)
		  Y=R*(u*diagm( 0 => 1 ./sqrt.(s)))
		  X=Q*(v*diagm( 0 => 1 ./sqrt.(s)))
		  return (X,Y)
	 end
	 return feastNS_core(Tf,hTf,integrand,hintegrand,facto, hfacto,FactoType,rrsolvens,biortho,x0,y0,nc,emid,ra,rb,eps,maxit; log=log, insideEps=insideEps)

end

#generalized inexact FEAST algorithm with gmres
function ifeast_linear(A,B,x0,alpha,isMaxit,nc,emid,ra,rb,eps,maxit;log=false,insideEps=1e-2, verbose=1)
	 (n,m0)=size(x0)

	 function integrand(z,x,lest,data,resvecs)
		  nc=data[:shiftIndex][z]
		  M=(z*B-A)
		  #b=(B*x*diagm(lest)-A*x)
		  #normb=maximum(sqrt.(sum(abs.(b).^2,1)))

		  int=zeros(ComplexF64,n,m0)
		  maxits=0
		  rhs=convert(Array{ComplexF64,2},resvecs)

		  for i in 1:m0
				#(int[:,i],history)=bicgstabl(M,resvecs[:,i],1,max_mv_products=isMaxit,tol=alpha,initial_zero=true,log=true)
				#(int[:,i],history)=idrs(M,rhs[:,i];maxiter=isMaxit,tol=alpha,log=true)
				(int[:,i],history)=gmres(M,resvecs[:,i],restart=isMaxit,tol=alpha,initially_zero=true,maxiter=isMaxit,log=true)
				#(int[:,i],history)=minres(M,resvecs[:,i],tol=alpha,initially_zero=true,maxiter=isMaxit,log=true)
				nlinits=size(history[:resnorm],1)
				data[:linIts][i,nc,data[:iterations]]=nlinits
				#data[:linResiduals][i,nc,data[:iterations]]=history[:resnorm][nlinits]
				data[:linResiduals][i,nc,data[:iterations]]=norm(resvecs[:,i]-M*int[:,i])/norm(resvecs[:,i])

				if(nlinits>maxits)
					 maxits=nlinits
				end
		  end
		  #int=\(M,resvecs)
		  #println("      linits=$maxits")
		  #int=zbicgstabBlock(M,resvecs,zeros(n,m0),isMaxit,alpha)
		  return (x-int)*spdiagm(0 => 1 ./(z.-lest))
	 end

	 function rrsolve(Q)
		  Aq=Q'*A*Q
		  Bq=Q'*B*Q
		  le,xe=eigen(Aq,Bq)
		  return (le,Q*xe)
	 end

	 Tf(l,x)=(B*x*spdiagm(0 => l)-A*x)

	 return feast_core(Tf,integrand,rrsolve,x0,nc,emid,ra,rb,eps,maxit;log=log,insideEps=insideEps, verbose=verbose)

end

#generalized ifeast with bicgstab
function ifeast_linearBicgstab(A,B,x0,alpha,isMaxit,nc,emid,ra,rb,eps,maxit;log=false,insideEps=1e-2,verbose=1)
	 (n,m0)=size(x0)

	 function integrand(z,x,lest,data,resvecs)
		  nc=data[:shiftIndex][z]
		  M1=(z*B-A)
		  P=speye(n)
		  M=P*M1*P
		  rhs=P*resvecs
		  #b=(B*x*diagm(lest)-A*x)
		  #normb=maximum(sqrt.(sum(abs.(b).^2,1)))

		  int=zeros(ComplexF64,n,m0)

		  for i in 1:m0
				(int[:,i],history)=bicgstabl(M,resvecs[:,i],1,max_mv_products=isMaxit,tol=alpha,initial_zero=true,log=true)
				nlinits=size(history[:resnorm],1)
				data[:linIts][i,nc,data[:iterations]]=nlinits
				data[:linResiduals][i,nc,data[:iterations]]=history[:resnorm][nlinits]
		  end

		  return (x-int)*spdiagm(0 => 1 ./(z.-lest))
	 end

	 function rrsolve(Q)
		  Aq=Q'*A*Q
		  Bq=Q'*B*Q
		  le,xe=eigen(Aq,Bq)
		  return (le,Q*xe)
	 end

	 Tf(l,x)=(B*x*spdiagm(0 => l)-A*x)

	 return feast_core(Tf,integrand,rrsolve,x0,nc,emid,ra,rb,eps,maxit;log=log,insideEps=insideEps,verbose=verbose)

end

#inexact two-sided ifeast with my own bicgstab implementation
function ifeastNS_linear(A,B,x0,y0,alpha,isMaxit,nc,emid,ra,rb,eps,maxit; log=false,insideEps=1e-2)
	 (n,m0)=size(x0)

	 #integrand(z,x,lest,data,resvecs)=\(z*B-A,B*x)
	 #hintegrand(z,y,lest,hdata,hresvecs)=\((z*B-A)',B'*y)
	 function integrand(z,x,lest,data,resvecs)
		  M=z*B-A
		  rhs=B*x*spdiagm(0 => lest)-A*x
		  int0=zeros(ComplexF64,n,m0)
		  int=zeros(ComplexF64,n,m0)
		  int=zbicgstabBlock(M,resvecs,int0,isMaxit,alpha)
		  for i in 1:m0
				#(int[:,i],history)=gmres(M,rhs[:,i],restart=isMaxit,tol=alpha,initially_zero=true,maxiter=isMaxit,log=true)
				#(int[:,i],history)=bicgstabl(M,rhs[:,i],2,max_mv_products=isMaxit,tol=alpha,initial_zero=true,log=true)
		  end

		  #int=\(M,resvecs)
		  return (x-int)*spdiagm(0 => 1 ./(z.-(lest)))
	 end

	 function hintegrand(z,y,lest,hdata,hresvecs)
		  M=(z*B-A)'
		  rhs=B'*y*spdiagm(0 => lest)'-A'*y
		  int0=zeros(ComplexF64,n,m0)
		  int=zeros(ComplexF64,n,m0)
		  int=zbicgstabBlock(M,hresvecs,int0,isMaxit,alpha)
		  for i in 1:m0
				#(int[:,i],history)=gmres(M,rhs[:,i],restart=isMaxit,tol=alpha,initially_zero=true,maxiter=isMaxit,log=true)

				#(int[:,i],history)=bicgstabl(M,rhs[:,i],2,max_mv_products=isMaxit,tol=alpha,initial_zero=true,log=true)
		  end
		  #int=\(M,hresvecs)
		  return (y-int)*spdiagm(0 => 1 ./(z'.-conj.(lest)))
	 end

	 function rrsolvens(Q,R)
		  Aq=R'*A*Q
		  Bq=R'*B*Q
		  (ye,le,xe)=nseig(Aq,Bq)
		  return (R*ye,le,Q*xe)
	 end

	 Tf(l,x)=(B*x*spdiagm(0 => l)-A*x)
	 hTf(l,y)=B'*y*spdiagm(0 => l)'-A'*y

	 #use SVD to B-biorthogonalize subspaces
	 function biortho(Q,R)
		  Bq=R'*B*Q
		  (u,s,v)=svd(Bq)
		  Y=R*u*diagm( 0 => 1 ./sqrt.(s))
		  X=Q*v*diagm( 0 => 1 ./sqrt.(s))
		  return (X,Y)
	 end

	 return feastNS_core(Tf,hTf,integrand,hintegrand,rrsolvens,biortho,x0,y0,nc,emid,ra,rb,eps,maxit; log=log, insideEps=insideEps)

end

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
const ConcreteFloat = Union{Float64,Float32,Float16,BigFloat}
const ConcreteInt = Union{Bool,Int8,Int16,Int32,Int64,Int128,BigInt}

abstract type AbstractParticles{T,N} <: Real end
"""
    struct Particles{T, N} <: AbstractParticles{T, N}

This type represents uncertainty using a cloud of particles.
# Constructors:
- `Particles()`
- `Particles(N::Integer)`
- `Particles([rng::AbstractRNG,] d::Distribution)`
- `Particles([rng::AbstractRNG,] N::Integer, d::Distribution; permute=true, systematic=true)`
- `Particles(v::Vector{T} where T)`
- `Particles(m::Matrix{T} where T)`: Creates multivariate particles (Vector{Particles})
"""
struct Particles{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
end

"""
    struct StaticParticles{T, N} <: AbstractParticles{T, N}

See `?Particles` for help. The difference between `StaticParticles` and `Particles` is that the `StaticParticles` store particles in a static vecetor. This makes runtimes much shorter, but compile times longer. See the documentation for some benchmarks. Only recommended for sample sizes of ≲ 300-400
"""
struct StaticParticles{T,N} <: AbstractParticles{T,N}
    particles::SArray{Tuple{N}, T, 1, N}
end

DNP(PT) = PT === Particles ? DEFAULT_NUM_PARTICLES : DEFAULT_STATIC_NUM_PARTICLES

const ParticleSymbols = (:Particles, :StaticParticles)

for PT in ParticleSymbols
    for D in (2,3,4,5)
        @eval function $PT{T,N}(m::AbstractArray{T,$D}) where {T,N}
            size(m, 1) == N || throw(ArgumentError("The first dimension of the array must be the same as the number N of particles."))
            inds = CartesianIndices(axes(m)[2:end])
            map(inds) do ind
                $PT{T,N}(@view(m[:,ind]))
            end
        end

        @eval function $PT(m::AbstractArray{T,$D}) where T
            N = size(m, 1)
            inds = CartesianIndices(axes(m)[2:end])
            map(inds) do ind
                $PT{T,N}(@view(m[:,ind]))
            end
        end
    end

    @eval begin
        $PT(v::Vector) = $PT{eltype(v),length(v)}(v)

        function $PT{T,N}(n::Real) where {T,N} # This constructor is potentially dangerous, replace with convert?
            if n isa AbstractParticles
                return convert($PT{T,N}, n)
            end
            v = fill(n,N)
            $PT{T,N}(v)
        end

        $PT{T,N}(p::$PT{T,N}) where {T,N} = p

        function $PT(rng::AbstractRNG, N::Integer=DNP($PT), d::Distribution{<:Any,VS}=Normal(0,1); permute=true, systematic=VS==Continuous) where VS
            if systematic
                v = systematic_sample(rng,N,d; permute=permute)
            else
                v = rand(rng, d, N)
            end
            $PT{eltype(v),N}(v)
        end
        function $PT(N::Integer=DNP($PT), d::Distribution{<:Any,VS}=Normal(0,1); kwargs...) where VS
            return $PT(Random.GLOBAL_RNG, N, d; kwargs...)
        end

        function $PT(::Type{T}, N::Integer=DNP($PT), d::Distribution=Normal(T(0),T(1)); kwargs...) where {T <: Real}
            eltype(d) == T || throw(ArgumentError("Element type of the provided distribution $d does not match $T. The element type of a distribution is the element type of a value sampled from it. Some distributions, like `Gamma(0.1f0)` generated `Float64` random number even though it appears like it should generate `Float32` numbers."))
            return $PT(Random.GLOBAL_RNG, N, d; kwargs...)
        end

        function $PT(rng::AbstractRNG, N::Integer, d::MultivariateDistribution)
            v = rand(rng,d,N)' |> copy # For cache locality
            $PT(v)
        end
        $PT(N::Integer, d::MultivariateDistribution) = $PT(Random.GLOBAL_RNG, N, d)

        nakedtypeof(p::$PT{T,N}) where {T,N} = $PT
        nakedtypeof(::Type{$PT{T,N}}) where {T,N} = $PT
    end
end

function Particles(rng::AbstractRNG, d::Distribution;kwargs...)
    Particles(rng, DEFAULT_NUM_PARTICLES, d; kwargs...)
end
Particles(d::Distribution; kwargs...) = Particles(Random.GLOBAL_RNG, d; kwargs...)

function StaticParticles(rng::AbstractRNG, d::Distribution;kwargs...)
    StaticParticles(rng, DEFAULT_STATIC_NUM_PARTICLES, d; kwargs...)
end
StaticParticles(d::Distribution;kwargs...) = StaticParticles(Random.GLOBAL_RNG, d; kwargs...)

const MvParticles = Vector{<:AbstractParticles} # This can not be AbstractVector since it causes some methods below to be less specific than desired
const ParticleArray = AbstractArray{<:AbstractParticles}
const SomeKindOfParticles = Union{<:AbstractParticles, ParticleArray}

Particles(p::StaticParticles{T,N}) where {T,N} = Particles{T,N}(p.particles)


text:
function pseudo_marginal_full(prng, data, choices, scale, σ², ϵ², ℓ², u)
#=
    Filippone, Maurizio, and Mark Girolami.
    "Pseudo-marginal Bayesian inference for Gaussian processes."
    IEEE Transactions on Pattern Analysis and Machine Intelligence (2014)
=##
    try
        K          = compute_gram_matrix(data, σ², ϵ², ℓ²)
        K          = PDMats.PDMat(K)
        μ, Σ, a, B = laplace_approximation(K, choices, zeros(length(u)), scale;
                                           verbose=false)
        Σ          = PDMats.PDMat(Σ)
        q          = MvNormal(μ, Σ)
        f_sample   = PDMats.unwhiten(q.Σ, u) + q.μ
        joint      = logjointlike(K, choices, f_sample, scale)
        logpml     = joint - logpdf(q, f_sample)
        return logpml, f_sample, q, K, Σ, a
    catch err
        # if(isa(err, LinearAlgebra.PosDefException))
        #     @warn "Cholesky failed. Rejecting proposal"
        #     return -Inf, nothing, nothing, nothing, nothing, nothing
        # else
        #     throw(err)
        # end
        @warn "Cholesky or determinant failed. Rejecting proposal"
        return -Inf, nothing, nothing, nothing, nothing, nothing
    end
end

@inline function pseudo_marginal_partial(prng, choices, q, K, scale, u)
#=
    Filippone, Maurizio, and Mark Girolami.
    "Pseudo-marginal Bayesian inference for Gaussian processes."
    IEEE Transactions on Pattern Analysis and Machine Intelligence (2014)
=##
    f_sample = PDMats.unwhiten(q.Σ, u) + q.μ
    joint    = logjointlike(K, choices, f_sample, scale)
    logpml   = joint - logpdf(q, f_sample)
    logpml
end

function ess_transition(prng, loglike, prev_θ, prev_like, prior)
#=
    Murray, Iain, Ryan Adams, and David MacKay.
    "Elliptical slice sampling."
    Artificial Intelligence and Statistics. 2010.
=##
    @label propose
    chol_fail = 0
    ν     = rand(prng, prior)
    μ     = mean(prior)
    u     = rand(prng)
    logy  = prev_like + log(u)
    ϵ     = rand(prng)*2*π
    ϵ_min = ϵ - 2*π
    ϵ_max = deepcopy(ϵ)

    proposals = 1
    while(true)
        cosϵ   = cos(ϵ)
        sinϵ   = sin(ϵ)
        a      = 1 - (cosϵ + sinϵ)

        prop_θ    = @. cosϵ*prev_θ + sinϵ*ν + a*μ
        prop_logp = loglike(prop_θ)

        if(isinf(prop_logp))
             chol_fail += 1
             if(chol_fail > 10)
                 throw(LinearAlgebra.PosDefException(1))
             end
             @goto propose
        end

        if(prop_logp > logy)
            return prop_θ, prop_logp, proposals
        else
            if(ϵ < 0)
                ϵ_min = deepcopy(ϵ)
            else
                ϵ_max = deepcopy(ϵ)
            end
            ϵ = rand(prng, Uniform(ϵ_min, ϵ_max))

            proposals += 1
        end
    end
end

function pm_ess(prng,
                samples::Int,
                warmup::Int,
                thinning::Int,
                initial_θ::Vector{<:Real},
                initial_f::Vector{<:Real},
                prior,
                scale::Real,
                data::Matrix{<:Real},
                choices::Matrix{<:Int})
#=
    Murray, Iain, and Matthew Graham.
    Pseudo-marginal slice sampling.
    Artificial Intelligence and Statistics. 2016.
=##
    θ_samples = zeros(length(initial_θ), samples)
    f_samples = zeros(length(initial_f), samples)
    a_samples = zeros(length(initial_f), samples)
    K_samples = Array{PDMats.PDMat}(undef, samples)
    prev_logpml = -Inf

    u_prior = MvNormal(length(initial_f), 1.0)

    θ_map, f, Σ, a, B, K = map_laplace(data, choices, initial_θ, scale, prior;
                                       verbose=true)
    θ      = θ_map
    u      = rand(prng, u_prior)
    q      = MvNormal(f, Σ)
    logpml = pseudo_marginal_partial(prng, choices, q, K, scale, u)

    u_acc_mavg = OnlineStats.Mean()
    θ_acc_mavg = OnlineStats.Mean()
    prog       = ProgressMeter.Progress(samples+warmup)
    for i = 1:(samples+warmup)
        Lu = u_in->begin
            pseudo_marginal_partial(prng, choices, q, K, scale, u_in)
        end
        u, logpml, u_nprop = ess_transition(prng, Lu, u, logpml, u_prior)

        Lθ = θ_in->begin
            logpml, f, q, K, Σ, a = pseudo_marginal_full(
                prng, data, choices, scale,
                θ_in[1], θ_in[2], θ_in[3:end], u)
            logpml
        end
        θ, logpml, θ_nprop = ess_transition(prng, Lθ, θ, logpml, prior)

        if(i > warmup)
            θ_samples[:,i-warmup] = θ
            f_samples[:,i-warmup] = f
            a_samples[:,i-warmup] = a
            K_samples[i-warmup]   = K
        end
        u_acc = 1/u_nprop
        θ_acc = 1/θ_nprop
        OnlineStats.fit!(u_acc_mavg, u_acc)
        OnlineStats.fit!(θ_acc_mavg, θ_acc)
        ProgressMeter.next!(
            prog; showvalues=[(:iteration, i),
                              (:pseudo_marginal_loglikelihood, logpml),
                              (:u_acceptance, u_acc),
                              (:θ_acceptance, θ_acc),
                              (:u_average_acceptance, u_acc_mavg.μ),
                              (:θ_average_acceptance, θ_acc_mavg.μ)
                              ])
    end
    θ_samples = θ_samples[:,1:thinning:end]
    f_samples = f_samples[:,1:thinning:end]
    a_samples = a_samples[:,1:thinning:end]
    K_samples = K_samples[1:thinning:end]

    θ_samples, f_samples, a_samples, K_samples
end


text:
"""
Performs crossover step for each particle pt in the chain
* `model`: model containing a likelihood function with data and priors
* `de`: differential evolution object
* `group`: a group of particles
"""
function crossover!(model, de, group)
    for pt in group
        # generate the proposal
        proposal = generate_proposal(de, pt, group)
        # compute the weight of the proposal: prior loglikelihood + data loglikelihood
        proposal.weight = priorlike(model, proposal) + model.model(proposal.Θ)
        # accept proposal according to Metropolis-Hastings rule
        update_particle!(de, pt, proposal)
    end
end

"""
Generate proposal according to θ' = θt + γ1(θm − θn) + γ2(θb − θt) + b
γ2=0 after burnin
* `de`: differential evolution object
* `Pt`: current particle
* `group`: a group of particles
"""
# 2.38/sqrt(2d)
function generate_proposal(de, Pt, group)
    Np = length(Pt.Θ)
    # select the base particle θb
    Pb = select_base(group)
    # group without Pb
    group_diff = setdiff(group, [Pt])
    # sample particles for θm and θn
    Pm,Pn = sample(group_diff, 2, replace=false)
    # sample gamma weights
    γ₁ = rand(Uniform(.5, 1))
    # set γ₂ = 0 after burnin
    de.iter > de.burnin ? γ₂=0.0 : γ₂=rand(Uniform(.5, 1))
    # sample noise for each parameter
    b = Uniform(-de.ϵ, de.ϵ)
    # compute proposal value
    Θp = Pt + γ₁*(Pm - Pn) + γ₂*(Pb - Pt) + b
    # enforce the parameter boundaries
    enforce_bounds!(de.bounds, Θp)
    # reset each parameter to previous value with probability (1-κ)
    recombination!(de, Pt, Θp)
    return Θp
end

"""
Selects base particle θb with probability proportional to weight
* `group`: a group of particles
"""
function select_base(group)
    w = map(x->x.weight, group)
    θ = exp.(w)/sum(exp.(w))
    p = sample(group, Weights(θ))
    return p
end

"""
Resets parameters of proposal to previous value with probability
(1-κ).
* `pp`: proposal particle
* `pt`: current partical
"""
function recombination!(de, pt, pp)
    de.κ == 1.0 ? (return nothing) : nothing
    N = length(pt.Θ)
    for i in 1:N
        if isa(pp.Θ[i], Array)
            recombination!(de, pt, pp, i)
        else
            rand() <= (1-de.κ) ? pp.Θ[i] = pt.Θ[i] : nothing
        end
    end
    return nothing
end

# Handles elements within an array
function recombination!(de, pt, pp, idx)
    N = length(pt.Θ[idx])
    for i in 1:N
        rand() <= (1-de.κ) ? pp.Θ[idx][i] = pt.Θ[idx][i] : nothing
    end
    return nothing
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RenewableDispatch <: RenewableGen
        name::String
        available::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        rating::Float64
        primemover::PrimeMovers.PrimeMover
        reactivepowerlimits::Union{Nothing, Min_Max}
        powerfactor::Float64
        op_cost::TwoPartCost
        basepower::Float64
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `activepower::Float64`
- `reactivepower::Float64`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: (0, nothing), action if invalid: error
- `primemover::PrimeMovers.PrimeMover`: PrimeMover Technology according to EIA 923
- `reactivepowerlimits::Union{Nothing, Min_Max}`
- `powerfactor::Float64`, validation range: (0, 1), action if invalid: error
- `op_cost::TwoPartCost`: Operation Cost of Generation [`TwoPartCost`](@ref)
- `basepower::Float64`: Base power of the unit in system base per unit, validation range: (0, nothing), action if invalid: warn
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct RenewableDispatch <: RenewableGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "PrimeMover Technology according to EIA 923"
    primemover::PrimeMovers.PrimeMover
    reactivepowerlimits::Union{Nothing, Min_Max}
    powerfactor::Float64
    "Operation Cost of Generation [`TwoPartCost`](@ref)"
    op_cost::TwoPartCost
    "Base power of the unit in system base per unit"
    basepower::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function RenewableDispatch(name, available, bus, activepower, reactivepower, rating, primemover, reactivepowerlimits, powerfactor, op_cost, basepower, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    RenewableDispatch(name, available, bus, activepower, reactivepower, rating, primemover, reactivepowerlimits, powerfactor, op_cost, basepower, services, dynamic_injector, ext, forecasts, InfrastructureSystemsInternal(), )
end

function RenewableDispatch(; name, available, bus, activepower, reactivepower, rating, primemover, reactivepowerlimits, powerfactor, op_cost, basepower, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    RenewableDispatch(name, available, bus, activepower, reactivepower, rating, primemover, reactivepowerlimits, powerfactor, op_cost, basepower, services, dynamic_injector, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function RenewableDispatch(::Nothing)
    RenewableDispatch(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        rating=0.0,
        primemover=PrimeMovers.OT,
        reactivepowerlimits=nothing,
        powerfactor=1.0,
        op_cost=TwoPartCost(nothing),
        basepower=0.0,
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::RenewableDispatch) = value.name
"""Get RenewableDispatch available."""
get_available(value::RenewableDispatch) = value.available
"""Get RenewableDispatch bus."""
get_bus(value::RenewableDispatch) = value.bus
"""Get RenewableDispatch activepower."""
get_activepower(value::RenewableDispatch) = value.activepower
"""Get RenewableDispatch reactivepower."""
get_reactivepower(value::RenewableDispatch) = value.reactivepower
"""Get RenewableDispatch rating."""
get_rating(value::RenewableDispatch) = value.rating
"""Get RenewableDispatch primemover."""
get_primemover(value::RenewableDispatch) = value.primemover
"""Get RenewableDispatch reactivepowerlimits."""
get_reactivepowerlimits(value::RenewableDispatch) = value.reactivepowerlimits
"""Get RenewableDispatch powerfactor."""
get_powerfactor(value::RenewableDispatch) = value.powerfactor
"""Get RenewableDispatch op_cost."""
get_op_cost(value::RenewableDispatch) = value.op_cost
"""Get RenewableDispatch basepower."""
get_basepower(value::RenewableDispatch) = value.basepower
"""Get RenewableDispatch services."""
get_services(value::RenewableDispatch) = value.services
"""Get RenewableDispatch dynamic_injector."""
get_dynamic_injector(value::RenewableDispatch) = value.dynamic_injector
"""Get RenewableDispatch ext."""
get_ext(value::RenewableDispatch) = value.ext

InfrastructureSystems.get_forecasts(value::RenewableDispatch) = value.forecasts
"""Get RenewableDispatch internal."""
get_internal(value::RenewableDispatch) = value.internal

InfrastructureSystems.set_name!(value::RenewableDispatch, val::String) = value.name = val
"""Set RenewableDispatch available."""
set_available!(value::RenewableDispatch, val::Bool) = value.available = val
"""Set RenewableDispatch bus."""
set_bus!(value::RenewableDispatch, val::Bus) = value.bus = val
"""Set RenewableDispatch activepower."""
set_activepower!(value::RenewableDispatch, val::Float64) = value.activepower = val
"""Set RenewableDispatch reactivepower."""
set_reactivepower!(value::RenewableDispatch, val::Float64) = value.reactivepower = val
"""Set RenewableDispatch rating."""
set_rating!(value::RenewableDispatch, val::Float64) = value.rating = val
"""Set RenewableDispatch primemover."""
set_primemover!(value::RenewableDispatch, val::PrimeMovers.PrimeMover) = value.primemover = val
"""Set RenewableDispatch reactivepowerlimits."""
set_reactivepowerlimits!(value::RenewableDispatch, val::Union{Nothing, Min_Max}) = value.reactivepowerlimits = val
"""Set RenewableDispatch powerfactor."""
set_powerfactor!(value::RenewableDispatch, val::Float64) = value.powerfactor = val
"""Set RenewableDispatch op_cost."""
set_op_cost!(value::RenewableDispatch, val::TwoPartCost) = value.op_cost = val
"""Set RenewableDispatch basepower."""
set_basepower!(value::RenewableDispatch, val::Float64) = value.basepower = val
"""Set RenewableDispatch services."""
set_services!(value::RenewableDispatch, val::Vector{Service}) = value.services = val
"""Set RenewableDispatch ext."""
set_ext!(value::RenewableDispatch, val::Dict{String, Any}) = value.ext = val

InfrastructureSystems.set_forecasts!(value::RenewableDispatch, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set RenewableDispatch internal."""
set_internal!(value::RenewableDispatch, val::InfrastructureSystemsInternal) = value.internal = val


text:
using CImGui
using CImGui.CSyntax
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using CImGuiExtensions
using CImGui: ImVec2, ImVec4, IM_COL32, ImU32
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using DataFrames # For this example, don't forget to add this package.
using Printf

function create_labelled_interval_context(data::AbstractArray)
    padding = (8,8,0,0)
    # The plot context draws an overview of the data.
    model = PlotlinesModel(data)
    outline_layout = RectangularLayout(pos = ImVec2(0,0), width = Cfloat(800), height = Cfloat(80))
    display_properties = PlotlinesDisplayProperties(id = "###overview"; caption = x -> "", createwindow = false, layout = outline_layout, xaxis = Axis(false, Cfloat(1)), yaxis = Axis(false, Cfloat(1)), padding = padding)
    plot_outline_context = PlotContext(PlotlinesControl(true), model, display_properties)

    # The nested interval keeps track of which portion of the data the user wants to select.
    nested_interval = NestedInterval(start = Float64(1.0), stop = Float64(length(data)), interval = 1:length(data))

    # The label interval context allows the user to label and mark a selected interval.
    label_control = LabelIntervalControl(true, nested_interval, "\0"^64*"\0")
    label_display_properties = LabelIntervalDisplayProperties(id = "###event markers" ; caption = "Conditions", createwindow =  false, col = ImVec4(0.0, 0, 0.9, 0.1), textcol = ImVec4(1.0, 1.0, 0.2, 1), layout = outline_layout, padding = padding)
    dict = Dict{String, LabelledInterval}()
    labelled_intervals = LabelledIntervals("Conditions", dict)
    label_interval_context = LabelIntervalContext(label_control, labelled_intervals, label_display_properties, plot_outline_context)
    return label_interval_context
end

function label_interval_example()

    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1280, 720, "Demo")
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    CImGui.StyleColorsDark()

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)

    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    data = [rand()*0.2 for i = 1:100]
    label_interval_context = create_labelled_interval_context(data)

    while !GLFW.WindowShouldClose(window)

        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        CImGui.Begin("Demo: Label Interval")
        @cstatic i1=Cint(1) i2=Cint(100) begin

        if isrunning(label_interval_context)
            !isnothing(label_interval_context) ? label_interval_context() : nothing
            control = label_interval_context.control
            i1 = Cint(get_start(control))
            i2 = Cint(get_stop(control))
        end
        CImGui.PushItemWidth(50)
        CImGui.SameLine()
        @c CImGui.DragInt("drag start", &i1, 1, 1, 100, "%d")
        CImGui.SameLine()
        @c CImGui.DragInt("drag stop", &i2, 1, 1, 100, "%d")
        i1 = i1 < i2 ? i1 : i2
        i2 = i2 > i1 ? i2 : i1
        if !isnothing(label_interval_context)
            control = label_interval_context.control
            set_start!(control, i1)
            set_stop!(control, i2)
        end
        CImGui.PopItemWidth()
        end

        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

    # cleanup
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end

# Run the demo
 label_interval_example()


text:
module LambdaCalculus

include("common_functions.jl")
include("naming_context.jl")

include("Named/Named.jl")
include("DeBruijn/DeBruijn.jl")
include("LocallyNameless/LocallyNameless.jl")

include("conversions.jl")

end # module LambdaCalculus


text:
# Elementary Givens rotation
using LinearAlgebra: Givens

# for reference: Julia LinearAlgebra.Givens
# immutable Givens{T} <: AbstractRotation{T}
#     i1::Int
#     i2::Int
#     c::T
#     s::T
# end

function LinearAlgebra.rmul!(b::OrthonormalBasis{T}, G::Givens) where {T}
    if T isa AbstractArray && IndexStyle(T) isa IndexLinear
        return _rmul_linear!(b, G)
    else
        return _rmul!(b, G)
    end
end

@fastmath function _rmul_linear!(b::OrthonormalBasis{<:AbstractArray}, G::Givens)
    q1, q2 = b[G.i1], b[G.i2]
    c = G.c
    s = G.s
    @inbounds @simd for i in 1:length(q1)
        q1[i], q2[i] = c * q1[i] - conj(s) * q2[i], s * q1[i] + c * q2[i]
    end
    return b
end

function _rmul!(b::OrthonormalBasis, G::Givens)
    q1, q2 = b[G.i1], b[G.i2]
    q1old = mul!(similar(q1), q1, true)
    q1 = axpby!(-conj(G.s), q2, G.c, q1)
    q2 = axpby!(G.s, q1old, G.c, q2)
    return b
end

# New types for discarding or for storing successive Givens transformations
struct NoVecs end
const novecs = NoVecs()
rmulc!(::NoVecs, ::Any) = novecs


text:
# Copied/adapted from code under:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    AbstractOpticalSystem{T<:Real}

Abstract type for any optical system, must parameterized by the datatype of
entities within the system `T` (numeric).

The OpticalSystem will hold relevant global parameters, as well as holding
all of the objects in the scene.
"""
abstract type AbstractOpticalSystem{T <: Real} end

# ============================================================================ #
# ====                        CSGOpticalSystem                            ==== #
# ============================================================================ #

"""
    CSGOpticalSystem{T,D<:Real,S<:Surface{T},L<:LensAssembly{T}} <: AbstractOpticalSystem{T}

An optical system containing a lens assembly with all optical elements and a
detector surface with associated image. The system can be at a specified
temperature and pressure.

There are two number types in the type signature.
The `T` type parameter is the numeric type for geometry in the optical system,
the `D` type parameter is the numeric type of the pixels in the detector image.
This way you can have `Float64` geometry, where high precision is essential,
but the pixels in the detector can be `Float32` since precision is much less
critical for image data, or Complex if doing wave optic simulations.

The detector can be any [`Surface`](@ref)
which implements [`uv`](@ref), [`uvtopix`](@ref) and [`onsurface`](@ref),
typically this is one of
[`Rectangle`](@ref), [`Ellipse`](@ref) or [`SphericalCap`](@ref).

```julia
CSGOpticalSystem(
    assembly::LensAssembly,
    detector::Surface,
    detectorpixelsx = 1000,
    detectorpixelsy = 1000, ::Type{D} = Float32;
    temperature = OpticSim.GlassCat.TEMP_REF,
    pressure = OpticSim.GlassCat.PRESSURE_REF
)
```
"""
struct CSGOpticalSystem{
    T,
    D <: Number,
    S <: Surface{T}, # should allow for more than one detector?
    L <: LensAssembly{T}, # need a new type for 'objects_in_scene' type thing
} <: AbstractOpticalSystem{T}
    assembly::L
    detector::S
    detectorimage::HierarchicalImage{D}
    temperature::T
    pressure::T

    function CSGOpticalSystem(
        assembly::L,
        detector::S,
        detectorpixelsx::Int = 1000,
        detectorpixelsy::Int = 1000,
        ::Type{D} = Float32;
        temperature::Union{T, Unitful.Temperature} = convert(T, TEMP_REF),
        pressure::T = convert(T, PRESSURE_REF),
    ) where {T <: Real, S <: Surface{T}, L <: LensAssembly{T}, D <: Number}
        @assert hasmethod(uv, (S, SVector{3, T})) "Detector must implement uv()"
        @assert hasmethod(uvtopix, (S, SVector{2, T}, Tuple{Int, Int})) "Detector must implement uvtopix()"
        @assert hasmethod(onsurface, (S, SVector{3, T})) "Detector must implement onsurface()"
        opticalinterface = interface(detector)
        @assert insidematerialid(opticalinterface) ==
                outsidematerialid(opticalinterface) "Detector must have same material either side"
        @assert interface(detector) !== NullInterface(T) "Detector can't have null interface"
        image = HierarchicalImage{D}(detectorpixelsy, detectorpixelsx)
        if temperature isa Unitful.Temperature
            temperature = Unitful.ustrip(T, °C, temperature)
        end
        return new{T, D, S, L}(
            assembly,
            detector,
            image,
            temperature,
            convert(T, pressure),
        )
    end
end

function Base.copy(a::CSGOpticalSystem)
    return CSGOpticalSystem(
        a.assembly,
        a.detector,
        size(a.detectorimage)...,
        temperature = (a.temperature)°C,
        pressure = a.pressure,
    )
end

# added this show method because the type of CSGOpticalSystem is gigantic and printing it in the REPL can crash the
# system
function Base.show(io::IO, a::CSGOpticalSystem{T}) where {T}
    return print(
        io,
        "CSGOpticalSystem{$T}($(temperature(a))," *
        " $(pressure(a)), $(assembly(a)), $(detector(a)))",
    )
end

"""
    assembly(system::AbstractOpticalSystem{T}) -> LensAssembly{T}

Get the [`LensAssembly`](@ref) of `system`.
"""
assembly(system::CSGOpticalSystem{T}) where {T <: Real} = system.assembly

detector(system::CSGOpticalSystem) = system.detector

"""
    detectorimage(system::AbstractOpticalSystem{T}) -> HierarchicalImage{D}

Get the detector image of `system`.
`D` is the datatype of the detector image and is not necessarily the same as
the datatype of the system `T`.
"""
detectorimage(system::CSGOpticalSystem) = system.detectorimage

detectorsize(system::CSGOpticalSystem) = size(system.detectorimage)

"""
    temperature(system::AbstractOpticalSystem{T}) -> T

Get the temperature of `system` in °C.
"""
temperature(system::CSGOpticalSystem{T}) where {T <: Real} = system.temperature

"""
    pressure(system::AbstractOpticalSystem{T}) -> T

Get the pressure of `system` in Atm.
"""
pressure(system::CSGOpticalSystem{T}) where {T <: Real} = system.pressure
"""
    resetdetector!(system::AbstractOpticalSystem{T})

Reset the deterctor image of `system` to zero.
"""
function resetdetector!(system::CSGOpticalSystem{T}) where {T <: Real}
    return reset!(system.detectorimage)
end

# not sure what the below does
# Base.Float32(a::T) where {T<:ForwardDiff.Dual} = Float32(ForwardDiff.value(a))

# ============================================================================ #
# ====                              Trace                                 ==== #
# ============================================================================ #

"""
    trace(system::AbstractOpticalSystem{T}, ray::OpticalRay{T}; trackrays = nothing, test = false)

Traces `system` with `ray`, if `test` is enabled then fresnel reflections are disabled and the power distribution will
not be correct. Returns either a [`LensTrace`](@ref) if the ray hits the detector or `nothing` otherwise.

`trackrays` can be passed an empty vector to accumulate the `LensTrace` objects at each intersection of `ray` with a
surface in the system.
"""
function trace(
    system::CSGOpticalSystem{T, D},
    r::OpticalRay{T, N};
    trackrays::Union{Nothing, Vector{LensTrace{T, N}}} = nothing,
    test::Bool = false,
) where {T <: Real, N, D <: Number}
    if power(r) < POWER_THRESHOLD
        return nothing
    end

    result = trace(
        system.assembly,
        r,
        temperature(system),
        pressure(system),
        trackrays = trackrays,
        test = test,
    )

    if result === nothing || result === nopower
        emptyintervalpool!(T)
        return nothing
    else #ray intersected lens assembly so continue to see if ray intersects detector
        intsct = surfaceintersection(detector(system), ray(result))
        if intsct === nothing # no intersection of final ray with detector
            emptyintervalpool!(T)
            return nothing
        end

        detintsct = closestintersection(intsct)
        if detintsct === nothing
            emptyintervalpool!(T)
            return nothing
        else
            # need to modify power and path length accordingly for the intersection with the detector
            surfintsct = point(detintsct)
            nml = normal(detintsct)
            opticalinterface = interface(detintsct)::FresnelInterface{T}
            λ = wavelength(r)

            # Optical path length is measured in mm
            # in this case the result ray is exact so no correction for RAY_OFFSET is needed
            geometricpathlength = norm(surfintsct - origin(ray(result)))
            opticalpathlength = geometricpathlength
            pow = power(result)

            m = outsidematerialid(opticalinterface)
            # compute updated power based on absorption coefficient of material using Beer's law
            # this will almost always not apply as the detector will be in air, but it's possible that the detector is
            # not in air, in which case this is necessary
            if !isair(m)
                mat::Glass = glassforid(m)
                nᵢ = index(
                    mat,
                    λ,
                    temperature = temperature(system),
                    pressure = pressure(system),
                )::T
                α = absorption(
                    mat,
                    λ,
                    temperature = temperature(system),
                    pressure = pressure(system),
                )::T
                if α > zero(T)
                    internal_trans = exp(-α * geometricpathlength)
                    if rand() >= internal_trans
                        return nothing
                    end
                    pow = pow * internal_trans
                end
                opticalpathlength = nᵢ * geometricpathlength
            end

            temp = LensTrace{T, N}(
                OpticalRay(
                    ray(ray(result)),
                    pow,
                    wavelength(result),
                    opl = pathlength(result) + opticalpathlength,
                    nhits = nhits(result) + 1,
                    sourcenum = sourcenum(r),
                    sourcepower = sourcepower(r),
                ),
                detintsct,
            )
            if trackrays !== nothing
                push!(trackrays, temp)
            end

            # increment the detector image
            pixu, pixv = uvtopix(
                detector(system),
                uv(detintsct),
                size(system.detectorimage),
            )
            system.detectorimage[pixv, pixu] += convert(D, sourcepower(r)) # TODO will need to handle different detector
            #      image types a bit better than this

            # should be okay to assume intersection will not be a DisjointUnion for all the types of detectors we will
            # be using
            emptyintervalpool!(T)
            return temp
        end
    end
end


text:
function createnodedf(nodes)

  #=
  mutable struct Node
    bharadius::Float64
    weight::Float64
    fc::Float64
  end
  =#

  node_types = [Float64, Float64, Float64]
  node_names = [:bharadius, :weight, :fc]

  nodedf = DataFrame(node_types, node_names, 0)
  for i in 1:size(nodes, 1)
      append!(nodedf,
        DataFrame(
          bharadius=nodes[i].bharadius,
          weight=nodes[i].weight,
          fc=nodes[i].fc
        )
      )
  end

  nodedf

end