text:
 using RvSpectML
 using DataFrames, Query
 using Statistics
 # order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
 using Dates

make_plots = true
#include("neid_solar_1_read.jl")
# order_list_timeseries = RvSpectML.make_order_list_timeseries(solar_data)
# RvSpectML.normalize_spectra!(order_list_timeseries,solar_data);

espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks","G2.espresso.mas")
  espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
  #lambda_range_with_data = (min = maximum(d->minimum(d.λ),solar_data), max = minimum(d->maximum(d.λ),solar_data) )
  line_list_df = espresso_df |>
      #@filter(lambda_range_with_data.min <= _.lambda ) |>
      #@filter( _.lambda < lambda_range_with_data.max) |>
      DataFrame

inst = RvSpectML.TheoreticalInstrument.TheoreticalInstrument1D()
 spec = RvSpectML.TheoreticalInstrument.generate_spectrum(line_list_df, inst)
 lambda_range_with_data = minimum(spec.λ):maximum(spec.λ)
 line_list_df = espresso_df |>
     @filter(minimum(lambda_range_with_data) <= _.lambda ) |>
     @filter( _.lambda < maximum(lambda_range_with_data) ) |>
     DataFrame

period = 30
 times = range(0.0, stop=0.5*period, length=9)
 rv_mean = -825
 Δ_rvs_true = 3.0 .* cos.(2π.*times./period)
 rvs_true = rv_mean .+ Δ_rvs_true
 @time spectra = RvSpectML.TheoreticalInstrument.generate_spectra_timeseries(times, line_list_df, inst, rvs_true, snr_per_pixel=1000)

lambda_range_with_data = (min = maximum(d->minimum(d.λ),spectra), max = minimum(d->maximum(d.λ),spectra) )
 espresso_mask_df = RvSpectML.read_mask_espresso(espresso_filename)
 chunk_list_df = espresso_mask_df |>
  @filter(lambda_range_with_data.min <= _.lambda ) |>
  @filter( _.lambda < lambda_range_with_data.max) |>
  DataFrame

find_overlapping_chunks(chunk_list_df)
chunk_list_df = RvSpectML.merge_chunks(chunk_list_df)
  #@assert find_overlapping_chunks(chunk_list_df) == nothing

order_list_timeseries = RvSpectML.make_chunk_list_timeseries(spectra, chunk_list_df)

order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
#RvSpectML.normalize_spectra!(order_list_timeseries,spectra);

# Setup to run CCF
mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=1.8) # 1.6)
  line_list = RvSpectML.CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
  ccf_plan = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list)
  v_grid = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan)

# Compute CCF's & measure RVs
tstart = now()
 ccfs = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
 println("# CCF runtime: ", now()-tstart)

#plot(v_grid,ccfs)
using Plots
rvs_ccf_gauss = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i],fit_type = "gaussian") for i in 1:length(order_list_timeseries) ]
 plot(times, rvs_ccf_gauss.-mean(rvs_ccf_gauss))
 rvs_ccf_quad = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i],fit_type = "quadratic") for i in 1:length(order_list_timeseries) ]
 println(" # RV mean: ", mean(rvs_ccf_gauss), "  max-min: ", maximum(rvs_ccf_gauss)-minimum(rvs_ccf_gauss) )
 plot!(times, rvs_ccf_quad.-mean(rvs_ccf_quad))

# Store estimated RVs in metadata
rvs_comp = rvs_ccf_gauss
 oversample_factor = 1
 oversample_fac_orders = 1
 map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_comp[i]-mean(rvs_comp), 1:length(order_list_timeseries) )

@time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2 )  = RvSpectML.make_template_spectra(order_list_timeseries)

order_grids = map(c->RvSpectML.make_grid_for_chunk(order_list_timeseries,c,oversample_factor=oversample_fac_orders, remove_rv_est=false), 1:num_chunks(order_list_timeseries) )
 start = now()
 @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2 )  = RvSpectML.pack_chunk_list_timeseries_to_matrix(order_list_timeseries,order_grids, alg=:TemporalGP ) # :Linear)
 println("# Pack into matrix runtime: ", now()-tstart)
 order_grids_2 = map(c->RvSpectML.make_grid_for_chunk(order_list_timeseries,c,oversample_factor=oversample_fac_orders, remove_rv_est=false), 1:num_chunks(order_list_timeseries) )
 @time ( spectral_orders_matrix_2, f_mean_2, var_mean_2, deriv_2, deriv2_2 ) = RvSpectML.pack_chunk_list_timeseries_to_matrix(order_list_timeseries,order_grids_2, alg=:Linear) # TemporalGP, smooth_factor=1.0)
 #order_grids_2 = map(c->RvSpectML.make_grid_for_chunk(order_list_timeseries,c,oversample_factor=oversample_fac_orders, remove_rv_est=true), 1:num_chunks(order_list_timeseries) )
 #@time spectral_orders_matrix_2 = RvSpectML.pack_shifted_chunk_list_timeseries_to_matrix(order_list_timeseries,order_grids_2, alg=:Linear)
 order_grids_3 = map(c->RvSpectML.make_grid_for_chunk(order_list_timeseries,c,oversample_factor=oversample_fac_orders, remove_rv_est=true), 1:num_chunks(order_list_timeseries) )
 @time ( spectral_orders_matrix_3, f_mean_3, var_mean_3, deriv_3, deriv2_3 ) = RvSpectML.pack_shifted_chunk_list_timeseries_to_matrix(order_list_timeseries,order_grids_3, alg=:TemporalGP, smooth_factor=4.0)

 #f_mean = calc_mean_spectrum(spectral_orders_matrix.flux,spectral_orders_matrix.var)
 #deriv = calc_mean_dfluxdlnlambda(spectral_orders_matrix.flux,spectral_orders_matrix.var,spectral_orders_matrix.λ,spectral_orders_matrix.chunk_map)
(rvs_1, σ_rvs_1) = RvSpectML.calc_rvs_from_taylor_expansion(spectral_orders_matrix,mean=f_mean,deriv=deriv)
 #f_mean_2 = calc_mean_spectrum(spectral_orders_matrix_2.flux,spectral_orders_matrix_2.var)
 #deriv_2 = calc_mean_dfluxdlnlambda(spectral_orders_matrix_2.flux,spectral_orders_matrix_2.var,spectral_orders_matrix_2.λ,spectral_orders_matrix_2.chunk_map)
(rvs_2, σ_rvs_2) = RvSpectML.calc_rvs_from_taylor_expansion(spectral_orders_matrix_2,mean=f_mean_2,deriv=deriv_2)
 #f_mean_3 = calc_mean_spectrum(spectral_orders_matrix_3.flux,spectral_orders_matrix_3.var)
 #deriv_3 = calc_mean_dfluxdlnlambda(spectral_orders_matrix_3.flux,spectral_orders_matrix_3.var,spectral_orders_matrix_3.λ,spectral_orders_matrix_3.chunk_map)
 #deriv2_3 = calc_mean_d2fluxdlnlambda2(spectral_orders_matrix_3.flux,spectral_orders_matrix_3.var,spectral_orders_matrix_3.λ,spectral_orders_matrix_3.chunk_map)
 (rvs_3, σ_rvs_3) = RvSpectML.calc_rvs_from_taylor_expansion(spectral_orders_matrix_3,mean=f_mean_3,deriv=deriv_3)#,deriv2=deriv2_3)
 #println("rvs_3 = [ ", rvs_3[1:3], ", ... ", "σ_rvs_3 = [ ", σ_rvs_3[1:3], ", ... ", )
 flush(stdout)
 println("oversample factor: ", oversample_fac_orders, "  RMS RVs_1: ", std(rvs_1), "   RMS RV_2: ", std(rvs_2), "  RMS RV_3: ", std(rvs_3) )
 println("oversample factor: ", oversample_fac_orders,  "  RMS RVs_1: ", std(rvs_1.-mean(rvs_1).-Δ_rvs_true),
                                                        "   RMS RV_2: ", std(rvs_2.-mean(rvs_2).-Δ_rvs_true),
                                                        "  RMS RV_3: ", std(rvs_3.-mean(rvs_3).-0.0.*Δ_rvs_true) )

plot(rvs_1.-mean(rvs_1).-0.0.*Δ_rvs_true, label="Δrv 1")

plot!((rvs_2.-mean(rvs_2)).-Δ_rvs_true, label="Δrv 2")
 plot!(rvs_3.-mean(rvs_3).-0.0.*Δ_rvs_true, label="Δrv 3")
 plot!((rvs_ccf_gauss.-mean(rvs_ccf_gauss)).-Δ_rvs_true, label="Δrv CCF gauss")
 #plot!((rvs_ccf_quad.-mean(rvs_ccf_quad)).-Δ_rvs_true, label="Δrv CCF quad")

idx_plot = 1000:1500
plot(f_mean[idx_plot])
plot(f_mean_2[idx_plot].-f_mean[idx_plot])
scatter!(f_mean_3[idx_plot].-f_mean[idx_plot], markersize=1.5)

scatter(deriv_2[idx_plot],deriv[idx_plot].-deriv_2[idx_plot], markersize=1.5)
scatter(deriv2_2[idx_plot],deriv2[idx_plot], markersize=1.5)
plot!(deriv_3[idx_plot].-deriv[idx_plot], markersize=1.5)

length(f_mean), length(f_mean_2), length(f_mean_3)

#plot!(f_mean_3)
f_mean_2
scatter(times,Δ_rvs_true,label="RV_true")
 plot!(times,rvs_1.-mean(rvs_1),label="RVs 1")
 scatter!(times,rvs_2.-mean(rvs_2),label="RVs 2")
 plot!(times,rvs_3.-mean(rvs_3),label="RVs 3")

using Plots
scatter(times,rvs_1.+Δ_rvs_true,label="RV_comp - RV_true")

plot()
 scatter!(Δ_rvs_true,rvs_comp.-mean(rvs_comp).-Δ_rvs_true)
 plot!(Δ_rvs_true,rvs_3, label="Δrv")

plot()
 xlims!(1,500)
 #plot((f_mean_3.-mean(f_mean_3))./std(f_mean_3), label="mean")
 #plot!(deriv_3./std(deriv_3), label="dfdlnλ")
 #plot!(deriv2_3./std(deriv2_3), label="d2fdlnλ2")

deriv2_3

plot(spectral_orders_matrix_2.λ[1:300],spectral_orders_matrix_2.flux[1:300,1],markersize=1.5)
plot!(spectral_orders_matrix.λ[1:300],spectral_orders_matrix.flux[1:300,1],markersize=1.5)

spectral_orders_matrix_2.flux
spectral_orders_matrix_2.λ[1:1000]

using Plots
size(spectral_orders_matrix.flux)
idx_plt = 52000:53000
idx_t = 5
plot(RvSpectML.get_λs(order_grids,idx_plt),spectral_orders_matrix.flux[idx_plt,idx_t])
plot!(RvSpectML.get_λs(order_grids_2,idx_plt),spectral_orders_matrix_2.flux[idx_plt,idx_t])
0
maximum(abs.(RvSpectML.get_λs(order_grids,idx_plt).-RvSpectML.get_λs(order_grids_2,idx_plt)))*RvSpectML.speed_of_light_mps
0
plot(RvSpectML.get_λs(order_grids,idx_plt),(RvSpectML.get_λs(order_grids_2,idx_plt).-RvSpectML.get_λs(order_grids,idx_plt))./

spectral_orders_matrix.flux

0
#=

 if 2 <= num_spectra_to_bin <=20
     spectral_orders_matrix = RvSpectML.bin_consecutive_spectra(spectral_orders_matrix,num_spectra_to_bin)
   end
   f_mean_orders = calc_mean_spectrum(spectral_orders_matrix.flux,spectral_orders_matrix.var)
   deriv_orders = calc_mean_dfluxdlnlambda(spectral_orders_matrix.flux,spectral_orders_matrix.var,spectral_orders_matrix.λ,spectral_orders_matrix.chunk_map)

 if make_plots
   plt_order = 40
   plt_order_pix = 3501:3800
   idx_plt = spectral_orders_matrix.chunk_map[plt_order][plt_order_pix]
   plt_λ = RvSpectML.get_λs(order_grids, idx_plt)
      mean_in_plt = mean(f_mean_orders[idx_plt])
      std_in_plt = stdm(f_mean_orders[idx_plt],mean_in_plt)
      local plt = plot()
      plot!(plt,plt_λ,(f_mean_orders[idx_plt].-mean_in_plt)./std_in_plt,label=:none,linecolor=:black)
      plot!(plt,plt_λ, deriv_orders[idx_plt]./std(deriv_orders[idx_plt]),label=:none, linecolor=:green)
      #scatter!(plt,plt_λ,(spectral_orders_matrix.flux[idx_plt,:].-mean_in_plt)./std_in_plt,markersize=1, label=:none)
      xlabel!(plt,"λ (Å)")
      ylabel!(plt,"Mean & Deriv")
      local plt3 = scatter(plt_λ,spectral_orders_matrix.flux[idx_plt,:].-f_mean_orders[idx_plt],markersize=1, label=:none)
      ylabel!(plt3,"Residuals")
      local pltall = plot(plt,plt3,layout=(2,1))
      display(pltall)
 end

 println("Computing RVs using dflux/dlnlambda from orders.")
   order_rvs = RvSpectML.calc_chunk_rvs_from_taylor_expansion(spectral_orders_matrix,mean=f_mean_orders,deriv=deriv_orders)
   ave_order_rvs = vec( sum(mapreduce(c ->order_rvs[c].rv./order_rvs[c].σ_rv.^2, hcat, 1:length(order_rvs)),dims=2) ./
                      sum(mapreduce(c ->(1.0./order_rvs[c].σ_rv.^2), hcat, 1:length(order_rvs)),dims=2) )
   rms_order_rvs = map(order->std(order.rv), order_rvs)
   mean_order_σrvs = map(order->mean(order.σ_rv), order_rvs)
   flush(stdout)
   println("# rms(ave_order_rvs)/√N = ", std(ave_order_rvs)/sqrt(length(ave_order_rvs)), "  <RMS RVs> = ", mean(rms_order_rvs)#=/sqrt(length(ave_order_rvs))=#, "  <σ_RVs> = ", mean(mean_order_σrvs) )
   if make_plots
     plt4 = histogram(abs.(ave_order_rvs),bins=40,label="|Ave Order RVs|", alpha=0.75)
     histogram!(plt4,rms_order_rvs,bins=40,label="RMS Order RVs", alpha=0.75)
     histogram!(plt4,mean_order_σrvs,bins=40,label="σ Order RVs", alpha=0.75)
     xlabel!("(m/s)")
     ylabel!("Counts")
   end

 chunk_rms_cut_off = quantile(rms_order_rvs,0.9)
   idx_good_chunks = findall(x-> x<= chunk_rms_cut_off, rms_order_rvs)
   @assert length(idx_good_chunks)>=1
   ave_good_chunks_rvs = vec( sum(mapreduce(c ->order_rvs[c].rv./order_rvs[c].σ_rv.^2, hcat, idx_good_chunks),dims=2) ./
                       sum(mapreduce(c ->(1.0./order_rvs[c].σ_rv.^2), hcat, idx_good_chunks),dims=2) )
   sigma_good_chunks_rvs = sqrt.(vec( sum(mapreduce(c ->1.0./order_rvs[c].σ_rv.^2, hcat, idx_good_chunks),dims=2) ./
                       sum(mapreduce(c ->(1.0./order_rvs[c].σ_rv.^4), hcat, idx_good_chunks),dims=2) ))

  rms_rvs_ave_good_orders = std(ave_good_chunks_rvs)
  mean_sigma_good_chunks = mean(sigma_good_chunks_rvs)
  flush(stdout)
  println("# rms(RVs_good_orders) = ", rms_rvs_ave_good_orders,  "  <σ_RV good orders> = ", mean_sigma_good_chunks, "   N_obs = ", length(ave_good_chunks_rvs) )
  #map(o->plot!(plt,chunk_list_timeseries.times,order_rvs[o].rv,yerr=order_rvs[o].σ_rv,markersize=1,label=:none),idx_good_chunks )

 if make_plots
    plt = plot()
    scatter!(plt,plt_times, rvs_1, yerr=σ_rvs_1, label="Equal weighted chunks", markersize=3, color=:blue, legend=:topleft)
    scatter!(plt,plt_times, ave_order_rvs,label="Ave Order RVs", markersize=3, color=:green)
    scatter!(plt,plt_times, ave_good_chunks_rvs, yerr=sigma_good_chunks_rvs, label="Ave good chunks", markersize=3, color=:red)
    xlabel!("Time (hours)")
    ylabel!("RV (m/s)")
    display(plt)
 end
=#


text:
using Documenter, Devito

makedocs(sitename="Devito", modules=[Devito])

deploydocs(
    repo = "github.com/ChevronETC/Devito.jl.git",
)

ext: jl


text:
struct Epoch{S<:TimeScale, T} <: Dates.AbstractDateTime
    scale::S
    second::Int64
    fraction::T
    error::T
    function Epoch{S}(second::Int64, fraction::T, error::T=zero(T)) where {S<:TimeScale, T<:AbstractFloat}
        return new{S, T}(S(), second, fraction, error)
    end
end

Epoch{S,T}(ep::Epoch{S,T}) where {S,T} = ep

function Epoch{S}(ep::Epoch{S}, Δt) where {S<:TimeScale}
    second, fraction, error = apply_offset(ep.second, ep.fraction, ep.error, Δt)
    return Epoch{S}(second, fraction, error)
end

Base.show(io::IO, ep::Epoch) = print(io, DateTime(ep), " ", timescale(ep))

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
#! format: off

using SIIPExamples

pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(pkgpath, "test", "2_PowerSystems_examples", "02_parse_matpower.jl"))

ybus = Ybus(sys)

ptdf = PTDF(sys)

lodf = LODF(sys)

ptdf["bus3-bus4-i_6", 3]

buses = collect(get_components(Bus, sys))
ybus[buses[1], buses[2]]

busname2num = get_components(Bus, sys) |> (c -> Dict(zip(get_name.(c), get_number.(c))))
ptdf["bus3-bus4-i_6", busname2num["bus3"]]

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


text:
function force_field(sys::System,field::NoField,comm::MPI.Comm)
end

function clear_field!(sys::System)
    for i in 1:sys.atoms.num_atomtype[1]
        fill!(sys.field.local_field[i] , 0.0)
        fill!(sys.field.global_field[i] , 0.0)
        fill!(sys.field.field_gradient[i] , 0.0)
    end
end

function SerializeCellIndex(i::Int64,j::Int64,k::Int64,celllength::Int64)
    return 1+(i-1)+(j-1)*celllength+(k-1)*celllength*celllength
end

function SerializeCellGradIndex(i::Int64,j::Int64,k::Int64,axis::Int64,celllength::Int64)
    return 1+(i-1)+(j-1)*celllength+(k-1)*celllength*celllength+(axis-1)*celllength*celllength*celllength
end

function particle2cell_vertex(cell_length::Float64,delta_r_x::Float64,delta_r_y::Float64,delta_r_z::Float64,index::Int64)
    if index==1
        return (cell_length-delta_r_x) * (cell_length-delta_r_y) * (cell_length-delta_r_z)
    elseif index==2
        return (cell_length-delta_r_x) * (delta_r_y) * (cell_length-delta_r_z)
    elseif index==3
        return (cell_length-delta_r_x) * (delta_r_y) * (delta_r_z)
    elseif index==4
        return (cell_length-delta_r_x) * (cell_length-delta_r_y) * (delta_r_z)
    elseif index==5
        return (delta_r_x) * (cell_length-delta_r_y) * (cell_length-delta_r_z)
    elseif index ==6
        return (delta_r_x) * (cell_length-delta_r_y) * (delta_r_z)
    elseif index==7
        return (delta_r_x) * (delta_r_y) * (cell_length-delta_r_z)
    elseif index ==8
        return (delta_r_x) * (delta_r_y) * (delta_r_z)
    end
end

function perparticle2cell!( coords::Array{Float64,1},
                            particle2cell::Array{Float64,1},
                            cell_length::Float64,
                            V_cell::Float64)
    delta_r_x=coords[1] - floor(coords[1]/cell_length) * cell_length
    delta_r_y=coords[2] - floor(coords[2]/cell_length) * cell_length
    delta_r_z=coords[3] - floor(coords[3]/cell_length) * cell_length
    particle2cell[1]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,1)/V_cell
    particle2cell[2]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,2)/V_cell
    particle2cell[3]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,3)/V_cell
    particle2cell[4]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,4)/V_cell
    particle2cell[5]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,5)/V_cell
    particle2cell[6]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,6)/V_cell
    particle2cell[7]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,7)/V_cell
    particle2cell[8]=particle2cell_vertex(cell_length,delta_r_x,delta_r_y,delta_r_z,8)/V_cell
end

function perparticle2density!(  local_field::Array{Array{Float64,1}},
                                types::Int64,
                                cellindex::Array{Int64,2},
                                particle2cell::Array{Float64,1},
                                num_cells::Int64)
    for j in 1:8
        local_field[types][SerializeCellIndex(cellindex[j,1],cellindex[j,2],cellindex[j,3],num_cells)] += particle2cell[j]
    end
end

function pbc_mesh(index::Int64,num_cells::Int64)
    if index>0
        return index-floor(Int64,index/(num_cells-1))*(num_cells-1)
    else
        return index-floor(Int64,index/(num_cells-1))*(num_cells-1)
    end
end

function update_particlecellindex!(coords::Array{Float64,1},cellindex::Array{Int64,2},cell_length::Float64,num_cells::Int64)

    cell_index_x=floor(coords[1] / cell_length) # 0 to num_cells-1, PBC: 0=num_cells-1
    cell_index_y=floor(coords[2] / cell_length) # 0 to num_cells-1, PBC: 0=num_cells-1
    cell_index_z=floor(coords[3] / cell_length) # 0 to num_cells-1, PBC: 0=num_cells-1

    cell_index_x=cell_index_x-floor(Int64,cell_index_x/((num_cells-1)))*(num_cells-1)
    cell_index_y=cell_index_y-floor(Int64,cell_index_y/((num_cells-1)))*(num_cells-1)
    cell_index_z=cell_index_z-floor(Int64,cell_index_z/((num_cells-1)))*(num_cells-1)

    cell_index_plus_x=cell_index_x +1-floor(Int64,(cell_index_x+1)/((num_cells-1)))*(num_cells-1)
    cell_index_plus_y=cell_index_y +1-floor(Int64,(cell_index_y+1)/((num_cells-1)))*(num_cells-1)
    cell_index_plus_z=cell_index_z +1-floor(Int64,(cell_index_z+1)/((num_cells-1)))*(num_cells-1)

    cell_index_x=cell_index_x+1
    cell_index_y=cell_index_y+1
    cell_index_z=cell_index_z+1
    cell_index_plus_x=cell_index_plus_x+1
    cell_index_plus_y=cell_index_plus_y+1
    cell_index_plus_z=cell_index_plus_z+1

    cellindex[1,1]=cell_index_x
    cellindex[1,2]=cell_index_y
    cellindex[1,3]=cell_index_z

    cellindex[2,1]=cell_index_x
    cellindex[2,2]=cell_index_plus_y
    cellindex[2,3]=cell_index_z

    cellindex[3,1]=cell_index_x
    cellindex[3,2]=cell_index_plus_y
    cellindex[3,3]=cell_index_plus_z

    cellindex[4,1]=cell_index_x
    cellindex[4,2]=cell_index_y
    cellindex[4,3]=cell_index_plus_z

    cellindex[5,1]=cell_index_plus_x
    cellindex[5,2]=cell_index_y
    cellindex[5,3]=cell_index_z

    cellindex[6,1]=cell_index_plus_x
    cellindex[6,2]=cell_index_y
    cellindex[6,3]=cell_index_plus_z

    cellindex[7,1]=cell_index_plus_x
    cellindex[7,2]=cell_index_plus_y
    cellindex[7,3]=cell_index_z

    cellindex[8,1]=cell_index_plus_x
    cellindex[8,2]=cell_index_plus_y
    cellindex[8,3]=cell_index_plus_z
end

function grad_field_vertex!(coord_cell::Array{Int64,1},
                            global_field::Array{Float64,1},
                            field_gradient::Array{Float64,1},
                            num_cells::Int64)
    icell_plus=pbc_mesh(coord_cell[1]+1,num_cells)+1
    jcell_plus=pbc_mesh(coord_cell[2]+1,num_cells)+1
    kcell_plus=pbc_mesh(coord_cell[3]+1,num_cells)+1

    icell_minus=pbc_mesh(coord_cell[1]-1,num_cells)+1
    jcell_minus=pbc_mesh(coord_cell[2]-1,num_cells)+1
    kcell_minus=pbc_mesh(coord_cell[3]-1,num_cells)+1

    icell_=coord_cell[1]+1
    jcell_=coord_cell[2]+1
    kcell_=coord_cell[3]+1

    g_x=global_field[SerializeCellIndex(icell_plus,jcell_,kcell_,num_cells)]-global_field[SerializeCellIndex(icell_minus,jcell_,kcell_,num_cells)]
    g_y=global_field[SerializeCellIndex(icell_,jcell_plus,kcell_,num_cells)]-global_field[SerializeCellIndex(icell_,jcell_minus,kcell_,num_cells)]
    g_z=global_field[SerializeCellIndex(icell_,jcell_,kcell_plus,num_cells)]-global_field[SerializeCellIndex(icell_,jcell_,kcell_minus,num_cells)]

    field_gradient[SerializeCellGradIndex(icell_,jcell_,kcell_,1,num_cells)] += 0.5*g_x
    field_gradient[SerializeCellGradIndex(icell_,jcell_,kcell_,2,num_cells)] += 0.5*g_y
    field_gradient[SerializeCellGradIndex(icell_,jcell_,kcell_,3,num_cells)] += 0.5*g_z
end

function force_field_particle!(forces::Array{Float64,1},
                                types::Int64,
                                particle2cell::Array{Float64,1},
                                cellindex::Array{Int64,2},
                                num_atomtype::Int64,
                                num_cells::Int64,
                                field::Field,
                                field_gradient::Array{Array{Float64,1},1})
    sum_graident_field_x=0.0
    sum_graident_field_y=0.0
    sum_graident_field_z=0.0
    for atomtype in 1:num_atomtype
        field_gradient_x=0.0
        field_gradient_y=0.0
        field_gradient_z=0.0
        for j in 1:8
            field_gradient_x+=field_gradient[atomtype][SerializeCellGradIndex(cellindex[j,1],cellindex[j,2],cellindex[j,3],1,num_cells)] * particle2cell[j]
            field_gradient_y+=field_gradient[atomtype][SerializeCellGradIndex(cellindex[j,1],cellindex[j,2],cellindex[j,3],2,num_cells)] * particle2cell[j]
            field_gradient_z+=field_gradient[atomtype][SerializeCellGradIndex(cellindex[j,1],cellindex[j,2],cellindex[j,3],3,num_cells)] * particle2cell[j]
        end
        forces[1]+=-field.χ[types,atomtype]*field_gradient_x
        forces[2]+=-field.χ[types,atomtype]*field_gradient_y
        forces[3]+=-field.χ[types,atomtype]*field_gradient_z
        sum_graident_field_x+=field_gradient_x
        sum_graident_field_y+=field_gradient_y
        sum_graident_field_z+=field_gradient_z
    end
    forces[1] += -1.0/field.κ * sum_graident_field_x
    forces[2] += -1.0/field.κ * sum_graident_field_y
    forces[3] += -1.0/field.κ * sum_graident_field_z
end

function force_field!(sys::System,field::CanonicalField,comm::MPI.Comm)

    if sys.first_step[1] || sys.current_step[1] % field.τ==0
        clear_field!(sys)
        perparticle2cell!.(sys.atoms.coords,sys.field.particle2cell,sys.field.mesh.cell_length,sys.field.mesh.V_cell)
        update_particlecellindex!.(sys.atoms.coords,sys.field.cellindex,sys.field.mesh.cell_length,sys.field.mesh.num_cells)
        perparticle2density!.(Ref(sys.field.local_field),sys.atoms.types,sys.field.cellindex,sys.field.particle2cell,sys.field.mesh.num_cells)
        for i in 1:sys.atoms.num_atomtype[1]
            MPI.Allreduce!(sys.field.local_field[i],sys.field.global_field[i],sum_mesh,comm)
        end
        for i in 1:sys.atoms.num_atomtype[1]
            grad_field_vertex!.(sys.field.mesh.mesh_index,Ref(sys.field.global_field[i]),Ref(sys.field.field_gradient[i]),sys.field.mesh.num_cells)
        end
    else
        perparticle2cell!.( sys.atoms.coords,
                            sys.field.particle2cell,
                            sys.field.mesh.cell_length,
                            sys.field.mesh.V_cell)

        update_particlecellindex!.( sys.atoms.coords,
                                    sys.field.cellindex,
                                    sys.field.mesh.cell_length,
                                    sys.field.mesh.num_cells)
    end

    force_field_particle!.( sys.atoms.forces,
                            sys.atoms.types,
                            sys.field.particle2cell,
                            sys.field.cellindex,
                            sys.atoms.num_atomtype[1],
                            sys.field.mesh.num_cells,
                            sys.field,
                            Ref(sys.field.field_gradient))

end


text:
using ConvDiffMIPDECO
using jInv.Mesh
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.LinearSolvers
using MUMPSjInv
using MAT

filename= joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples","Sources3D.mat")
file = matread(filename)
wtrue = vec(file["W"]) ## Sources on finest mesh
domain = file["domain"]
m      = vec(file["m"])

M      = getRegularMesh(domain,m)

rec = rand(100,2) .*[domain[2] domain[4] ]

Ainv = getMUMPSsolver()

# build linear interpolation matrix from nodes to receiver locations
x1c,x2c,x3c = getCellCenteredAxes(M)
rec3D = [kron(ones(length(x3c)),rec) kron(x3c,ones(size(rec,1)))]

x1,x2,x3 = getNodalAxes(M)
P = interpmat(x1,x2,x3, rec3D);

# build param that holds forward problem
pFor = getConvDiffFEMParam(M,Ainv=Ainv,P=P);
dtrue, pFor = getData(vec(wtrue),pFor);
utrue = pFor.Fields

file["utrue"]  = reshape(utrue,tuple(M.n.+1...));
file["dtrue"]  = dtrue
file["rec"]    = rec
file["sig"]    = pFor.sig
file["v"]      = pFor.v

matwrite(filename,file)

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
abstract type ConicSection end

struct Circle <: ConicSection
    radius
    center
end

struct Ellipse <: ConicSection
end

struct Parabola <: ConicSection
end

struct Hyperbola <: ConicSection
end


text:
# ------------ Vectorized utilities ------------ #

function keyset(sel::L, n_len::Int) where L <: ConstrainedSelection
    keyset = Set{Int}()
    min = n_len
    for (k, v) in sel.tree
        k isa Int && k > 0 && k <= n_len && !isempty(v) && begin
            push!(keyset, k)
            if k < min
                min = k
            end
        end
    end
    return min, keyset
end

function keyset(sel::ConstrainedEmptySelection, n_len::Int)
    keyset = Set{Int}()
    return n_len, keyset
end

function keyset(sel::L, n_len::Int) where L <: UnconstrainedSelection
    keyset = Set{Int}()
    min = n_len
    for (k, _) in sel.tree
        k isa Int && k > 0 && k <= n_len && begin
            push!(keyset, k)
            if k < min
                min = k
            end
        end
    end
    return min, keyset
end

function keyset(sel::UnconstrainedEmptySelection, n_len::Int)
    keyset = Set{Int}()
    return n_len, keyset
end


text:
# backward pass of the SDDP algorithm
function fBuild_D(td, ωd, currentPath, τ, Δt, T, qpopt = false, solveOpt = true, hardened = [])
    prevtpInd = maximum([i for i in 1:length(currentPath) if currentPath[i][2] < td]);
    currentSol = currentPath[prevtpInd][1];
    # precalculate data
    Rdict = Dict();
    Xdict = Dict();
    for k in fData.brList
        Rdict[k] = fData.g[k]/(fData.g[k]^2 + fData.b[k]^2);
        Xdict[k] = -fData.b[k]/(fData.g[k]^2 + fData.b[k]^2);
    end
    Ω = [ω for ω in keys(pDistr.ωDistrn)];
    Bparams = Dict();
    for t in td:T
        # create B parameters
        for k in fData.brList
            # if the line is disrupted and it is within disruption time
            if (((k[1],k[2]) in ωd)|((k[2],k[1]) in ωd))&(t <= td + τ)
                if ((k[1],k[2]) in hardened) | ((k[2],k[1]) in hardened)
                    Bparams[k,t] = 1;
                else
                    Bparams[k,t] = 0;
                end
            else
                Bparams[k,t] = 1;
            end
        end
        for i in fData.genIDList
            if (i in ωd)&(t <= td + τ)
                if i in hardened
                    Bparams[i,t] = 1;
                else
                    Bparams[i,t] = 0;
                end
            else
                Bparams[i,t] = 1;
            end
        end
    end
    for i in fData.genIDList
        Bparams[i,td - 1] = 1;
    end

    if qpopt
        mp = Model(solver = IpoptSolver(print_level = 0, linear_solver = "ma27"));
    else
        mp = Model(solver = GurobiSolver(GUROBI_ENV, OutputFlag = 0, QCPDual = 1,Threads = 1,NumericFocus = 3));
    end

    # set up the variables
    @variable(mp, fData.Pmin[i]*Bparams[i,t] <= sp[i in fData.genIDList,t in (td - 1):T] <= fData.Pmax[i]*Bparams[i,t]);
    @variable(mp, fData.Qmin[i]*Bparams[i,t] <= sq[i in fData.genIDList,t in td:T] <= fData.Qmax[i]*Bparams[i,t]);
    sphatsum = Dict();
    for t in td:T
        for i in fData.IDList
            sphatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sphatsum[i,t] += sp[j,t];
                end
            end
        end
    end
    sqhatsum = Dict();
    for t in td:T
        for i in fData.IDList
            sqhatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sqhatsum[i,t] += sq[j,t];
                end
            end
        end
    end

    @variable(mp, p[k in fData.brList, t in td:T]);
    @variable(mp, q[k in fData.brList, t in td:T]);
    @variable(mp, fData.Vmin[i]^2 <= v[i in fData.IDList, t in td:T] <= fData.Vmax[i]^2);
    @variable(mp, 0 <= w[i in bData.IDList, t in (td - 1):T] <= bData.cap[i]);
    @variable(mp, y[i in bData.IDList, t in td:T]);
    @variable(mp, zp[i in bData.IDList, t in td:T]);
    @variable(mp, zq[i in bData.IDList, t in td:T]);
    @variable(mp, lpp[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lqp[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lpm[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lqm[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, 0 <= u[i in bData.IDList] <= bData.uCap[i]);
    @variable(mp, θ[tp in (td + τ + 1):T, ω in Ω] >= 0);

    # set up the constraints
    @constraint(mp, pBalance[i in fData.IDList, t in td:T], sum(zp[b,t] for b in bData.IDList if bData.Loc[b] == i) + lpp[i,t] - lpm[i,t] +
        sphatsum[i,t] - dData.pd[i][t] == sum(p[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, qBalance[i in fData.IDList, t in td:T], sum(zq[b,t] for b in bData.IDList if bData.Loc[b] == i) + lqp[i,t] - lqm[i,t] +
        sqhatsum[i,t] - dData.qd[i][t] == sum(q[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, pequal[k in fData.brList, t in td:T], p[k,t] == -p[(k[2],k[1],k[3]),t]);
    @constraint(mp, qequal[k in fData.brList, t in td:T], q[k,t] == -q[(k[2],k[1],k[3]),t]);
    if qpopt
        @constraint(mp, lineThermal1[k in fData.brList, t in td:T;Bparams[k,t] == 1], p[k,t]^2 + q[k,t]^2 <= fData.rateA[k]^2);
        @constraint(mp, bThermal[i in bData.IDList, t in td:T], zp[i,t]^2 + zq[i,t]^2 <= u[i]^2);
    else
        @constraint(mp, lineThermal1[k in fData.brList, t in td:T;Bparams[k,t] == 1], norm([p[k,t],q[k,t]]) <= fData.rateA[k]);
        @constraint(mp, bThermal[i in bData.IDList, t in td:T], norm([zp[i,t],zq[i,t]]) <= u[i]);
    end
    @constraint(mp, lineThermal2[k in fData.brList, t in td:T;Bparams[k,t] == 0], p[k,t] == 0);
    @constraint(mp, lineThermal3[k in fData.brList, t in td:T;Bparams[k,t] == 0], q[k,t] == 0);
    @constraint(mp, powerflow1[k in fData.brList, t in td:T;Bparams[k,t] == 1], v[k[2],t] == v[k[1],t] - 2*(Rdict[k]*p[k,t] + Xdict[k]*q[k,t]));
    @constraint(mp, rampUp[i in fData.genIDList, t in td:T; Bparams[i,t] == 1], sp[i,t] - sp[i,t - 1] <= fData.RU[i]);
    @constraint(mp, rampDown[i in fData.genIDList, t in td:T; Bparams[i,t] == 1], sp[i,t] - sp[i,t - 1] >= fData.RD[i]);
    @constraint(mp, bInv[i in bData.IDList, t in td:T], w[i,t] == w[i,t-1] - y[i,t]*Δt);
    @constraint(mp, bThermal1[i in bData.IDList, t in td:T], zp[i,t] <= u[i]);
    @constraint(mp, bThermal2[i in bData.IDList, t in td:T], zq[i,t] <= u[i]);
    @constraint(mp, bEfficient[i in bData.IDList, l in 1:length(bData.ηα[i]), t in td:T], zp[i,t] <= bData.ηα[i][l]*y[i,t] + bData.ηβ[i][l]);
    @constraint(mp, bInvIni[i in bData.IDList], w[i,td - 1] == currentSol.w[i,td - 1]);
    @constraint(mp, spIni[i in fData.genIDList], sp[i,td - 1] == currentSol.sp[i,td - 1]);
    @constraint(mp, uIni[i in bData.IDList], u[i] == currentSol.u[i]);

    # set up the cuts, here tp is the disruption time
    for tp in (td + τ + 1):T
        for ω in Ω
            if (tp,ω) in keys(cutDict)
                for l in 1:length(cutDict[tp,ω])
                    @constraint(mp, θ[tp,ω] >= cutDict[tp,ω][l].vhat +
                        sum(cutDict[tp,ω][l].λ[i]*sp[i,tp - 1] for i in fData.genIDList) +
                        sum(cutDict[tp,ω][l].γ[i]*w[i,tp - 1] + cutDict[tp,ω][l].μ[i]*u[i] for i in bData.IDList));
                end
            end
        end
    end

    # set up the objective function
    objExpr = calObj(mp, td, τ, T, fData, bData, pDistr, sp, lpp, lqp, lpm, lqm, θ);
    @objective(mp, Min, objExpr);

    if solveOpt
        # solve the problem
        # optimize!(mp, with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0,
        #     QCPDual = 1, NumericFocus = 3, BarQCPConvTol = 1e-9, FeasibilityTol = 1e-9));
        #optimize!(mp, with_optimizer(Ipopt.Optimizer, linear_solver = "ma27", print_level = 0, acceptable_tol = 1e-8, max_iter = 10000));
        statusMp = solve(mp);
        println(statusMp, " ", td, " ", ωd);
        # obtain the primal solutions & obj value
        vhat = getobjectivevalue(mp);
        # obtain the solutions
        solSp = Dict();
        solw = Dict();
        solu = Dict();
        for i in fData.genIDList
            solSp[i] = getvalue(sp[i,td - 1]);
        end
        for i in bData.IDList
            solu[i] = getvalue(u[i]);
            solw[i] = getvalue(w[i,td - 1]);
        end
        # obtain the dual solutions
        dsolλ = Dict();
        dsolγ = Dict();
        dsolμ = Dict();
        for i in fData.genIDList
            dsolλ[i] = getdual(spIni[i]);
            vhat -= dsolλ[i]*solSp[i];
            if abs(dsolλ[i]) < 1e-4
                if dsolλ[i] < 0
                    vhat += dsolλ[i]*fData.Pmax[i];
                end
                dsolλ[i] = 0;
            end
        end
        for i in bData.IDList
            dsolγ[i] = getdual(bInvIni[i]);
            vhat -= dsolγ[i]*solw[i];
            if abs(dsolγ[i]) < 1e-4
                if dsolγ[i] < 0
                    vhat += dsolγ[i]*bData.cap[i];
                end
                dsolγ[i] = 0;
            end
            dsolμ[i] = getdual(uIni[i]);
            vhat -= dsolμ[i]*solu[i];
            if abs(dsolμ[i]) < 1e-4
                if dsolμ[i] < 0
                    vhat += dsolμ[i]*bData.uCap[i];
                end
                dsolμ[i] = 0;
            end
        end
        cutTemp = cutData(statusMp,dsolλ,dsolγ,dsolμ,vhat);
        return cutTemp;
    else
        return mp;
    end
end

function dfBuild_D(td, ωd, currentPath, τ, Δt, T, qpopt = false, solveOpt = true, hardened = [])
    # trying to find the previous solution that can be used for this td
    prevtpInd = maximum([i for i in 1:length(currentPath) if currentPath[i][2] < td]);
    currentSol = currentPath[prevtpInd][1];
    # precalculate data
    Rdict = Dict();
    Xdict = Dict();
    for k in fData.brList
        Rdict[k] = fData.g[k]/(fData.g[k]^2 + fData.b[k]^2);
        Xdict[k] = -fData.b[k]/(fData.g[k]^2 + fData.b[k]^2);
    end
    Ω = [ω for ω in keys(pDistr.ωDistrn)];
    Bparams = Dict();
    for t in td:T
        # create B parameters
        for k in fData.brList
            # if the line is disrupted and it is within disruption time
            if (((k[1],k[2]) in ωd)|((k[2],k[1]) in ωd))&(t <= td + τ)
                if ((k[1],k[2]) in hardened) | ((k[2],k[1]) in hardened)
                    Bparams[k,t] = 1;
                else
                    Bparams[k,t] = 0;
                end
            else
                Bparams[k,t] = 1;
            end
        end
        for i in fData.genIDList
            if (i in ωd)&(t <= td + τ)
                if i in hardened
                    Bparams[i,t] = 1;
                else
                    Bparams[i,t] = 0;
                end
            else
                Bparams[i,t] = 1;
            end
        end
    end
    for i in fData.genIDList
        Bparams[i,td - 1] = 1;
    end

    if qpopt
        dp = Model(solver = IpoptSolver(print_level = 0, linear_solver = "ma27"));
    else
        dp = Model(solver = GurobiSolver(GUROBI_ENV,OutputFlag = 0,Threads = 1));
    end

    # simple bounds dual variables
    @variable(dp, λspu[i in fData.genIDList,t in td:T] >= 0);
    @variable(dp, λspl[i in fData.genIDList,t in td:T] <= 0);
    @variable(dp, λsqu[i in fData.genIDList,t in td:T] >= 0);
    @variable(dp, λsql[i in fData.genIDList,t in td:T] <= 0);
    @variable(dp, λvu[i in fData.IDList,t in td:T] >= 0);
    @variable(dp, λvl[i in fData.IDList,t in td:T] <= 0);
    @variable(dp, λwu[i in bData.IDList,t in td:T] >= 0);
    @variable(dp, λuu[i in bData.IDList] >= 0);

    # cuts dual variables
    λcuts = Dict();
    spCuts = Dict();
    wCuts = Dict();
    for tp in (td-1):T
        spCuts[tp] = [];
        wCuts[tp] = [];
    end
    for tp in (td + τ + 1):T
        for ω in Ω
            if (tp,ω) in keys(cutDict)
                for l in 1:length(cutDict[tp,ω])
                    λcuts[tp,ω,l] = @variable(dp, upperbound = 0);
                    push!(spCuts[tp - 1],(tp,ω,l));
                    push!(wCuts[tp - 1],(tp,ω,l));
                end
            end
        end
    end

    # linear constraints dual variables
    @variable(dp, λpb[i in fData.IDList, t in td:T]);
    @variable(dp, λqb[i in fData.IDList, t in td:T]);
    @variable(dp, λpe[k in fData.brList, t in td:T]);
    @variable(dp, λqe[k in fData.brList, t in td:T]);
    @variable(dp, λlt2[k in fData.brList, t in td:T]);
    @variable(dp, λlt3[k in fData.brList, t in td:T]);
    @variable(dp, λpf[k in fData.brList, t in td:T]);
    @variable(dp, λru[i in fData.genIDList, t in td:T] >= 0);
    @variable(dp, λrd[i in fData.genIDList, t in td:T] <= 0);
    @variable(dp, λbInv[i in bData.IDList, t in td:T]);
    @variable(dp, λbCoeff[i in bData.IDList, l in 1:length(bData.ηα[i]), t in td:T] >= 0);
    @variable(dp, λsp[i in fData.genIDList]);
    @variable(dp, λw[i in bData.IDList]);
    @variable(dp, λu[i in bData.IDList]);
    @variable(dp, λgc2[i in fData.genIDList, t in td:T]);
    @variable(dp, λgcAux1[i in fData.genIDList, t in td:T]);
    @variable(dp, λgcAux2[i in fData.genIDList, t in td:T]);
    @variable(dp, λgcAux3[i in fData.genIDList, t in td:T]);

    # SOCP constraints dual variables
    @variable(dp,μ1[k in fData.brList, t in td:T, j in 1:2]);
    @variable(dp,μ2[i in bData.IDList, t in td:T, j in 1:2]);
    @variable(dp,μ3[i in fData.genIDList, t in td:T, j in 1:2]);
    @variable(dp,ν1[k in fData.brList, t in td:T] >= 0);
    @variable(dp,ν2[i in bData.IDList, t in td:T] >= 0);
    @variable(dp,ν3[i in fData.genIDList, t in td:T] >= 0);

    # dual constraints:
    # sp
    bool3 = Dict();
    for i in fData.genIDList
        if fData.cp[i].n == 3
            bool3[i] = 1;
        else
            bool3[i] = 0;
        end
    end
    fDict,lDict,θDict = calDualC(td, τ, T, fData, pDistr);
    @constraint(dp, spConstraint[i in fData.genIDList, t in td:(T - 1)],
        λspu[i,t] + λspl[i,t] + λpb[fData.Loc[i],t] + λru[i,t]*Bparams[i,t] -
        λru[i,t+1]*Bparams[i,t+1] + λrd[i,t]*Bparams[i,t] - λrd[i,t+1]*Bparams[i,t+1] -
        λgcAux3[i,t]*sqrt(fData.cp[i].params[1])*bool3[i] - λgc2[i,t]*fData.cp[i].params[1]*(1 - bool3[i]) -
        sum(λcuts[item]*cutDict[(item[1],item[2])][item[3]].λ[i] for item in spCuts[t] if t >= td + τ) == 0);
    @constraint(dp, spConstraintT[i in fData.genIDList],
        λspu[i,T] + λspl[i,T] + λpb[fData.Loc[i],T] + λru[i,T]*Bparams[i,T] + λrd[i,T]*Bparams[i,T] -
        λgcAux3[i,T]*sqrt(fData.cp[i].params[1])*bool3[i] - λgc2[i,T]*fData.cp[i].params[1]*(1 - bool3[i]) == 0);# for all time period after td
    @constraint(dp, spConstraint1[i in fData.genIDList], -λru[i,td]*Bparams[i,td] - λrd[i,td]*Bparams[i,td] + λsp[i] == 0); # for td-1
    # sq
    @constraint(dp, sqConstraint[i in fData.genIDList, t in td:T],
        λsqu[i,t] + λsql[i,t] + λqb[fData.Loc[i],t] == 0);
    # p
    @constraint(dp, pConstraint[k in fData.brList, t in td:T], -λpb[k[1],t] + λpe[k,t] + λpe[(k[2],k[1],k[3]),t] +
        (1 - Bparams[k,t])*λlt2[k,t] + 2*Rdict[k]*λpf[k,t]*Bparams[k,t] - μ1[k,t,1]*Bparams[k,t] == 0);
    # q
    @constraint(dp, qConstraint[k in fData.brList, t in td:T], -λqb[k[1],t] + λqe[k,t] + λqe[(k[2],k[1],k[3]),t] +
        (1 - Bparams[k,t])*λlt3[k,t] + 2*Xdict[k]*λpf[k,t]*Bparams[k,t] - μ1[k,t,2]*Bparams[k,t] == 0);
    # v
    @constraint(dp, vConstraint[i in fData.IDList, t in td:T], λvu[i,t] + λvl[i,t] + sum(λpf[k,t]*Bparams[k,t] for k in fData.branchDict2[i]) +
        sum(-λpf[k,t]*Bparams[k,t] for k in fData.branchDict1[i]) == 0);
    # lpp
    @constraint(dp, lppConstraint[i in fData.IDList, t in td:T], λpb[i,t] + lDict[i,t] >= 0);
    # lpm
    @constraint(dp, lpmConstraint[i in fData.IDList, t in td:T], -λpb[i,t] + lDict[i,t] >= 0);
    # lqp
    @constraint(dp, lqpConstraint[i in fData.IDList, t in td:T], λqb[i,t] + lDict[i,t] >= 0);
    # lqm
    @constraint(dp, lqmConstraint[i in fData.IDList, t in td:T], -λqb[i,t] + lDict[i,t] >= 0);

    # w
    @constraint(dp, wConstraint[i in bData.IDList, t in td:(T - 1)], λwu[i,t] + λbInv[i,t] - λbInv[i,t+1] -
        sum(λcuts[item]*cutDict[(item[1],item[2])][item[3]].γ[i] for item in wCuts[t] if t >= td + τ) >= 0);
    @constraint(dp, wConstraintT[i in bData.IDList], λwu[i,T] + λbInv[i,T] >= 0);# time after td
    @constraint(dp, wConstraint1[i in bData.IDList], -λbInv[i,td] + λw[i] >= 0);   # td-1
    # u
    @constraint(dp, uConstraint[i in bData.IDList], λuu[i] + λu[i] - sum(ν2[i,t] for t in td:T) -
        sum(λcuts[item]*cutDict[(item[1],item[2])][item[3]].μ[i] for item in keys(λcuts) if item[1] >= td + τ + 1) >= 0);
    # zp
    @constraint(dp, zpConstraint1[i in bData.IDList, t in td:T], λpb[bData.Loc[i],t] + sum(λbCoeff[i,l,t] for l in 1:length(bData.ηα[i])) -
        μ2[i,t,1] == 0);
    # zq
    @constraint(dp, zqConstraint1[i in bData.IDList, t in td:T], λqb[bData.Loc[i],t] - μ2[i,t,2] == 0);
    # y
    @constraint(dp, yConstraint1[i in bData.IDList, t in td:T], λbInv[i,t]*Δt - sum(bData.ηα[i][l]*λbCoeff[i,l,t] for l in 1:length(bData.ηα[i])) == 0);

    # θ
    @constraint(dp, θConstraint1[tp in (td + τ + 1):T, ω in Ω; (tp,ω) in keys(cutDict)],
        sum(λcuts[tp,ω,l] for l in 1:length(cutDict[tp,ω])) + θDict[tp,ω] >= 0);

    # f
    @constraint(dp, fConstraint3[i in fData.genIDList, t in td:T;fData.cp[i].n == 3],
        -λgcAux1[i,t] - λgcAux2[i,t] + fDict[i,t] >= 0);
    @constraint(dp, fConstraint2[i in fData.genIDList, t in td:T;fData.cp[i].n == 2], λgc2[i,t] + fDict[i,t] >= 0);

    # tAux1
    @constraint(dp, tAux1Constraint[i in fData.genIDList, t in td:T;fData.cp[i].n == 3],
        -ν3[i,t] + λgcAux1[i,t] == 0);
    # tAux2
    @constraint(dp, tAux2Constraint[i in fData.genIDList, t in td:T;fData.cp[i].n == 3],
        -μ3[i,t,1] + λgcAux2[i,t] == 0);
    # tAux3
    @constraint(dp, tAux3Constraint[i in fData.genIDList, t in td:T;fData.cp[i].n == 3],
        -μ3[i,t,2] + λgcAux3[i,t] == 0);

    # SOCP
    if qpopt
        for t in td:T
            for k in fData.brList
                if Bparams[k,t] == 1
                    @constraint(dp, μ1[k,t,1]^2 + μ1[k,t,2]^2 <= ν1[k,t]^2);
                end
            end
            for i in bData.IDList
                @constraint(dp, μ2[i,t,1]^2 + μ2[i,t,2]^2 <= ν2[i,t]^2);
            end
            for i in fData.genIDList
                if fData.cp[i].n == 3
                    @constraint(dp, μ3[i,t,1]^2 + μ3[i,t,2]^2 <= ν3[i,t]^2);
                end
            end
        end
    else
        for t in td:T
            for k in fData.brList
                if Bparams[k,t] == 1
                    μ1List = [μ1[k,t,1],μ1[k,t,2]];
                    @constraint(dp,norm(μ1List) <= ν1[k,t]);
                end
            end
            for i in bData.IDList
                μ2List = [μ2[i,t,1],μ2[i,t,2]];
                @constraint(dp,norm(μ2List) <= ν2[i,t]);
            end
            for i in fData.genIDList
                if fData.cp[i].n == 3
                    μ3List = [μ3[i,t,1],μ3[i,t,2]];
                    @constraint(dp,norm(μ3List) <= ν3[i,t]);
                end
            end
        end
    end

    # objective function
    @objective(dp, Max, sum(sum(-λspu[i,t]*fData.Pmax[i]*Bparams[i,t] - λspl[i,t]*fData.Pmin[i]*Bparams[i,t] -
            λsqu[i,t]*fData.Qmax[i]*Bparams[i,t] - λsql[i,t]*fData.Qmin[i]*Bparams[i,t] -
            λru[i,t]*fData.RU[i]*Bparams[i,t] - λrd[i,t]*fData.RD[i]*Bparams[i,t] for i in fData.genIDList) for t in td:T) +
            sum(sum(-λvu[i,t]*fData.Vmax[i]^2 - λvl[i,t]*fData.Vmin[i]^2 - λpb[i,t]*dData.pd[i][t] - λqb[i,t]*dData.qd[i][t]
            for i in fData.IDList) for t in td:T) + sum(sum(-ν1[k,t]*fData.rateA[k]*Bparams[k,t] for k in fData.brList) +
            sum(-λwu[i,t]*bData.cap[i] - sum(λbCoeff[i,l,t]*bData.ηβ[i][l] for l in 1:length(bData.ηα[i]))
            for i in bData.IDList) for t in td:T) + sum(-λsp[i]*currentSol.sp[i,td - 1] for i in fData.genIDList) +
            sum(-λuu[i]*bData.uCap[i] - λu[i]*currentSol.u[i] - λw[i]*currentSol.w[i,td - 1] for i in bData.IDList) +
            sum(sum(-λgcAux1[i,t]*((fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1])) -
            λgcAux2[i,t]*((-fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1])) -
            λgcAux3[i,t]*(fData.cp[i].params[2]/(2*sqrt(fData.cp[i].params[1]))) for i in fData.genIDList if fData.cp[i].n == 3) for t in td:T) -
            sum(sum(sum(λcuts[tp,ω,l]*cutDict[tp,ω][l].vhat for l in 1:length(cutDict[tp,ω])) for ω in Ω if (tp,ω) in keys(cutDict)) for tp in (td + τ + 1):T)
            );

    # return the cut or the dual problem
    if solveOpt
        statusDp = solve(dp);
        println(statusDp, " ", td, " ", ωd);
        # obtain the primal solutions & obj value
        vhat = getobjectivevalue(dp);

        solSp = Dict();
        solw = Dict();
        solu = Dict();
        for i in fData.genIDList
            solSp[i] = currentSol.sp[i,td - 1];
        end
        for i in bData.IDList
            solu[i] = currentSol.u[i];
            solw[i] = currentSol.w[i,td - 1];
        end
        # obtain the dual solutions
        dsolλ = Dict();
        dsolγ = Dict();
        dsolμ = Dict();
        for i in fData.genIDList
            dsolλ[i] = -getvalue(λsp[i]);
            vhat -= dsolλ[i]*solSp[i];
            if abs(dsolλ[i]) < 1e-4
                if dsolλ[i] < 0
                    vhat += dsolλ[i]*fData.Pmax[i];
                end
                dsolλ[i] = 0;
            end
        end
        for i in bData.IDList
            dsolγ[i] = -getvalue(λw[i]);
            vhat -= dsolγ[i]*solw[i];
            if abs(dsolγ[i]) < 1e-4
                if dsolγ[i] < 0
                    vhat += dsolγ[i]*bData.cap[i];
                end
                dsolγ[i] = 0;
            end
            dsolμ[i] = -getvalue(λu[i]);
            vhat -= dsolμ[i]*solu[i];
            if abs(dsolμ[i]) < 1e-4
                if dsolμ[i] < 0
                    vhat += dsolμ[i]*bData.uCap[i];
                end
                dsolμ[i] = 0;
            end
        end

        cutTemp = cutData(statusDp,dsolλ,dsolγ,dsolμ,vhat);
        return cutTemp;
    else
        return dp;
    end

    # SOCP constraints
    # @constraint(mp, lineThermal1[k in fData.brList, t in td:T;Bparams[k,t] == 1], p[k,t]^2 + q[k,t]^2 <= fData.rateA[k]^2);
    # @constraint(mp, bThermal[i in bData.IDList, t in td:T], zp[i,t]^2 + zq[i,t]^2 <= u[i]^2);
    # @constraint(mp,gcAux1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux1[i,t] == f[i,t] +
    #     (fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
    # @constraint(mp,gcAux2[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux2[i,t] == f[i,t] +
    #     (-fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
    # @constraint(mp,gcAux3[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux3[i,t] == sqrt(fData.cp[i].params[1])*sp[i,t] +
    #     fData.cp[i].params[2]/(2*sqrt(fData.cp[i].params[1])));
    # @constraint(mp, genCost1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3],
    #     norm([tAux2[i,t],tAux3[i,t]]) <= tAux1[i,t]);
end

function constructBackwardM(td, T, Δt, trialPaths, matchedTrial, qpopt = false, hardened = [])
    # construct the math program given the state variables and current stage
    Ω = [ω for ω in keys(pDistr.ωDistrn)];
    paraSet = Iterators.product(Ω,matchedTrial);

    cutCurrentData = pmap(item -> dfBuild_D(td, pDistr.ωDict[item[1]], trialPaths[item[2]], pDistr.ωτ[item[1]], Δt, T, qpopt, true, hardened), paraSet);
    for j in procs()
        remotecall_fetch(cutUpdate,j,td,Ω,paraSet,cutCurrentData);
    end

    # for ω in Ω
    #     # solve the later stage problem
    #     cutCurrent = fBuild_D(td, ω, prevSol, τ, Δt, T, fData, bData, dData, pDistr, cutDict);
    #     if (td,ω) in keys(cutDict)
    #         push!(cutDict[td,ω],cutCurrent);
    #     else
    #         cutDict[td,ω] = [cutCurrent];
    #     end
    # end
    # return cutDict;
end

function exeBackward(T, Δt, trialPaths, qpopt = false, hardened = [])
    # execution of forward pass
    # input: trialPaths: the collection of trial points
    #        cutDict: previously generated cuts (preset in every core)
    # output: update the cutDict
    tpDict = Dict();
    for n in keys(trialPaths)
        tpDict[n] = [trialPaths[n][i][2] for i in 1:length(trialPaths[n])];
    end
    for t in T:-1:2
        matchedTrial = [];
        possiblePath = [];
        for n in keys(trialPaths)
            pathTemp = [(trialPaths[n][i][2],trialPaths[n][i][3],trialPaths[n][i][4]) for i in 2:length(trialPaths[n]) if trialPaths[n][i][2] < t];
            if t in tpDict[n]
                if !(pathTemp in possiblePath)
                    push!(matchedTrial,n);
                    push!(possiblePath,pathTemp);
                end
            end
        end
        if possiblePath != []
            constructBackwardM(t, T, Δt, trialPaths, matchedTrial, qpopt, hardened);
        end
        println("Time $(t) Passed");
    end
end

function exeBackward_last(T, Δt, trialPaths, qpopt = false, hardened = [])
    # execution of forward pass
    # input: trialPaths: the collection of trial points
    #        cutDict: previously generated cuts (preset in every core)
    # output: update the cutDict
    tpDict = Dict();
    for n in keys(trialPaths)
        tpDict[n] = [trialPaths[n][i][2] for i in 1:length(trialPaths[n])];
    end
    τMin = minimum([i for i in values(pDistr.ωτ)]);
    for t in T:-1:(T-τMin)
        matchedTrial = [];
        possiblePath = [];
        for n in keys(trialPaths)
            pathTemp = [(trialPaths[n][i][2],trialPaths[n][i][3],trialPaths[n][i][4]) for i in 2:length(trialPaths[n]) if trialPaths[n][i][2] < t];
            if t in tpDict[n]
                if !(pathTemp in possiblePath)
                    push!(matchedTrial,n);
                    push!(possiblePath,pathTemp);
                end
            end
        end
        if possiblePath != []
            constructBackwardM(t, T, Δt, trialPaths, matchedTrial, qpopt, hardened);
        end
        println("Time $(t) Passed");
    end
end

function exeBackwardAll(T, Δt, trialPaths, qpopt = false, hardened = [])
    # execution of forward pass
    # input: trialPaths: the collection of trial points
    #        cutDict: previously generated cuts (preset in every core)
    # output: update the cutDict
    tpDict = Dict();
    for n in keys(trialPaths)
        possibleTList = [t for t in 1:T];
        for i in 2:length(trialPaths[n])
            t = trialPaths[n][i][2];
            for tp in 1:trialPaths[n][i][4]
                delInd = findfirst(x -> x == t+tp, possibleTList);
                if delInd != nothing
                    deleteat!(possibleTList,delInd);
                end
            end
        end
        tpDict[n] = possibleTList;
    end
    for t in T:-1:2
        possiblePath = [];
        matchedTrial = [];
        for n in keys(trialPaths)
            if t in tpDict[n]
                pathTemp = [(trialPaths[n][i][2],trialPaths[n][i][3],trialPaths[n][i][4]) for i in 2:length(trialPaths[n]) if trialPaths[n][i][2] < t];
                if !(pathTemp in possiblePath)
                    push!(matchedTrial,n);
                    push!(possiblePath,pathTemp);
                end
            end
        end
        if possiblePath != []
            constructBackwardM(t, T, Δt, trialPaths, matchedTrial, qpopt, hardened);
        end
        println("Time $(t) Passed");
    end
end


text:
"""
`AskTellOptimizer` that utilizes randomization to generate the candidates.
"""
abstract type StochasticApproximationOptimizer <: AskTellOptimizer end

const SPSADefaultParameters = ParamsDict(
    :Alpha => 0.602,  # The optimal value is 1.0 but values down to 0.602 often can give faster convergence
    :Gamma => 0.101,  # The optimal value is 1/6 but values down to 0.101 often can give faster convergence
    :a     => 0.0017,
    :c     => 1.9, # Recommendation is value 1.9 but that assumes noisy function, otherwise should be low
    :A     => 10
)

mutable struct SimultaneousPerturbationSA2{E<:EmbeddingOperator} <: StochasticApproximationOptimizer
    embed::E # embed candidate into search space
    parameters::Parameters
    k::Int
    n::Int
    theta::Individual
    delta_ck::Individual

    function SimultaneousPerturbationSA2(problem::OptimizationProblem, embed::E, parameters::Parameters) where {E<:EmbeddingOperator}
        ss = search_space(problem)
        n = numdims(ss)
        new{E}(embed, chain(SPSADefaultParameters, parameters),
               0, n, rand_individual(ss), zeros(Float64, n))
    end
end

# by default use RandomBound embedder
SimultaneousPerturbationSA2(problem::OptimizationProblem, parameters::Parameters) =
    SimultaneousPerturbationSA2(problem, RandomBound(search_space(problem)), parameters)

name(spsa::SimultaneousPerturbationSA2) = "SPSA2 (Simultaneous Perturbation Stochastic Approximation, 1st order, 2 samples)"

sample_bernoulli_vector(n::Int) = [2.0*round(rand()) - 1.0 for _ in 1:n]

function ask(spsa::SimultaneousPerturbationSA2)
    delta = sample_bernoulli_vector(spsa.n)
    ck = spsa.parameters[:c]/(spsa.k + 1)^spsa.parameters[:Gamma]
    spsa.delta_ck = ck * delta

    theta_plus = spsa.theta + spsa.delta_ck
    theta_minus = spsa.theta - spsa.delta_ck

    [Candidate{Float64}(theta_plus, 1),
     Candidate{Float64}(theta_minus, 2)]
end

function tell!(spsa::SimultaneousPerturbationSA2, rankedCandidates::AbstractVector{<:Candidate})
    # Use index of rank to get right values for yplus and yminus, respectively.
    if rankedCandidates[1].index == 1
        yplus = rankedCandidates[1].fitness
        yminus = rankedCandidates[2].fitness
    else
        yplus = rankedCandidates[2].fitness
        yminus = rankedCandidates[1].fitness
    end

    # Estimate gradient.
    ghat = (yplus - yminus) ./ (2.0 * spsa.delta_ck)

    # Calc new estimate of theta based on estimate of gradient, ghat.
    ak = spsa.parameters[:a]/(spsa.k + 1 + spsa.parameters[:A])^spsa.parameters[:Alpha]
    theta_new = spsa.theta - ak * ghat
    apply!(spsa.embed, theta_new, spsa.theta)
    spsa.theta = theta_new
    spsa.k += 1

    return 0
end


text:
# Julia 1.7 doesn't infer this:
has_boxed_variables_f(f::F) where {F} =
    _any(ntuple(i -> fieldtype(typeof(f), i) <: Core.Box, Val(nfields(f)))...)

@generated has_boxed_variables_g(::F) where {F} = any(t -> t <: Core.Box, fieldtypes(F))

if VERSION < v"1.7-"
    const has_boxed_variables = has_boxed_variables_f
else
    const has_boxed_variables = has_boxed_variables_g
end

function verify_no_boxes(f::F, context) where {F}
    has_boxed_variables(f) && handle_boxed_variables(f, context())
    return
end

baremodule AssistantMode
using Base: @enum, UInt8
@enum Kind::UInt8 begin
    ignore
    warn
    warn_always
    error
end
end

const ASSISTANT_MODE = Ref(AssistantMode.warn)

struct SetAssistantResult
    old::AssistantMode.Kind
    new::AssistantMode.Kind
end

"""
    FLoops.assistant(mode::Symbol)
    FLoops.assistant(enable::Bool)

Set assistant mode; i.e., what to do when FLoops.jl finds a problematic usage
pattern.

Assistant modes:

* `:ignore`: do nothing
* `:warn`: print warning once
* `:warn_always`: print warning always
* `:error`: throw an error

`FLoops.assistant(false)` is equivalent to `FLoops.assistant(:ignore)` and
`FLoops.assistant(true)` is equivalent to `FLoops.assistant(:warn)`.
"""
assistant
function assistant(mode::Symbol)
    m = nothing
    if isdefined(AssistantMode, mode)
        m = getproperty(AssistantMode, mode)
    end
    if !(m isa AssistantMode.Kind)
        error(
            "invalid mode: ",
            mode,
            " (must be one of `:ignore`, `:warn`, `:warn_always`, and `:error`)",
        )
    end
    p = ASSISTANT_MODE[]
    ASSISTANT_MODE[] = m
    return SetAssistantResult(p, m)
end

assistant(enable::Bool) = enable ? assistant(:warn) : assistant(:ignore)

function Base.show(io::IO, ::MIME"text/plain", result::SetAssistantResult)
    printstyled(io, "FLoops.assistant"; color = :blue)
    print(io, ":")
    print(io, "\n  old mode: ", result.old)
    print(io, "\n  new mode: ", result.new)
end

function handle_boxed_variables(f, context)
    mode = ASSISTANT_MODE[]
    mode == AssistantMode.ignore && return

    err = HasBoxedVariableError(f)
    if mode == AssistantMode.error
        throw(err)
    elseif mode == AssistantMode.warn || mode == AssistantMode.warn_always
        ctx = context.ctx::MacroContext
        @warn(
            "Correctness and/or performance problem detected",
            error = err,
            _module = ctx.module_,
            _file = string(something(ctx.source.file, "none")),
            _line = ctx.source.line,
            _id = Symbol(context.id, mode == AssistantMode.warn ? :_once : :_always),
            maxlog = mode == AssistantMode.warn ? 1 : nothing,
        )
    else
        error("unknown mode: ", mode)
    end
    return
end

struct HasBoxedVariableError <: Exception
    f::Any
end

function Base.showerror(io::IO, err::HasBoxedVariableError)
    f = err.f
    varnames = [
        fieldname(typeof(f), i) for i in 1:nfields(f) if fieldtype(typeof(f), i) <: Core.Box
    ]

    print(io, "HasBoxedVariableError: Closure ")
    printstyled(io, nameof(f); color = :cyan)
    print(io, " (defined in ")
    printstyled(io, parentmodule(f); color = :cyan)
    print(io, ")")
    print(io, " has ", length(varnames), " boxed ")
    if length(varnames) == 1
        print(io, "variable:")
    else
        print(io, "variables:")
    end

    isfirst = true
    for v in varnames
        if isfirst
            isfirst = false
            print(io, ' ')
        else
            print(io, ", ")
        end
        printstyled(io, v; bold = true, color = :red)
    end

    println(io)
    print(io, "See: ")
    printstyled(
        io,
        "https://juliafolds.github.io/FLoops.jl/dev/howto/avoid-box/";
        color = :blue,
    )

    println(io)
    printstyled(io, "NOTE:", bold = true, color = :light_black)
    printstyled(
        io,
        " To disable this ",
        ASSISTANT_MODE[] == AssistantMode.error ? "error" : "warning",
        ", call `";
        color = :light_black,
    )
    printstyled(io, "FLoops.assistant(false)"; bold = true, color = :light_black)
    printstyled(io, "`."; color = :light_black)
end

function _make_closure_with_a_box()
    local a
    set(b) = (z = a; a = b; z)
    return set
end

function _make_closure_without_a_box(a = 1)
    add(b) = a + b
    return add
end

_box_detection_works() =
    has_boxed_variables(_make_closure_with_a_box()) &&
    !has_boxed_variables(_make_closure_without_a_box())

if try
    _box_detection_works()
catch err
    @error(
        "Error in `has_boxed_variables`. Disabling boxing detection.",
        exception = (err, catch_backtrace())
    )
    false
end
    const _verify_no_boxes = verify_no_boxes
else
    # Since `Core.Box` is internal, we can't rely on that it exists.
    const _verify_no_boxes = donothing
end


text:
"""
    ODESolvers

Ordinary differential equation solvers
"""
module ODESolvers

using GPUifyLoops
using StaticArrays
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

using ..SpaceMethods
using ..LinearSolvers
using ..MPIStateArrays: device, realview

export solve!, updatedt!, gettime

abstract type AbstractODESolver end
"""
    gettime(solver::AbstractODESolver)

Returns the current simulation time of the ODE solver `solver`
"""
gettime(solver::AbstractODESolver) = solver.t

"""
    getdt(solver::AbstractODESolver)

Returns the current simulation time step of the ODE solver `solver`
"""
getdt(solver::AbstractODESolver) = solver.dt

function dostep! end

"""
    updatedt!(solver::AbstractODESolver, dt)

Change the time step size to `dt` for the ODE solver `solver`.
"""
updatedt!(solver::AbstractODESolver, dt) =
  error("Variable time stepping not implemented for $(typeof(solver))")

"""
    updatetime!(solver::AbstractODESolver, time)

Change the current time to `time` for the ODE solver `solver`.
"""
updatetime!(solver::AbstractODESolver, time) =
  error("Variable time stepping not implemented for $(typeof(solver))")

isadjustable(solver::AbstractODESolver) = true

# {{{ run!
"""
    solve!(Q, solver::AbstractODESolver; timeend,
           stopaftertimeend=true, numberofsteps, callbacks)

Solves an ODE using the `solver` starting from a state `Q`. The state `Q` is
updated inplace. The final time `timeend` or `numberofsteps` must be specified.

A series of optional callback functions can be specified using the tuple
`callbacks`; see [`GenericCallbacks`](@ref).
"""
function solve!(Q, solver::AbstractODESolver, p=nothing; timeend::Real=Inf,
                adjustfinalstep=true, numberofsteps::Integer=0, callbacks=())

  @assert isfinite(timeend) || numberofsteps > 0
  if adjustfinalstep && !isadjustable(solver)
    error("$solver does not support time step adjustments. Can only be used with `adjustfinalstep=false`.")
  end
  t0 = gettime(solver)

  # Loop through an initialize callbacks (if they need it)
  foreach(callbacks) do cb
    try
      cb(true)
    catch
    end
  end

  step = 0
  time = t0
  while time < timeend
    step += 1

    time = dostep!(Q, solver, p, timeend, adjustfinalstep)

    # FIXME: Determine better way to handle postcallback behavior
    # Current behavior:
    #   retval = 1 exit after all callbacks
    #   retval = 2 exit immediately
    retval = 0
    for (i, cb) in enumerate(callbacks)
      # FIXME: Consider whether callbacks need anything, or if function closure
      #        can be used for everything
      thisretval = cb()
      thisretval = (thisretval == nothing) ? 0 : thisretval
      !(thisretval in (0, 1, 2)) &&
      error("callback #$(i) returned invalid value. It should return either:
            `nothing` (continue time stepping)
            `0`       (continue time stepping)
            `1`       (stop time stepping after all callbacks)
            `2`       (stop time stepping immediately)")
      retval = max(thisretval, retval)
      retval == 2 && return gettime(solver)
    end
    retval == 1 && return gettime(solver)

    # Figure out if we should stop
    if numberofsteps == step
      return gettime(solver)
    end
  end
  gettime(solver)
end
# }}}

include("LowStorageRungeKuttaMethod.jl")
include("StrongStabilityPreservingRungeKuttaMethod.jl")
include("AdditiveRungeKuttaMethod.jl")
include("MultirateInfinitesimalStepMethod.jl")
include("MultirateRungeKuttaMethod.jl")

end # module

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
"""Contains example usage of the features in the OpticSim.jl package."""
module Examples
using ..OpticSim
using ..OpticSim.Vis
# using ..OpticSim.GlassCat use this if you want to type SCHOTT.N_BK7 rather than OpticSim.GlassCat.SCHOTT.N_BK7
using StaticArrays
using DataFrames
using Images
using Unitful
using Plots
using LinearAlgebra

# Create a geometric hemisphere
function hemisphere()::CSGTree
    sph = Sphere(10.0)
    pln = Plane(0.0, 0.0, -1.0, 0.0, 0.0, 0.0)
    csgintersection(sph, pln)() #csg operations create a csggenerator which instantiates the csg tree after applying a rigid body transformation. This allows you to make as many instances of the object as you want with different transformations. We just want the CSGTree object rather than a generator.
end

# Create an optical hemisphere that has optical material properties so it will reflect and refract light. In the previous example the hemisphere object had optical properties of OpticSim.GlassCat.Air, which is the default optical interface, so it won't refract or reflect light.
function opticalhemisphere()::CSGOpticalSystem
    sph = Sphere(10.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))
    pln = Plane(0.0, 0.0, -1.0, 0.0, 0.0, 0.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))
    assy = LensAssembly{Float64}(csgintersection(sph, pln)())
    return CSGOpticalSystem(assy, Rectangle(1.0, 1.0, SVector{3,Float64}(0.0, 0.0, 1.0), SVector{3,Float64}(0.0, 0.0, -11.0)))
end

#! format: off

cooketriplet(::Type{T} = Float64, detpix::Int = 1000) where {T<:Real} = AxisymmetricOpticalSystem{T}(
    DataFrame(Surface = [:Object, 1, 2, 3, :Stop, 5, 6, :Image],
              Radius = [Inf, 26.777, 66.604, -35.571, 35.571, 35.571, -26.777, Inf],
              OptimizeRadius = [false,true,true,true,true,true,true,false],
              Thickness = [Inf, 4.0, 2.0, 4.0, 2.0, 4.0, 44.748, missing],
              OptimizeThickness = [false,true,true,true,true,true,true,false],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SK16, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SF2, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SK16, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf, 8.580, 7.513, 7.054, 6.033, 7.003, 7.506, 15.0]), detpix, detpix)
export cooketriplet

#no longer works
cooketripletlensonly(::Type{T} = Float64) where {T<:Real} = AxisymmetricLens{T}(
    DataFrame(Surface = [:Object, 1, 2, 3, :Stop, 5, 6, :Image],
              Radius = [Inf, 26.777, 66.604, -35.571, 35.571, 35.571, -26.777, Inf],
              OptimizeRadius = [false,true,true,true,true,true,true,false],
              Thickness = [Inf, 4.0, 2.0, 4.0, 2.0, 4.0, 44.748, missing],
              OptimizeThickness = [false,true,true,true,true,true,true,false],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SK16, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SF2, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SK16, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf, 8.580, 7.513, 7.054, 6.033, 7.003, 7.506, 15.0]))
export cooketripletlensonly

cooketripletfirstelement(::Type{T} = Float64) where {T<:Real} = AxisymmetricOpticalSystem(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [Inf, -35.571, 35.571, Inf],
              Thickness = [Inf, 4.0, 44.748, missing],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SK16, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf, 7.054, 6.033, 15.0]))

convexplano(::Type{T} = Float64) where {T<:Real} = AxisymmetricOpticalSystem{T}(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [Inf, 60.0, Inf, Inf],
              Thickness = [Inf, 10.0, 57.8, missing],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf, 9.0, 9.0, 15.0]))

doubleconvex(frontradius::T,rearradius::T) where{T<:Real} =
AxisymmetricOpticalSystem{T}(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [T(Inf64), frontradius, rearradius, T(Inf64)],
              OptimizeRadius = [false,true,true,false],
              Thickness = [T(Inf64), T(10.0), T(57.8), missing],
              OptimizeThickness = [false,false,false,false],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [T(Inf64), T(9.0), T(9.0), T(15.0)]))

doubleconvexconic(::Type{T} = Float64) where {T<:Real} =
              AxisymmetricOpticalSystem{T}(
                  DataFrame(Surface = [:Object, 1, 2, :Image],
                            Radius = [Inf64, 60, -60, Inf64],
                            OptimizeRadius = [false,true,true,false],
                            Thickness = [Inf64, 10.0, 57.8, missing],
                            OptimizeThickness = [false,false,false,false],
                            Conic = [missing, 0.01, 0.01, missing],
                            OptimizeConic = [false, true, true, false],
                            Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
                            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]))

doubleconvexlensonly(frontradius::T,rearradius::T) where{T<:Real} =
AxisymmetricLens{T}(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [T(Inf64), frontradius, rearradius, T(Inf64)],
              OptimizeRadius = [false,true,true,false],
              Thickness = [T(Inf64), T(10.0), T(57.8), missing],
              OptimizeThickness = [false,false,false,false],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [T(Inf64), T(9.0), T(9.0), T(15.0)]))
export doubleconvexlensonly

doubleconvexprescription() = DataFrame(
    Surface = [:Object, 1, 2, :Image],
    Radius = [Inf64, 60, -60, Inf64],
    OptimizeRadius = [false,true,true,false],
    Thickness = [Inf64, 10.0, 57.8, missing],
    OptimizeThickness = [false,true,true,false],
    Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
    SemiDiameter = [Inf64, 9.0, 9.0, 15.0])

doubleconvex(::Type{T} = Float64; temperature::Unitful.Temperature = OpticSim.GlassCat.TEMP_REF_UNITFUL, pressure::T = T(OpticSim.GlassCat.PRESSURE_REF)) where {T<:Real} = AxisymmetricOpticalSystem{T}(doubleconvexprescription(),temperature = temperature, pressure = pressure)

doubleconcave(::Type{T} = Float64) where {T<:Real} = AxisymmetricOpticalSystem{T}(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [Inf64, -41.0, 41.0, Inf64],
              Thickness = [Inf64, 10.0, 57.8, missing],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf64, 9.0, 9.0, 15.0]))

planoconcaverefl(::Type{T} = Float64) where {T<:Real} = AxisymmetricOpticalSystem{T}(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [Inf64, Inf64, -41.0, Inf64],
              Thickness = [Inf64, 10.0, -57.8, missing],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf64, 9.0, 9.0, 25.0],
              Reflectance = [missing, missing, 1.0, missing]))

concaveplano(::Type{T} = Float64) where {T<:Real} = AxisymmetricOpticalSystem{T}(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [Inf64, -41.0, Inf64, Inf64],
              Thickness = [Inf64, 10.0, 57.8, missing],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf64, 9.0, 9.0, 15.0]))

planoplano(::Type{T} = Float64) where {T<:Real} = AxisymmetricOpticalSystem{T}(
    DataFrame(Surface = [:Object, 1, 2, :Image],
              Radius = [Inf64, Inf64, Inf64, Inf64],
              Thickness = [Inf64, 10.0, 57.8, missing],
              Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
              SemiDiameter = [Inf64, 9.0, 9.0, 15.0]))

#! format: on

end #module SphericalLenses

function autodrawrays(lens::AxisymmetricOpticalSystem = cooketriplet(), angle = 10; kwargs...)
    f1 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourcenum = 1)
    Vis.drawtracerays(lens, raygenerator = f1, test = true, trackallrays = true, colorbysourcenum = true; kwargs...)
    f2 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourceangle = angle / 180 * π, sourcenum = 2)
    Vis.drawtracerays!(lens, raygenerator = f2, test = true, trackallrays = true, colorbysourcenum = true; kwargs...)
end

function autospotdiag(lens::AxisymmetricOpticalSystem = cooketriplet(); kwargs...)
    f1 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourcenum = 1)
    f2 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourceangle = 5 / 180 * π, sourcenum = 2)
    f3 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourceangle = 10 / 180 * π, sourcenum = 3)
    Vis.spotdiaggrid(lens, [f1, f2, f3]; kwargs...)
end

# Display the spot diagram of a simple cooketriplet lens
function hexapolarspotdiagramexample(lens = cooketriplet(), numrings::Int = 5, angle = 0.0)
    Vis.spotdiag(lens, samples = numrings, sourceangle = angle)
end

function cartesiangridspotdiagramexample(lens = cooketriplet(), numsamples::Int = 5, angle = 0.0)
    Vis.spotdiag(lens, hexapolar = false, samples = numsamples, sourceangle = angle)
end

function SchmidtCassegrainTelescope()
    # glass entrance lens on telescope
    topsurf = Plane(SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air), vishalfsizeu = 12.00075, vishalfsizev = 12.00075)
    botsurf = AcceleratedParametricSurface(ZernikeSurface(12.00075, radius = -1.14659768e+4, aspherics = [(4, 3.68090959e-7), (6, 2.73643352e-11), (8, 3.20036892e-14)]), 17, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))
    coverlens = csgintersection(leaf(Cylinder(12.00075, 1.4)), csgintersection(leaf(topsurf), leaf(botsurf, RigidBodyTransform(OpticSim.rotmatd(0, 180, 0), SVector(0.0, 0.0, -0.65)))))
    # big mirror with a hole in it
    bigmirror = ConicLens(OpticSim.GlassCat.SCHOTT.N_BK7, -72.65, -95.2773500000134, 0.077235, Inf, 0.0, 0.2, 12.18263, frontsurfacereflectance = 1.0)
    bigmirror = csgdifference(bigmirror, leaf(Cylinder(4.0, 0.3, interface = opaqueinterface()), translation(0.0, 0.0, -72.75)))
    # small mirror supported on a spider
    smallmirror = SphericalLens(OpticSim.GlassCat.SCHOTT.N_BK7, -40.65, Inf, -49.6845, 1.13365, 4.3223859, backsurfacereflectance = 1.0)
    obscuration1 = Circle(4.5, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -40.649), interface = opaqueinterface())
    obscurations2 = Spider(3, 0.5, 12.0, SVector(0.0, 0.0, -40.65))
    # put it together with the detector
    la = LensAssembly(coverlens(), bigmirror(), smallmirror(), obscuration1, obscurations2...)
    det = Circle(3.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -92.4542988), interface = opaqueinterface())

    return CSGOpticalSystem(la, det)
end

drawSchmidt(; kwargs...) = Vis.drawtracerays(SchmidtCassegrainTelescope(), raygenerator = UniformOpticalSource(CollimatedSource(GridRectOriginPoints(5, 5, 10.0, 10.0, position = SVector(0.0, 0.0, 20.0))), 0.55), trackallrays = true, colorbynhits = true, test = true, numdivisions = 100; kwargs...)

function prism_refraction()
    # build the triangular prism
    int = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_SF14, OpticSim.GlassCat.Air)
    s = 2.0
    prism = csgintersection(leaf(Plane(SVector(0.0, -1.0, 0.0), SVector(0.0, -s, 0.0), interface = int, vishalfsizeu = 2 * s, vishalfsizev = 2 * s)), csgintersection(Plane(SVector(0.0, sind(30), cosd(30)), SVector(0.0, s * sind(30), s * cosd(30)), interface = int, vishalfsizeu = 2 * s, vishalfsizev = 2 * s), Plane(SVector(0.0, sind(30), -cosd(30)), SVector(0.0, s * sind(30), -s * cosd(30)), interface = int, vishalfsizeu = 2 * s, vishalfsizev = 2 * s)))
    sys = CSGOpticalSystem(LensAssembly(prism()), Rectangle(15.0, 15.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -20.0), interface = opaqueinterface()))
    # create some 'white' light
    rays = Vector{OpticalRay{Float64,3}}(undef, 0)
    for i in 0:7
        λ = ((i / 7) * 200 + 450) / 1000
        r = OpticalRay(SVector(0.0, -3.0, 10.0), SVector(0.0, 0.5, -1.0), 1.0, λ)
        push!(rays, r)
    end
    raygen = RayListSource(rays)
    # draw the result
    Vis.drawtracerays(sys, raygenerator = raygen, test = true, trackallrays = true)
end

function zoom_lens(pos = 1)
    if pos == 0
        stop = 2.89
        zoom = 9.48
        dist = 4.46970613
    elseif pos == 1
        stop = 3.99
        zoom = 4.48
        dist = 21.21
    else
        stop = 4.90
        zoom = 2.00
        dist = 43.81
    end
    #! format: off
    return AxisymmetricOpticalSystem{Float64}(
        DataFrame(Surface = [:Object, :Stop, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, :Image],
                  Radius = [Inf64, Inf64, -1.6202203499676E+01, -4.8875855327468E+01, 1.5666614444619E+01, -4.2955326460481E+01, 1.0869565217391E+02, 2.3623907394283E+01, -1.6059097478722E+01, -4.2553191489362E+02, -3.5435861091425E+01, -1.4146272457208E+01, -2.5125628140704E+02, -2.2502250225023E+01, -1.0583130489999E+01, -4.4444444444444E+01, Inf64],
                  Aspherics = [missing, missing, missing, missing, missing, [(4, 1.03860000000E-04), (6, 1.42090000000E-07), (8, -8.84950000000E-09), (10, 1.24770000000E-10), (12, -1.03670000000E-12), (14, 3.65560000000E-15)], missing, missing, [(4, 4.27210000000E-05), (6, 1.24840000000E-07), (8, 9.70790000000E-09), (10, -1.84440000000E-10), (12, 1.86440000000E-12), (14, -7.79750000000E-15)], [(4, 1.13390000000E-04), (6, 4.81650000000E-07), (8, 1.87780000000E-08), (10, -5.75710000000E-10), (12, 8.99940000000E-12), (14, -4.67680000000E-14)], missing, missing, missing, missing, missing, missing, missing],
                  Thickness = [Inf64, 0.0, 5.18, 0.10, 4.40, 0.16, 1.0, 4.96, zoom, 4.04, 1.35, 1.0, 2.80, 3.0, 1.22, dist, missing],
                  Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_LAH66, OpticSim.GlassCat.Air, OpticSim.GlassCat.NIKON.LLF6, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_TIH6, OpticSim.GlassCat.OHARA.S_FSL5, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_FSL5, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_LAL8, OpticSim.GlassCat.SCHOTT.S_FL4, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_LAH66, OpticSim.GlassCat.Air, missing],
                  SemiDiameter = [Inf64, stop, 3.85433218451, 3.85433218451, 4.36304692871, 4.36304692871, 4.72505505439, 4.72505505439, 4.72505505439, 4.45240784026, 4.45240784026, 4.50974054117, 4.50974054117, 4.50974054117, 4.76271114409, 4.76271114409, 15.0]))
    #! format: on

end

function fresnel(convex = true; kwargs...)
    lens = FresnelLens(OpticSim.GlassCat.SCHOTT.N_BK7, 0.0, convex ? 15.0 : -15.0, 1.0, 8.0, 0.8, conic = 0.1)
    sys = CSGOpticalSystem(LensAssembly(lens()), Rectangle(15.0, 15.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -25.0), interface = opaqueinterface()))
    Vis.drawtracerays(sys; test = true, trackallrays = true, numdivisions = 30, kwargs...)
end

function grating(; period = 1.0, θ = 0.0, λ = 0.55, kwargs...)
    int = ThinGratingInterface(SVector(0.0, 1.0, 0.0), period, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, minorder = -2, maxorder = 2, reflectance = [0.0, 0.0, 0.1, 0.0, 0.0], transmission = [0.05, 0.1, 0.4, 0.1, 0.05])
    grating = ThinGratingSurface(Rectangle(5.0, 5.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0)), int)
    back = Rectangle(30.0, 30.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 25.0))
    sys = CSGOpticalSystem(LensAssembly(grating, back), Rectangle(30.0, 30.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -25.0), interface = opaqueinterface()))
    Vis.drawtracerays(sys; raygenerator = UniformOpticalSource(CollimatedSource(OriginPoint{Float64}(100, position = SVector(0.0, 0.0, 10.0), direction = SVector(0.0, sind(θ), -cosd(θ)))), λ), trackallrays = true, rayfilter = nothing, kwargs...)
end

function reflgrating(; period = 1.0, θ = 0.0, λ = 0.55, kwargs...)
    int = ThinGratingInterface(SVector(0.0, 1.0, 0.0), period, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, minorder = -2, maxorder = 2, transmission = [0.0, 0.0, 0.1, 0.0, 0.0], reflectance = [0.05, 0.1, 0.4, 0.1, 0.05])
    grating = ThinGratingSurface(Rectangle(5.0, 5.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0)), int)
    back = Rectangle(30.0, 30.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -25.0))
    sys = CSGOpticalSystem(LensAssembly(grating, back), Rectangle(30.0, 30.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 25.0), interface = opaqueinterface()))
    Vis.drawtracerays(sys; raygenerator = UniformOpticalSource(CollimatedSource(OriginPoint{Float64}(100, position = SVector(0.0, 0.0, 10.0), direction = SVector(0.0, sind(θ), -cosd(θ)))), λ), trackallrays = true, rayfilter = nothing, kwargs...)
end

function HOE(refl = false, firstorderonly = false; kwargs...)
    rect = Rectangle(5.0, 5.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0))
    if refl
        int = HologramInterface(SVector(0.0, -10.0, 20.0), ConvergingBeam, SVector(0.0, 0.0, -200), ConvergingBeam, 0.55, 9.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, !firstorderonly)
    else
        int = HologramInterface(SVector(0.0, -10.0, -20.0), ConvergingBeam, SVector(0.0, 0.0, -200), ConvergingBeam, 0.55, 5.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, !firstorderonly)
    end
    obj = HologramSurface(rect, int)
    back = Rectangle(50.0, 50.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 25.0))
    sys = CSGOpticalSystem(LensAssembly(obj, back), Rectangle(50.0, 50.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -25.0), interface = opaqueinterface()))
    Vis.drawtracerays(sys; raygenerator = UniformOpticalSource(GridSource(OriginPoint{Float64}(10, position = SVector(0.0, 0.0, 10.0), direction = SVector(0.0, 0.0, -1.0)), 1, 15, 0.0, π / 6), 0.55), trackallrays = true, rayfilter = nothing, kwargs...)
end

function HOEfocus(; kwargs...)
    rect = Rectangle(5.0, 5.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0))
    int = HologramInterface(SVector(0.0, -3.0, -20.0), ConvergingBeam, SVector(0.0, 0.0, -1.0), CollimatedBeam, 0.55, 9.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, false)
    obj = HologramSurface(rect, int)
    sys = CSGOpticalSystem(LensAssembly(obj), Rectangle(10.0, 10.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -25.0), interface = opaqueinterface()))
    Vis.drawtracerays(sys; raygenerator = UniformOpticalSource(CollimatedSource(GridRectOriginPoints(5, 5, 3.0, 3.0, position = SVector(0.0, 0.0, 10.0), direction = SVector(0.0, 0.0, -1.0))), 0.55), trackallrays = true, rayfilter = nothing, test = true, kwargs...)
end

function HOEcollimate(; kwargs...)
    rect = Rectangle(5.0, 5.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0))
    int = HologramInterface(SVector(0.1, -0.05, -1.0), CollimatedBeam, SVector(0.0, 0.0, 10), DivergingBeam, 0.55, 9.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, false)
    obj = HologramSurface(rect, int)
    sys = CSGOpticalSystem(LensAssembly(obj), Rectangle(10.0, 10.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -25.0), interface = opaqueinterface()))
    Vis.drawtracerays(sys; raygenerator = UniformOpticalSource(GridSource(OriginPoint{Float64}(1, position = SVector(0.0, 0.0, 10.0), direction = SVector(0.0, 0.0, -1.0)), 5, 5, π / 4, π / 4), 0.55), trackallrays = true, rayfilter = nothing, test = true, kwargs...)
end

function eyetrackHOE(nrays = 5000, det = false, showhead = true, zeroorder = false; kwargs...)
    # TODO update for new specs from Chris
    hoehalfwidth = 50.0 #25.0
    hoehalfheight = 50.0 #22.5
    hoecenter = SVector(-8.0 - 25.0, 0.0, -10.0 - 25.0)
    rect = Rectangle(hoehalfheight, hoehalfwidth, SVector(0.0, 1.0, 0.0), hoecenter)
    er = 15.0
    cornea_rad = 7.85
    corneavertex = SVector(0.0, er, 0.0)
    sourceloc = SVector(-33.0, er, 0.0)
    camloc = SVector(20.0, 3.0, -11.0)
    camdir = corneavertex - camloc
    camdir_norm = normalize(camdir)

    interfaces = []

    # offset = SVector(-5.0, 10.0, -10.0)
    # for θ in 0:(π / 6):(2π)
    #     ledloc = SVector(20 * cos(θ) + offset[1], 0 + offset[2], 15 * sin(θ) + offset[3])
    #     int = HologramInterface(ledloc, ConvergingBeam, sourceloc, DivergingBeam, 0.78, 100.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, zeroorder)
    #     push!(interfaces, int)
    # end

    dirs = [SVector(0.7713, 0.6350, -0.0437), SVector(0.5667, 0.8111, -0.1445), SVector(0.3400, 0.9349, -0.1017), SVector(0.1492, 0.9878, 0.0445), SVector(0.0249, 0.9686, 0.2474), SVector(-0.0184, 0.8855, 0.4643), SVector(0.0254, 0.7537, 0.6567), SVector(0.1548, 0.5964, 0.7876), SVector(0.3570, 0.4462, 0.8207), SVector(0.5959, 0.3470, 0.7242), SVector(0.7976, 0.3449, 0.4948), SVector(0.8680, 0.4555, 0.1978)]

    for d in dirs
        int = HologramInterface(normalize(d), CollimatedBeam, sourceloc, DivergingBeam, 0.78, 100.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, zeroorder)
        # int = HologramInterface(corneavertex - 10 * d, ConvergingBeam, sourceloc, DivergingBeam, 0.78, 100.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, zeroorder)
        push!(interfaces, int)
    end

    mint = MultiHologramInterface(interfaces...)
    obj = MultiHologramSurface(rect, mint)
    cornea = leaf(Sphere(cornea_rad, interface = FresnelInterface{Float64}(OpticSim.GlassCat.EYE.CORNEA, OpticSim.GlassCat.Air, reflectance = 1.0, transmission = 0.0)), translation(0.0, er + cornea_rad, 0.0))()

    # cam settings
    fnum = 2.0
    fov = 80
    sensorrad = 1.0
    barrellength = sensorrad / tand(fov / 2)
    aprad = barrellength / fnum / 2
    camrad = max(sensorrad, aprad)

    camap = Annulus(aprad, camrad, camdir_norm, camloc)
    distfromcamtoeye = norm(camdir)
    focallength = 1 / (1 / distfromcamtoeye + 1 / barrellength)
    camlens = ParaxialLensEllipse(focallength, aprad, aprad, -camdir_norm, camloc)
    barrelloc = camloc - barrellength / 2 * camdir_norm
    barreltop = Plane(camdir_norm, camloc)
    barrelbot = Plane(-camdir_norm, camloc - 3 * barrellength * camdir_norm)
    barrelrot = OpticSim.rotmatbetween(SVector(0.0, 0.0, 1.0), camdir_norm)
    cambarrel = csgintersection(barrelbot, csgintersection(barreltop, leaf(Cylinder(camrad, barrellength, interface = opaqueinterface(Float64)), RigidBodyTransform(barrelrot, barrelloc))))()
    camdet = Circle(sensorrad, camdir_norm, camloc - barrellength * camdir_norm, interface = opaqueinterface(Float64))

    # sourceleft = hoecenter[1] + hoehalfwidth - sourceloc[1]
    # sourceright = hoecenter[1] - hoehalfwidth - sourceloc[1]
    # sourceleftθ = atan(sourceleft, sourceloc[2])
    # sourcerightθ = atan(sourceright, sourceloc[2])
    # midθ = (sourceleftθ + sourcerightθ) / 2
    # sourcedir = normalize(SVector(er * tan(midθ), -er, 0.0))
    # sourceextentθ = abs(midθ - sourcerightθ)
    # source = CosineOpticalSource(RandomSource(OriginPoint{Float64}(1, position = sourceloc, direction = sourcedir), nrays, sourceextentθ), 1.0, 0.78)

    rays = Vector{OpticalRay{Float64,3}}(undef, nrays)
    @simd for i in 1:nrays
        p = point(rect, rand() * 2 - 1, rand() * 2 - 1)
        rays[i] = OpticalRay(sourceloc, p - sourceloc, 1.0, 0.78)
    end
    source = RayListSource(rays)

    sys = CSGOpticalSystem(LensAssembly(obj, cornea, camlens, cambarrel, camap), camdet, 800, 800)
    if det
        Vis.show(OpticSim.traceMT(sys, source))
    else
        Vis.drawtracerays(sys; raygenerator = source, trackallrays = true, kwargs...)
        # for θ in 0:(π / 6):(2π)
        #     ledloc = SVector(20 * cos(θ) + offset[1], 0 + offset[2], 15 * sin(θ) + offset[3])
        #     Vis.draw!(leaf(Sphere(1.0), translation(ledloc...)), color = :red)
        # end
        for d in dirs
            # Vis.draw!(leaf(Sphere(1.0), translation((corneavertex - 10 * d)...)), color = :red)
            Vis.draw!((corneavertex - 50 * d, corneavertex), color = :red)
        end
        if showhead
            Vis.draw!(joinpath(@__DIR__, "../../OBJ/glasses.obj"), scale = 100.0, transform = RigidBodyTransform(OpticSim.rotmatd(90, 0, 0), [27.0, 45.0, -8.0]), color = :black)
            Vis.draw!(joinpath(@__DIR__, "../../OBJ/femalehead.obj"), scale = 13.0, transform = RigidBodyTransform(OpticSim.rotmatd(0, 0, 180), [27.0, 105.0, -148.0]), color = :white)
        end
        Vis.display()
    end
end

function multiHOE(; kwargs...)
    rect = Rectangle(5.0, 5.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0))
    int1 = HologramInterface(SVector(-5.0, 0.0, -20.0), ConvergingBeam, SVector(0.0, -1.0, -1.0), CollimatedBeam, 0.55, 100.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, false)
    int2 = HologramInterface(SVector(5.0, 0.0, -20.0), ConvergingBeam, SVector(0.0, 1.0, -1.0), CollimatedBeam, 0.55, 100.0, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, false)
    mint = MultiHologramInterface(int1, int2)
    obj = MultiHologramSurface(rect, mint)
    sys = CSGOpticalSystem(LensAssembly(obj), Rectangle(10.0, 10.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -20.0), interface = opaqueinterface()))
    s1 = UniformOpticalSource(CollimatedSource(RandomRectOriginPoints(500, 3.0, 3.0, position = SVector(0.0, 3.0, 3.0), direction = SVector(0.0, -1.0, -1.0))), 0.55, sourcenum = 1)
    s2 = UniformOpticalSource(CollimatedSource(RandomRectOriginPoints(500, 3.0, 3.0, position = SVector(0.0, -3.0, 3.0), direction = SVector(0.0, 1.0, -1.0))), 0.55, sourcenum = 2)
    s3 = UniformOpticalSource(CollimatedSource(RandomRectOriginPoints(500, 3.0, 3.0, position = SVector(0.0, 0.0, 3.0), direction = SVector(0.0, 0.0, -1.0))), 0.55, sourcenum = 3)
    Vis.drawtracerays(sys; raygenerator = s1, trackallrays = true, colorbysourcenum = true, rayfilter = nothing, kwargs...)
    Vis.drawtracerays!(sys; raygenerator = s2, trackallrays = true, colorbysourcenum = true, rayfilter = nothing, drawgen = true, kwargs...)
    Vis.drawtracerays!(sys; raygenerator = s3, trackallrays = true, colorbysourcenum = true, rayfilter = nothing, drawgen = true, kwargs...)
end
export Examples


text:
# This file includes an example script to plot the basis functions of PSK modulator.

using DigitalCommunications
using Plots

# Simulation parameters
k = 2
M = 2^k
nsymbols = Int(1e6)
nbits = k * nsymbols
ebno = collect(0 : 10)          # Snr per bit
esno = ebno .+ 10 * log10(k)    # Snr ber symbol

# Communcation system components
gen = BitGenerator(nbits)
modulator = Modulator(FSK(M))

# Plot the basis of the modulator
t = 0 : modulator.tsample : modulator.pulse.duration
plt = plot(layout=(constelsize(modulator.scheme),1))
for (i, base) in enumerate(basis(modulator))
    plot!(t, base.(t), subplot=i)
end
display(plt)


text:
#MDDatasets: Tools to store/manipulate multi-dimensional data
#-------------------------------------------------------------------------------
__precompile__(true)
#=
TAGS:
	#WANTCONST, HIDEWARN_0.7
=#

module MDDatasets

import Base: values
using Statistics

#==Suggested scalar data types
	-Use concrete types of the largest size natively supported by processor
   -Eventually should move to 128-bit values, etc.==#
const DataFloat = Float64
const DataInt = Int64
const DataComplex = Complex{Float64}
#==NOTES:
Don't use Int.  If this is like c, "Int" represents the "native" integer
size/type.  That might not be the same as the largest integer that can be
handled with reasonable efficiency

No reason to have an alias for Bool.  Probably best to keep the default
representation.==#

#Type used to dispatch on a symbol & minimize namespace pollution:
#-------------------------------------------------------------------------------
struct DS{Symbol}; end; #Dispatchable symbol
DS(v::Symbol) = DS{v}()

#-------------------------------------------------------------------------------
include("functions.jl")
include("math.jl")
include("arrays.jl")
include("base.jl")
include("datahr.jl")
include("datars.jl")
include("datasetop.jl")
include("meas_cross.jl")
include("meas_cksignals.jl")
include("meas_binsignals.jl")
include("broadcast.jl")
include("broadcast_datahr.jl")
include("broadcast_datars.jl")
include("datasetop_reg.jl")
include("measinterface.jl")
include("show.jl")
include("doc.jl")

#==TODO: Watch out for value() being exported by multiple modules...
Maybe it can be defined in "Units"
==#

#Data types:
#-------------------------------------------------------------------------------
export DataFloat, DataInt, DataComplex
export DataMD #Prefered abstraction for high-level functions
export PSweep #Parameter sweep
export Index #A way to identify parameters as array indicies
export DataFloat, DataInt, DataComplex, DataF1 #Leaf data types
export DataHR, DataRS

#Support functions:
#-------------------------------------------------------------------------------
export ensure

#Scalar operations:
#-------------------------------------------------------------------------------
export sanitize #Clamp down infinite values & substitute NaNs

#Accessor functions:
#-------------------------------------------------------------------------------
export value #High-collision WARNING: other modules probably want to export "value".
export sweep #Access values from a particular parameter sweep.
export sweeps #Get the list of parameter sweeps in DataHR.
#Base.ndims #Return number of dimensions (parametric sweeps).
export dimension #Returns the dimension corresponding to the given parameter name.
export subscripts #Provides subscripts iterator to access each element in DataHR.
export getsubarray
export coordinates #Get parameter sweep coordinates corresponding to given subscripts.
export paramlist #Return a list of parameter values being swept.
export parameter #Get parameter values for a particular sweep.

#Basic dataset operations:
#-------------------------------------------------------------------------------
@doc "Gets a dataset with all the x-values as the y-values" xval
export xval #Gets a dataset with all the x-values as the y-values
#value(ds, x=value) (already exported): get y-value @ given x.
export clip #clips a dataset within an xrange
export sample #Samples a dataset
export delta, xshift, xscale
export yvsx
export deriv, iinteg #derivative, definite integral, indefinite integral
export integ, xmin, xmax

#Cross operations:
#-------------------------------------------------------------------------------
#==API
NOTE: using "{Event}" for type stability (I hope)

Returns x-values of d1 (up-to nmax) when d1 crosses 0 (nmax=0: get all crossings):
xcross(d1, nmax; xstart, allow::CrossType)
xcross(Event, d1, nmax; xstart, allow::CrossType)

Returns y-values of d2 (up-to nmax) when when d1 crosses d2 (nmax=0: get all crossings):
ycross(d1, d2, nmax; xstart, allow::CrossType)
ycross(Event, d1, d2, nmax; xstart, allow::CrossType)

Returns scalar x-value of d1 on n-th zero-crossing:
xcross1(d1; n, xstart, allow::CrossType)
xcross1(Event, d1; n, xstart, allow::CrossType)

Returns scalar y-value of d1 on n-th crossing of d1 & d2:
ycross1(d1, d2; n, xstart, allow::CrossType)
ycross1(Event, d1, d2; n, xstart, allow::CrossType)

Resultant x-values are x-values of d1 @ crossings:
(x/y)cross

Resultant x-values corresponds to number of current crossing event (1, 2, ..., n):
(x/y)cross{Event}
==#
export Event #Identifies result as having event count in x-axis
export CrossType #To filter out unwanted crossings
export xcross, xcross1 #Measure x @ crossing events
export ycross, ycross1 #Measure y @ crossing events
export measdelay #Measure delay between crossing events of two signals
export measperiod, measfreq, measduty, measckstats
export measck2q, measskew
export measrise, measfall

export meas #Use meas(:MEASTYPE, ...) to minimize namespace pollution.
#==Sample usage:
TODO: Deprecate direct call of x/ycross/1 & measdelay
ON CONDITION: Only if this interface does not mess up type stability.
meas(:xcross, Event, d1)
==#

#==Initialization
===============================================================================#
function __init__()
	show(:HACK_SHOWTOUNFREEZE) #Not sure why show unfreezes dependent modules (those using this one)
	println() #newline
	return
end

end #MDDatasets

#Last line