text:
module Lora

using Distributions
using Formatting
using StatsBase

import Base:
  ==,
  Dict,
  close,
  convert,
  copy!,
  eltype,
  flush,
  getindex,
  isequal,
  keys,
  mark,
  mean,
  mean!,
  open,
  rand,
  read!,
  read,
  reset,
  run,
  show,
  write,
  writemime

import Distributions:
  @check_args,
  @distr_support,
  failprob,
  logpdf,
  params,
  pdf,
  succprob

import ForwardDiff

import ReverseDiffSource

export
  ### Types
  AM,
  AMState,
  ARS,
  ARSState,
  AcceptanceRateMCTuner,
  BasicContMuvParameter,
  BasicContMuvParameterNState,
  BasicContMuvParameterState,
  BasicContParamIOStream,
  BasicContUnvParameter,
  BasicContUnvParameterNState,
  BasicContUnvParameterState,
  BasicDiscMuvParameter,
  BasicDiscMuvParameterNState,
  BasicDiscMuvParameterState,
  BasicDiscParamIOStream,
  BasicDiscUnvParameter,
  BasicDiscUnvParameterNState,
  BasicDiscUnvParameterState,
  BasicGibbsJob,
  BasicMCJob,
  BasicMCRange,
  BasicMCTune,
  BasicMavVariableNState,
  BasicMavVariableState,
  BasicMuvVariableNState,
  BasicMuvVariableState,
  BasicUnvVariableNState,
  BasicUnvVariableState,
  BasicVariableIOStream,
  Binary,
  Constant,
  ContMuvMarkovChain,
  ContUnvMarkovChain,
  ContinuousParameter,
  ContinuousParameterNState,
  ContinuousParameterState,
  Data,
  Dependence,
  DependenceVector,
  Deterministic,
  DiscMuvMarkovChain,
  DiscUnvMarkovChain,
  DiscreteParameter,
  DiscreteParameterNState,
  DiscreteParameterState,
  DualAveragingMCTuner,
  GenericModel,
  GibbsJob,
  HMC,
  HMCSampler,
  HMCSamplerState,
  HMCState,
  Hyperparameter,
  IntegerVector,
  LMCSampler,
  LMCSamplerState,
  MALA,
  MALAState,
  MCJob,
  MCRange,
  MCSampler,
  MCSamplerState,
  MCTuner,
  MCTunerState,
  MH,
  MHSampler,
  MHSamplerState,
  MHState,
  MarkovChain,
  MatrixvariateParameter,
  MatrixvariateParameterNState,
  MatrixvariateParameterState,
  MultivariateParameter,
  MultivariateParameterNState,
  MultivariateParameterState,
  MuvAMState,
  MuvAMWGState,
  MuvHMCState,
  MuvMALAState,
  MuvSMMALAState,
  Parameter,
  ParameterIOStream,
  ParameterNState,
  ParameterState,
  ParameterStateVector,
  ParameterVector,
  RAM,
  RAMState,
  Random,
  RealMatrix,
  RealVector,
  RobertsRosenthalMCTune,
  RobertsRosenthalMCTuner,
  SMMALA,
  SMMALAState,
  Sampleability,
  Transformation,
  UnivariateParameter,
  UnivariateParameterNState,
  UnivariateParameterState,
  UnvAMWGState,
  UnvHMCState,
  UnvMALAState,
  UnvSMMALAState,
  VanillaMCTuner,
  Variable,
  VariableIOStream,
  VariableNState,
  VariableState,
  VariableStateVector,
  VariableVector,

  ### Functions
  acceptance,
  add_dimension,
  add_edge!,
  add_vertex!,
  codegen,
  count!,
  covariance!,
  dataset,
  datasets,
  diagnostics,
  edge_index,
  edges,
  erf_rate_score,
  ess,
  examples,
  failprob,
  iact,
  in_degree,
  in_edges,
  in_neighbors,
  indices,
  is_directed,
  is_indexed,
  iterate,
  job2dot,
  likelihood_model,
  logistic,
  logistic_rate_score,
  lognormalise,
  logpdf,
  lzv,
  make_edge,
  mcse,
  mcvar,
  model2dot,
  multivecs,
  normalise,
  num_edges,
  num_vertices,
  out_degree,
  out_edges,
  out_neighbors,
  output,
  params,
  pdf,
  qzv,
  rate!,
  reset!,
  reset_burnin!,
  revedge,
  run,
  sampler_state,
  save!,
  save,
  setpdf!,
  setprior!,
  softabs,
  sort_by_index,
  source,
  succprob,
  target,
  tune!,
  tuner_state,
  vertex_index,
  vertex_key,
  vertices

include("base.jl")

include("format.jl")

include("data.jl")

include("codegen.jl")

include("stats/logistic.jl")

include("distributions/Binary.jl")
include("distributions/TruncatedNormal.jl")

include("autodiff/reverse.jl")
include("autodiff/forward.jl")

include("states/VariableStates.jl")
include("states/ParameterStates/ParameterStates.jl")
include("states/ParameterStates/BasicDiscUnvParameterState.jl")
include("states/ParameterStates/BasicDiscMuvParameterState.jl")
include("states/ParameterStates/BasicContUnvParameterState.jl")
include("states/ParameterStates/BasicContMuvParameterState.jl")

include("nstates/VariableNStates.jl")
include("nstates/ParameterNStates/ParameterNStates.jl")
include("nstates/ParameterNStates/BasicDiscUnvParameterNState.jl")
include("nstates/ParameterNStates/BasicDiscMuvParameterNState.jl")
include("nstates/ParameterNStates/BasicContUnvParameterNState.jl")
include("nstates/ParameterNStates/BasicContMuvParameterNState.jl")

include("iostreams/VariableIOStreams.jl")
include("iostreams/ParameterIOStreams/ParameterIOStreams.jl")
include("iostreams/ParameterIOStreams/BasicDiscParamIOStream.jl")
include("iostreams/ParameterIOStreams/BasicContParamIOStream.jl")

include("variables/variables.jl")
include("variables/parameters/parameters.jl")
include("variables/parameters/BasicDiscUnvParameter.jl")
include("variables/parameters/BasicDiscMuvParameter.jl")
include("variables/parameters/BasicContUnvParameter.jl")
include("variables/parameters/BasicContMuvParameter.jl")
include("variables/dependencies.jl")

include("models/GenericModel.jl")
include("models/generators.jl")

include("ranges/ranges.jl")
include("ranges/BasicMCRange.jl")

include("tuners/tuners.jl")
include("tuners/VanillaMCTuner.jl")
include("tuners/AcceptanceRateMCTuner.jl")
include("tuners/DualAveragingMCTuner.jl")

include("samplers/samplers.jl")
include("samplers/ARS.jl")
include("samplers/MH.jl")
include("samplers/AM.jl")
include("samplers/RAM.jl")
include("samplers/HMC.jl")
# include("samplers/NUTS.jl")
include("samplers/MALA.jl")
include("samplers/SMMALA.jl")

include("jobs/jobs.jl")
include("jobs/BasicMCJob.jl")
include("jobs/BasicGibbsJob.jl")

include("samplers/iterate/ARS.jl")
include("samplers/iterate/MH.jl")
include("samplers/iterate/AM.jl")
include("samplers/iterate/RAM.jl")
include("samplers/iterate/HMC.jl")
# include("samplers/iterate/NUTS.jl")
include("samplers/iterate/MALA.jl")
include("samplers/iterate/SMMALA.jl")
include("samplers/iterate/iterate.jl")

include("stats/acceptance.jl")

include("stats/mean.jl")

include("stats/variance/mcvar.jl")
include("stats/variance/zv.jl")

include("stats/covariance.jl")

include("stats/convergence/ess.jl")
include("stats/convergence/iact.jl")

include("stats/metrics.jl")

end


text:
function gc_to_comp_sites(sites::SiteParam,groups::GroupParam)
    #given some groups and some sites calculated over those groups
    #calculates "flattened" sites
    comps = groups.components
    sitenames = deepcopy(sites.sites)
    gcnames = groups.flattenedgroups
    for i in 1:length(sitenames)
        gci = gcnames[i]
        sitei = sitenames[i]
        for a in 1:length(sitei)
            sitei[a] = gci * '{' * sitei[a] * '}'
        end
    end

    flattened_comp_sitenames = collect(Iterators.flatten(sitenames))
    flattened_comp_nsites = collect(Iterators.flatten(sites.n_sites))
    siteidx = deepcopy(sites.i_sites)
    for i in 1:length(siteidx)
        sitei = siteidx[i]
        for a in 1:length(sitei)
            sitei[a] = i
        end
    end
    flattened_gcidx = collect(Iterators.flatten(siteidx))
    flattened_comp_isites = collect(Iterators.flatten(sites.i_sites))
    idxdict = Dict((flattened_gcidx[i],flattened_comp_isites[i]) => i for i in 1:length(flattened_gcidx) )

    n_flattenedsites = Vector{Vector{Int}}(undef,length(comps))
    flat_groups = groups.n_flattenedgroups
    for i in 1:length(n_flattenedsites)
        n_flattenedsites[i] = flat_groups[i][flattened_gcidx] .* flattened_comp_nsites
    end
    k = length(flattened_comp_sitenames)
    pairs = [(comps[i],[flattened_comp_sitenames[j] =>n_flattenedsites[i][j] for j in 1:k]) for i in 1:length(comps)]
    n_flattenedsites

    for pair in pairs
        pairname,pairvec = pair
        filter!(x->!iszero(last(x)),pairvec)
    end

    return SiteParam(pairs),idxdict
end

#returns a compressed assoc matrix corresponding to the indices
#of the old gc values, arranged by component instead.
#that is,
function gc_to_comp_assoc_idx(param::AssocParam,sites::SiteParam,idxdict)
    pvals = param.values
    vals,outer,inner = pvals.values, pvals.outer_indices,pvals.inner_indices
    ngc = length(vals)
    nsites = length(sites.flattenedsites)
    comps = sites.components
    ncomps = length(comps)
    site1 = Vector{Int64}(undef,ngc)
    site2 = Vector{Int64}(undef,ngc)
    assoc_idxdict = Dict{Tuple{Int,Int},Int}()

    for ii in 1:ngc
        i,j = outer[ii]
        a,b = inner[ii]
        s1 =  idxdict[(i,a)]
        s2 = idxdict[(j,b)]
        site1[ii] = s1
        site2[ii] = s2
        assoc_idxdict[(s1,s2)] = ii
        assoc_idxdict[(s2,s1)] = ii
    end

    #return site1,site2
    x = Matrix{Matrix{Int}}(undef,ncomps,ncomps)
    for i = 1:ncomps
        xi = zeros(Int,nsites,nsites)
        for (a,b) in Iterators.product(sites.i_sites[i],sites.i_sites[i])
            abval = get(assoc_idxdict,(a,b),0)

            xi[a,b] = abval
            xi[b,a] = abval

        end
        x[i,i] = xi
        for j = 1:i-1
            xi = zeros(Int,nsites,nsites)
            for (a,b) in Iterators.product(sites.i_sites[i],sites.i_sites[j])
                abval = get(assoc_idxdict,(a,b),0)

                xi[a,b] = abval
                xi[b,a] = abval
            end
            x[i,j] = xi
            x[j,i] = xi
        end
    end
    val =  Compressed4DMatrix(x)
    new_inneridx = val.inner_indices
    for (idx,(i,j),(a,b)) ∈ indices(val)
        _a = sites.i_flattenedsites[i][a]
        _b = sites.i_flattenedsites[j][b]
        new_inneridx[idx] = (_a,_b)
    end

    n = findall(!=((0,0)),new_inneridx)
    values = val.values[n]
    outer_indices = val.outer_indices[n]
    inner_indices = new_inneridx[n]
    outer_size = val.outer_size
    inner_size = val.inner_size
    newvals = Compressed4DMatrix(values,outer_indices,inner_indices,outer_size,inner_size)
end

ext: jl


text:
using XGBoost

dtrain = DMatrix("../data/agaricus.txt.train")
dtest = DMatrix("../data/agaricus.txt.test")

param = ["max_depth" => 2,
         "eta" => 1,
         "silent" => 0,
         "objective" => "binary:logistic"]
watchlist  = [(dtest,"eval"), (dtrain,"train")]
num_round = 3

bst = xgboost(dtrain, num_round, param = param, watchlist = watchlist)

print ("start testing prediction from first n trees\n")
labels = get_info(dtest, "label")

### predict using first 1 tree
pred1 = predict(bst, dtest, ntree_limit = 1)
# by default, we predict using all the trees
pred2 = predict(bst, dtest)

print("error of pred1= ", sum((pred1 .> 0.5) .!= labels) / float(size(pred1)[1]), "\n")
print("error of pred2= ", sum((pred2 .> 0.5) .!= labels) / float(size(pred2)[1]), "\n")


text:
@time @testset "MC: IsingModel Simulation" begin
    Random.seed!(123)
    m = IsingModel(dims=2, L=8);
    mc = MC(m, beta=0.35, sweeps=1000, thermalization=10);
    run!(mc, verbose=false);

    # Check measurements
    measured = measurements(mc)
    if VERSION.major == 1 && VERSION.minor == 5 # v1.5.x
        @test   25.87  ≈ measured[:Magn].M |> mean         atol=0.01
        @test    0.82  ≈ measured[:Magn].M |> std_error    atol=0.01
        @test  913.    ≈ measured[:Magn].M2 |> mean        atol=1.0
        @test   47.    ≈ measured[:Magn].M2 |> std_error   atol=1.0
        @test    0.404 ≈ measured[:Magn].m |> mean         atol=0.001
        @test    0.013 ≈ measured[:Magn].m |> std_error    atol=0.001
        @test    1.331 ≈ measured[:Magn].chi |> mean       atol=0.001

        @test  -59.55  ≈ measured[:Energy].E |> mean       atol=0.01
        @test    1.05  ≈ measured[:Energy].E |> std_error  atol=0.01
        @test 3896.    ≈ measured[:Energy].E2 |> mean      atol=1.0
        @test  136.    ≈ measured[:Energy].E2 |> std_error atol=1.0
        @test   -0.930 ≈ measured[:Energy].e |> mean       atol=0.001
        @test    0.016 ≈ measured[:Energy].e |> std_error  atol=0.001
        @test    0.66  ≈ measured[:Energy].C |> mean       atol=0.01
    else # assuming "v1.6"
        @test   25.61  ≈ measured[:Magn].M |> mean         atol=0.01
        @test    0.82  ≈ measured[:Magn].M |> std_error    atol=0.01
        @test  903.    ≈ measured[:Magn].M2 |> mean        atol=1.0
        @test   47.    ≈ measured[:Magn].M2 |> std_error   atol=1.0
        @test    0.400 ≈ measured[:Magn].m |> mean         atol=0.001
        @test    0.013 ≈ measured[:Magn].m |> std_error    atol=0.001
        @test    1.354 ≈ measured[:Magn].chi |> mean       atol=0.001

        @test  -59.10  ≈ measured[:Energy].E |> mean       atol=0.01
        @test    0.97  ≈ measured[:Energy].E |> std_error  atol=0.01
        @test 3838.    ≈ measured[:Energy].E2 |> mean      atol=1.0
        @test  125.    ≈ measured[:Energy].E2 |> std_error atol=1.0
        @test   -0.924 ≈ measured[:Energy].e |> mean       atol=0.001
        @test    0.015 ≈ measured[:Energy].e |> std_error  atol=0.001
        @test    0.66  ≈ measured[:Energy].C |> mean       atol=0.01
    end
    @test isempty(mc.configs) == true
end

@time @testset "DQMC: attractive HubbardModel Simulation" begin
    Random.seed!(123)
    m = HubbardModelAttractive(4, 2);
    N = 4*4
    mc = DQMC(m, beta=2.0, sweeps=1000, thermalization=1000, measure_rate=1);
    mc[:G]    = greens_measurement(mc, m)
    mc[:CDC]  = charge_density_correlation(mc, m)
    mc[:SDCx] = spin_density_correlation(mc, m, :x)
    mc[:SDCy] = spin_density_correlation(mc, m, :y)
    mc[:SDCz] = spin_density_correlation(mc, m, :z)
    mc[:PC]   = pairing_correlation(mc, m, K=5)
    @time run!(mc, verbose=false);

    # Check measurements
    measured = measurements(mc)
    atol = 0.05

    # hardcoded means from 10k + 100k sweeps at delta_tau = 0.05

    # Greens
    @test check([
        0.498957, -0.154395, 0.000053, -0.154416, -0.154154, 0.000141, 0.034410, 0.000210, 0.000108, 0.034333, -0.000200, 0.034306, -0.154648, -0.000083, 0.034208, 0.000085, -0.152890, 0.498639, -0.153695, 0.000091, 0.000273, -0.154036, -0.000050, 0.033964, 0.033928, 0.000006, 0.033996, -0.000136, 0.000074, -0.154045, 0.000229, 0.034132, 0.000128, -0.154022, 0.499580, -0.153634, 0.034171, 0.000025, -0.154078, 0.000146, -0.000004, 0.034092, 0.000131, 0.034072, 0.034245, -0.000105, -0.153440, -0.000119, -0.153696, 0.000080, -0.154157, 0.501049, 0.000086, 0.034123, -0.000070, -0.154196, 0.034125, -0.000081, 0.034150, 0.000055, -0.000041, 0.034279, 0.000088, -0.154186, -0.153264, 0.000371, 0.034060, 0.000245, 0.498036, -0.154372, 0.000155, -0.153932, -0.153608, 0.000092, 0.034082, 0.000201, 0.000232, 0.034094, -0.000321, 0.034104, 0.000262, -0.153699, 0.000081, 0.034058, -0.153528, 0.500891, -0.154033, 0.000206, 0.000059, -0.153755, 0.000289, 0.034010, 0.034177, -0.000105, 0.033938, -0.000036, 0.033740, 0.000082, -0.153234, -0.000050, 0.000050, -0.153561, 0.499599, -0.153333, 0.033883, -0.000045, -0.152781, -0.000036, -0.000049, 0.034155, 0.000227, 0.034133, 0.000148, 0.034225, -0.000064, -0.153766, -0.153904, 0.000017, -0.154348, 0.499199, -0.000076, 0.034137, 0.000148, -0.153761, 0.034247, 0.000029, 0.034081, 0.000019, 0.000052, 0.034184, 0.000035, 0.034163, -0.153797, 0.000044, 0.034215, 0.000005, 0.500607, -0.154214, 0.000131, -0.154038, -0.154267, -0.000143, 0.034096, -0.000078, 0.033862, -0.000138, 0.034106, -0.000017, 0.000144, -0.154080, 0.000050, 0.034031, -0.153451, 0.498722, -0.153177, 0.000093, -0.000084, -0.153991, 0.000129, 0.034177, -0.000130, 0.034281, 0.000006, 0.034198, 0.034250, 0.000166, -0.154725, 0.000197, 0.000084, -0.154542, 0.498484, -0.154350, 0.034178, -0.000020, -0.154121, -0.000026, 0.034039, -0.000125, 0.034244, -0.000080, 0.000145, 0.034210, 0.000123, -0.154070, -0.153918, 0.000131, -0.153565, 0.500541, 0.000031, 0.034289, 0.000287, -0.154012, -0.153298, 0.000064, 0.034199, 0.000129, 0.000237, 0.034193, -0.000130, 0.034159, -0.153641, -0.000098, 0.034004, 0.000039, 0.499920, -0.153995, 0.000098, -0.154059, 0.000059, -0.154089, -0.000253, 0.034143, 0.034101, 0.000048, 0.034257, 0.000005, -0.000129, -0.154156, 0.000003, 0.034052, -0.153945, 0.501709, -0.153576, -0.000164, 0.034110, 0.000164, -0.154301, 0.000288, -0.000214, 0.034268, 0.000260, 0.034275, 0.034270, 0.000093, -0.153669, 0.000238, 0.000065, -0.154439, 0.499976, -0.154191, 0.000088, 0.034170, -0.000178, -0.153559, 0.034077, 0.000021, 0.034207, -0.000004, -0.000100, 0.034081, -0.000005, -0.153590, -0.154030, -0.000257, -0.153332, 0.498969
    ], (measured[:G] |> mean)[:], atol)
    @test check([
        0.01903630965131304, 0.005367516479841765, 0.0024584912434881534, 0.004890991356523148, 0.006363727124689783, 0.002670798811073267, 0.0021066937139101443, 0.00301559391676348, 0.003194080063641578, 0.0013442384852690233, 0.0021844921675267036, 0.0017208532918349584, 0.004902256544938516, 0.003022457879188603, 0.001959177525541325, 0.0023162054411693, 0.004275338756144313, 0.016190447778031265, 0.00622439056707409, 0.0027891086921042144, 0.0024747787231799347, 0.004743929829605756, 0.0026022570507337954, 0.0018811623274672816, 0.0020296410229540124, 0.001917007923316469, 0.0017821822872064257, 0.0015874651085898946, 0.002073785126319732, 0.005203727105701843, 0.0024428864285170842, 0.0014482403743486525, 0.002732283686976738, 0.006241513316141949, 0.017947118144596073, 0.0066844610214952535, 0.002073388696305197, 0.0026748418574992702, 0.006343349608775392, 0.002119334209726534, 0.0019276798962566477, 0.001912086325670065, 0.0027078659347968587, 0.0017274987770035345, 0.001730730946599209, 0.002755975555112104, 0.0062626591213190445, 0.002538493010301427, 0.006098666836121245, 0.002351234708256203, 0.007537561825151221, 0.01482652736208729, 0.002226492682037637, 0.0018870261129214604, 0.0028813646224384037, 0.0056738286375110775, 0.001978694193790825, 0.0015194200298218336, 0.0018930490650442217, 0.002954208794234731, 0.002539403264143692, 0.002043787453315343, 0.002465116172071212, 0.00633723160781742, 0.005152182403326165, 0.00221980109044263, 0.00209307523188663, 0.002279203611722848, 0.01727462737987002, 0.006193953525747859, 0.0022950804203056424, 0.00605660516317818, 0.0056744842883902765, 0.0021411097877753606, 0.0015106572026585446, 0.002285965572601132, 0.0025851060623688572, 0.0019413813197430884, 0.0018293773252474807, 0.0018633913153026519, 0.0023424483738403476, 0.00560666781142671, 0.0028486833567268675, 0.0018996697083859427, 0.006778416349325631, 0.015529197915268668, 0.005802396675855519, 0.002577418620315891, 0.002537373388636465, 0.004267295023739163, 0.002216982496441896, 0.0014470288731124991, 0.0015448157761172826, 0.0026749909170203207, 0.0015960502414737114, 0.0018638971133579443, 0.002176172975333444, 0.0024661226628342214, 0.0059976529594891895, 0.0024275919659689023, 0.0024874460957389284, 0.004866152879901385, 0.016989888842258615, 0.00606786788834373, 0.0023650404715101145, 0.0022059530709299794, 0.006279985224570244, 0.0022591078220764883, 0.002244262826955961, 0.0017607727293324292, 0.002805151032378269, 0.0019203802814192405, 0.002808211640591048, 0.0014493790069284944, 0.0026815107294531825, 0.005244751689767846, 0.004457555716976892, 0.0027461862709045517, 0.007432682630957972, 0.01755180302734928, 0.0027546143757098953, 0.0013564156087778078, 0.002707184989403521, 0.005099207586692647, 0.00150711950287787, 0.0018469456135621182, 0.0021664965287707533, 0.0028573759197268573, 0.002257562058149153, 0.0018752987165948749, 0.0018474643382295432, 0.0015607159782713032, 0.004821113270875764, 0.002364393502281104, 0.0015914287961616067, 0.002017966865185779, 0.016678090766298847, 0.005532950615319388, 0.0025109463649847317, 0.00490999391748333, 0.005540189949851306, 0.0025778872950522782, 0.0015297393649903565, 0.0021526364720696263, 0.0015571786088503922, 0.0030715698656264088, 0.0019442247942314111, 0.001959609115393741, 0.002451806750645247, 0.00559182183037939, 0.002364438409286417, 0.002033163891428257, 0.007385122100593753, 0.016250001923804477, 0.00588029391580204, 0.0022664676364664562, 0.0022713083805696375, 0.004920463385631935, 0.002555556649254254, 0.0024253973749420706, 0.0016998977444620306, 0.0023730453964623094, 0.0026884831369355915, 0.002052298758092647, 0.0017405567697183519, 0.0025459146940258817, 0.006477775129878893, 0.0027559441428842125, 0.002507883713380407, 0.005579200483908438, 0.015149687637628374, 0.005716157182003588, 0.0020198143514409796, 0.002678143865896058, 0.005627834489234238, 0.0019505457843814742, 0.0015667358381430906, 0.0019693925268430904, 0.0017634657213137022, 0.0027309321641583004, 0.0023192160020215594, 0.001492816198492005, 0.002785975799730211, 0.005637866774905859, 0.007373646876499863, 0.0024441882113822207, 0.005057050695680436, 0.016981962144626644, 0.0026358465926791313, 0.0021629519326452178, 0.001978862212401592, 0.00559129466697721, 0.0053878998822222395, 0.002378047860973303, 0.0021696210935349316, 0.0021858413712564696, 0.002603168942043464, 0.0018142457233731126, 0.0019318206065561722, 0.0016916601775288252, 0.006462811051987628, 0.0023964206261950405, 0.00201234725665299, 0.0020856386696421355, 0.01818721557230035, 0.0078336450885551, 0.002879127598649167, 0.005559421935709471, 0.0025272247116726937, 0.006749213317453175, 0.002407881174684203, 0.001723770533529353, 0.001685814864352903, 0.002595017712103584, 0.001784867442122867, 0.0016512832001637893, 0.0029511346895967168, 0.004831990577677809, 0.002363895478273945, 0.0017365815749286698, 0.0061205221823493194, 0.017463498150406046, 0.006898041663060922, 0.002930549060657847, 0.0013097885169025025, 0.002477957003239124, 0.005041328409858466, 0.0022703149157202326, 0.0017091772829437746, 0.0014423544812366467, 0.0026124944984324892, 0.0017586921671544939, 0.0018143303377610813, 0.001981021736540454, 0.005538479554528282, 0.0015071193020268251, 0.0023646394774796474, 0.005869543468904877, 0.01732179843958681, 0.004853138482281587, 0.0019211890883876162, 0.0019921005132619825, 0.0029311927912516376, 0.005390881176746915, 0.0014190304988554389, 0.0018902190049541409, 0.0023298593140222215, 0.0024414631574547737, 0.0025022457477321563, 0.001662072259828581, 0.0022766455352105424, 0.0052142514926209705, 0.005934479064251144, 0.002908842935009184, 0.006626344451829513, 0.017266063739950332
    ], (measured[:G] |> std_error)[:], atol)

    # Charge Density Correlation
    @test check([
        1.559723, 0.947929, 0.947920, 0.947929, 0.947920, 1.002382, 1.002578, 1.002382, 1.002578, 1.002973, 1.002722, 0.998367, 0.998192, 0.998192, 0.998367, 1.002403
    ], mean(measured[:CDC]), atol)
    @test check([
        0.013869638657504228, 0.013470915194620922, 0.013282700278496358, 0.013470915194621087, 0.013282700278496358, 0.013550942726610091, 0.014307471098069393, 0.014492552221530538, 0.014129517405169136, 0.01449255222153069, 0.014307471098069393, 0.013948472512305249, 0.01394847251230509, 0.01389084520478554, 0.01389084520478546, 0.013133272694007561
    ], std_error(measured[:CDC]), atol)

    # Spin density correlations (x, y, z)
    @test check([
        0.441557, -0.043403, -0.043389, -0.043403, -0.043389, -0.000937, -0.000938, -0.000937, -0.000938, -0.000940, -0.000942, -0.002002, -0.002002, -0.002002, -0.002002, -0.000530
    ], mean(measured[:SDCx]), atol)
    @test check([
        0.0015350907401287547, 0.0002368878415500336, 0.0002607507287443457, 0.0002368878415500336, 0.00026075072874433736, 6.64636472973452e-5, 2.9316744933934093e-5, 3.521865202114811e-5, 6.282380159287888e-5, 3.521865202114811e-5, 2.931674493393406e-5, 2.5828399958572843e-5, 2.5828399958573172e-5, 2.7096518574704045e-5, 2.7096518574704045e-5, 2.605149773144738e-5
    ], std_error(measured[:SDCx]), atol)
    @test check([
        0.441557, -0.043403, -0.043389, -0.043403, -0.043389, -0.000937, -0.000938, -0.000937, -0.000938, -0.000940, -0.000942, -0.002002, -0.002002, -0.002002, -0.002002, -0.000530
    ], mean(measured[:SDCy]), atol)
    @test check([
        0.0015350907401287547, 0.0002368878415500336, 0.0002607507287443457, 0.0002368878415500336, 0.00026075072874433736, 6.64636472973452e-5, 2.9316744933934093e-5, 3.521865202114811e-5, 6.282380159287888e-5, 3.521865202114811e-5, 2.931674493393406e-5, 2.5828399958572843e-5, 2.5828399958573172e-5, 2.7096518574704045e-5, 2.7096518574704045e-5, 2.605149773144738e-5
    ], std_error(measured[:SDCy]), atol)
    @test check([
        0.441557, -0.043403, -0.043389, -0.043403, -0.043389, -0.000937, -0.000938, -0.000937, -0.000938, -0.000940, -0.000942, -0.002002, -0.002002, -0.002002, -0.002002, -0.000530
    ], mean(measured[:SDCz]), atol)
    @test check([
        0.0015350907401287547, 0.00023688784155001527, 0.00026075072874433736, 0.00023688784155001527, 0.0002607507287443457, 6.64636472973453e-5, 2.9316744933933842e-5, 3.521865202114802e-5, 6.282380159287893e-5, 3.521865202114805e-5, 2.931674493393381e-5, 2.58283999585735e-5, 2.5828399958572843e-5, 2.7096518574704045e-5, 2.7096518574704357e-5, 2.605149773144732e-5
    ], std_error(measured[:SDCz]), atol)

    # Pairing Correlations
    @test check([
        0.278901, 0.026748, 0.026731, 0.026735, 0.026738, 0.000645, 0.000645, 0.000644, 0.000645, 0.000646, 0.000649, 0.001480, 0.001480, 0.001480, 0.001479, 0.000368, -0.076903, -0.076920, -0.000011, -0.000007, -0.000010, -0.000008, 0.000002, 0.000003, -0.000009, -0.000007, 0.000003, 0.000001, 0.000001, -0.000003, 0.000002, -0.000003, -0.076893, -0.000007, -0.076900, -0.000010, -0.000011, -0.000005, -0.000011, 0.000003, 0.000002, 0.000002, -0.000012, 0.000001, 0.000001, 0.000003, -0.000003, -0.000003, -0.076920, -0.000007, -0.000008, -0.076903, -0.000009, 0.000001, -0.000011, -0.000010, 0.000002, -0.000007, 0.000001, 0.000002, -0.000003, 0.000003, 0.000003, -0.000003, -0.076900, -0.000005, -0.000012, -0.000011, -0.076893, 0.000001, 0.000003, -0.000010, -0.000007, 0.000001, -0.000011, -0.000003, 0.000002, 0.000003, 0.000002, -0.000003, -0.076905, -0.000012, -0.000009, -0.076881, -0.000014, 0.000001, -0.000009, -0.000013, 0.000003, -0.000010, 0.000003, 0.000003, -0.000003, 0.000003, 0.000004, -0.000003, 0.247193, 0.022671, 0.023848, 0.022669, 0.023849, 0.000062, 0.000063, 0.000062, 0.000062, 0.000219, -0.000080, 0.001106, 0.001156, 0.001156, 0.001106, 0.000061, -0.001087, -0.005200, 0.022667, 0.022664, -0.005202, 0.000063, -0.001086, 0.000062, -0.000011, 0.000063, 0.000062, -0.005200, 0.001156, -0.005203, 0.001156, -0.000011, -0.001095, 0.023853, -0.005200, 0.021701, -0.005200, -0.000076, 0.000216, 0.000216, -0.000078, -0.001093, -0.000009, -0.005200, 0.001215, 0.001106, -0.005200, -0.000009, -0.001098, -0.005201, -0.005202, 0.022667, 0.022666, -0.000010, 0.000061, -0.001088, 0.000063, 0.000062, 0.000061, 0.001156, 0.001156, -0.005203, -0.005200, -0.000010, -0.076892, -0.000006, -0.000013, -0.000012, -0.076889, 0.000001, 0.000002, -0.000015, -0.000011, 0.000002, -0.000012, -0.000003, 0.000002, 0.000003, 0.000003, -0.000003, -0.001087, 0.022671, -0.005204, -0.005202, 0.022667, 0.000061, -0.000011, 0.000061, -0.001084, 0.000062, 0.000061, 0.001156, -0.005202, 0.001155, -0.005202, -0.000011, 0.247187, 0.023860, 0.022668, 0.023852, 0.022665, 0.000061, 0.000061, 0.000061, 0.000061, -0.000080, 0.000217, 0.001156, 0.001106, 0.001106, 0.001155, 0.000061, -0.001096, -0.005202, -0.005202, 0.022670, 0.022668, -0.000010, 0.000061, -0.001083, 0.000063, 0.000062, 0.000062, 0.001156, 0.001156, -0.005203, -0.005200, -0.000010, -0.001096, -0.005203, 0.023855, -0.005203, 0.021695, -0.000079, -0.000079, 0.000216, 0.000216, -0.000010, -0.001094, 0.001214, -0.005203, -0.005202, 0.001105, -0.000010, -0.076881, -0.076905, -0.000009, -0.000010, -0.000013, -0.000009, 0.000003, 0.000004, -0.000014, -0.000012, 0.000003, 0.000001, 0.000003, -0.000003, 0.000003, -0.000003, -0.001093, 0.021701, -0.005200, 0.023853, -0.005200, 0.000216, -0.000076, -0.000078, 0.000216, -0.001095, -0.000009, -0.005200, 0.001106, 0.001215, -0.005200, -0.000009, -0.001088, 0.022666, 0.022667, -0.005200, -0.005203, -0.001098, 0.000062, -0.000010, 0.000061, 0.000063, 0.000061, -0.005201, -0.005202, 0.001156, 0.001156, -0.000010, 0.247193, 0.022671, 0.023848, 0.022669, 0.023849, 0.000062, 0.000063, 0.000062, 0.000062, 0.000219, -0.000080, 0.001106, 0.001156, 0.001156, 0.001106, 0.000061, -0.001086, 0.022667, -0.005203, -0.005200, 0.022664, 0.000062, -0.000011, 0.000063, -0.001087, 0.000063, 0.000062, 0.001156, -0.005202, 0.001156, -0.005200, -0.000011, -0.076889, -0.000011, -0.076892, -0.000015, -0.000012, -0.000006, -0.000012, 0.000003, 0.000002, 0.000003, -0.000013, 0.000002, 0.000001, 0.000002, -0.000003, -0.000003, -0.001083, 0.022668, 0.022670, -0.005200, -0.005203, -0.001096, 0.000062, -0.000010, 0.000062, 0.000063, 0.000061, -0.005202, -0.005202, 0.001156, 0.001156, -0.000010, -0.001094, -0.005203, 0.021695, -0.005202, 0.023855, 0.000216, 0.000216, -0.000079, -0.000079, -0.000010, -0.001096, 0.001105, -0.005203, -0.005203, 0.001214, -0.000010, -0.001084, -0.005202, 0.022671, 0.022667, -0.005202, 0.000062, -0.001087, 0.000061, -0.000011, 0.000061, 0.000061, -0.005202, 0.001156, -0.005204, 0.001155, -0.000011, 0.247187, 0.023860, 0.022668, 0.023852, 0.022665, 0.000061, 0.000061, 0.000061, 0.000061, -0.000080, 0.000217, 0.001156, 0.001106, 0.001106, 0.001155, 0.000061
    ], mean(measured[:PC])[:], atol)
    @test check([
        0.0032863968694560874, 0.00042772252675633826, 0.0003518557357542367, 0.0004743615772443992, 0.0003187801835120089, 4.807866876266873e-5, 2.7443632781155007e-5, 3.087921747910029e-5, 4.295229784612734e-5, 2.7131480773117706e-5, 2.3360145651071746e-5, 4.6038701058516214e-5, 5.688082298878186e-5, 3.477280028969865e-5, 3.0122652022979076e-5, 2.028509759307593e-5, 0.0006793221509329129, 0.0006754871522655642, 0.00011843270443530931, 0.00011834451604052362, 0.00010615690256295016, 0.0001243664358098925, 0.00011105578571808066, 0.00012160922087931815, 2.8535027637957004e-5, 2.7566801509094677e-5, 2.1371276153817884e-5, 2.6505038538425558e-5, 2.4270303829044734e-5, 2.7585721128230576e-5, 3.5578032465851165e-5, 3.792786779387757e-5, 0.0006820245874155875, 0.00011194387420314651, 0.0006431256216668431, 0.00010798281383960798, 0.00011839952161819387, 2.747026452012285e-5, 0.00012140889313203758, 2.6049298318322218e-5, 0.00011433387417751547, 0.00012556056981436842, 2.3905752587305103e-5, 2.509640204375193e-5, 3.635761848562791e-5, 2.465954105096387e-5, 2.5137728666196422e-5, 3.594115483339847e-5, 0.000675487152265577, 0.0001243664358098925, 0.00011105578571808066, 0.0006793221509329129, 0.00012160922087931817, 0.00011834451604052362, 2.6505038538425558e-5, 2.4270303829044734e-5, 2.7585721128230576e-5, 0.00011843270443530931, 0.00010615690256295016, 2.7566801509094677e-5, 2.1371276153817884e-5, 3.792786779387757e-5, 2.8535027637957007e-5, 3.5578032465851165e-5, 0.0006431256216668431, 0.00012140889313203759, 0.00011433387417751547, 0.00012556056981436842, 0.0006820245874155747, 2.509640204375193e-5, 2.465954105096387e-5, 0.00011194387420314651, 0.00011839952161819387, 2.5137728666196422e-5, 0.00010798281383960797, 3.594115483339846e-5, 2.7470264520122846e-5, 2.6049298318322218e-5, 2.3905752587305103e-5, 3.6357618485627904e-5, 0.0005712777277038281, 0.00011138201020643361, 0.00011547593940896772, 0.0006185946733805882, 0.0001250322089856034, 0.00011099025943433941, 2.2911566901869562e-5, 2.7069997775801493e-5, 2.547577645224039e-5, 0.00011394307751195347, 0.00011775034186131933, 2.516833737154611e-5, 2.5865813768100625e-5, 3.5584130541739836e-5, 2.4073549197964844e-5, 3.3867429395236943e-5, 0.003539859671041112, 0.00016362767390622782, 0.0002246191699138663, 0.00017845343617490725, 0.00021268203128181403, 3.8944094417145216e-5, 1.1728571613108348e-5, 1.1083582357649566e-5, 3.267501471471958e-5, 1.2586212552108007e-5, 1.0841571469820574e-5, 2.0426518109442387e-5, 2.0692330560425303e-5, 2.198816223405185e-5, 2.075131248696356e-5, 1.4831088536320838e-5, 0.0003515926250209984, 8.364142899254581e-5, 0.00014397768229978744, 0.00015017652713828134, 6.957343027475577e-5, 1.2452188285585125e-5, 1.239432966071291e-5, 1.170251699251169e-5, 1.2510294872646639e-5, 0.00035642112923309986, 1.118906729995006e-5, 7.789029219399845e-5, 2.226663792917963e-5, 2.036813610782931e-5, 6.491932854948363e-5, 1.3179372854946754e-5, 0.0003739027082114623, 0.00024321517489145862, 6.247968326554638e-5, 0.00011844392077501222, 6.85233902257041e-5, 0.0003876668716336736, 1.8169742082622794e-5, 2.0000420131148687e-5, 1.6304784801120634e-5, 1.8452438168822835e-5, 1.9947241100376264e-5, 6.670410552204984e-5, 6.65542070420358e-5, 2.4596544124544127e-5, 1.96976597592882e-5, 1.6180256909043756e-5, 0.0003515376320593395, 7.232255697701556e-5, 7.068637292329068e-5, 0.00015234503654893702, 0.00015095582137596026, 1.4494089296318829e-5, 1.227868057837655e-5, 1.3777969885874737e-5, 1.3804986335903444e-5, 1.1115952291683665e-5, 0.0003455496357607849, 2.1224866990386695e-5, 7.8317894232434e-5, 2.0719827570342864e-5, 6.534385175842493e-5, 1.1001207962411203e-5, 0.0006235602912179505, 0.00011871663414839694, 0.00012212270640311318, 0.00011953805525290922, 0.000620221577511302, 2.3188118242046084e-5, 2.7692702639765495e-5, 0.00011079105994197718, 0.0001158178029217961, 2.4528820820853455e-5, 0.00012000417025520125, 3.5920923142131966e-5, 2.500081224729451e-5, 2.866707037909041e-5, 2.594548700031861e-5, 3.595985463636785e-5, 0.0003759489283794092, 0.00014814065865731806, 7.202950072343839e-5, 7.270595014281034e-5, 0.00014942849039939674, 1.4941393663445203e-5, 1.4004128039011912e-5, 0.0003352754974847287, 1.2861788315205813e-5, 1.5430611673346766e-5, 1.2180076096398397e-5, 2.246155755059947e-5, 7.752875372782713e-5, 7.116360379326014e-5, 2.172539102064846e-5, 1.220029529846814e-5, 0.003582411455132999, 0.0002825484222289856, 0.00018099728215045593, 0.0003012407141657386, 0.00017876709707790486, 3.915089723405503e-5, 1.2658245820587652e-5, 1.2771479152914011e-5, 3.2790573680904847e-5, 1.297033650205429e-5, 1.0730709603946599e-5, 2.76191232435535e-5, 2.8418952423800635e-5, 2.143839940531424e-5, 2.000220367457017e-5, 1.6179065936505336e-5, 0.0003489193143297205, 7.455741065724246e-5, 6.928355385121493e-5, 0.00015452264028963485, 0.00015484581765932366, 1.2996421657299983e-5, 1.3936851912543852e-5, 1.2758762469247603e-5, 1.3025773586153131e-5, 1.288541383164654e-5, 0.00032266596637020817, 2.168642281079922e-5, 7.852068318219601e-5, 2.2004304032514555e-5, 6.674635580332466e-5, 1.307995089540432e-5, 0.00037025964275628185, 7.432900732372818e-5, 0.00023730774991964098, 7.765883516354168e-5, 0.00013037536437217285, 2.0332637110704405e-5, 1.87160378916926e-5, 2.0473022691921813e-5, 0.0003660098090677224, 2.0992174207137846e-5, 2.155902062952222e-5, 3.961712150627199e-5, 2.3168964349629955e-5, 7.398702065306552e-5, 7.433944300255508e-5, 1.855683938692942e-5, 0.0006185946733805882, 0.0005712777277038281, 0.00011394307751195347, 0.00011099025943433941, 0.00011775034186131933, 0.00011138201020643361, 0.00011547593940896772, 0.0001250322089856034, 2.4073549197964844e-5, 2.516833737154611e-5, 2.5865813768100625e-5, 2.2911566901869562e-5, 2.7069997775801493e-5, 2.547577645224039e-5, 3.3867429395236943e-5, 3.5584130541739836e-5, 0.00038766687163367355, 0.0001184439207750168, 6.670410552204831e-5, 0.00024321517489145862, 6.655420704203732e-5, 0.0003739027082114623, 1.845243816882282e-5, 1.9947241100376264e-5, 1.6180256909043756e-5, 1.8169742082622794e-5, 2.0000420131148687e-5, 6.2479683265548e-5, 6.852339022570557e-5, 1.96976597592882e-5, 2.4596544124543955e-5, 1.6304784801120634e-5, 0.0003455496357607849, 0.0001509558213759674, 0.00015234503654893702, 7.8317894232434e-5, 6.53438517584244e-5, 1.3777969885874739e-5, 0.0003515376320593395, 1.3804986335903442e-5, 1.1115952291683665e-5, 1.449408929631883e-5, 1.1001207962411203e-5, 7.232255697701509e-5, 2.0719827570343067e-5, 7.068637292328972e-5, 2.1224866990386895e-5, 1.227868057837655e-5, 0.003539859671041112, 0.00016362767390622782, 0.00022461916991386148, 0.0001784534361749194, 0.00021268203128181403, 3.894409441714523e-5, 1.1728571613108352e-5, 1.1083582357649566e-5, 3.2675014714719584e-5, 1.2586212552108007e-5, 1.0841571469820574e-5, 2.0426518109442387e-5, 2.0692330560425506e-5, 2.198816223405185e-5, 2.075131248696356e-5, 1.4831088536320838e-5, 0.00035642112923309986, 0.00014397768229979497, 6.491932854948468e-5, 7.789029219399845e-5, 0.00015017652713828134, 1.239432966071291e-5, 1.2510294872646642e-5, 0.0003515926250209984, 1.118906729995006e-5, 1.3179372854946755e-5, 1.2452188285585125e-5, 2.03681361078291e-5, 8.3641428992545e-5, 6.957343027475577e-5, 2.226663792917944e-5, 1.170251699251169e-5, 0.000620221577511302, 0.00011079105994197718, 0.0006235602912179505, 0.00012000417025520125, 0.0001158178029217961, 2.500081224729451e-5, 0.00011871663414839694, 2.866707037909041e-5, 0.00012212270640311318, 0.00011953805525290922, 2.594548700031861e-5, 2.3188118242046084e-5, 3.595985463636785e-5, 2.7692702639765495e-5, 2.4528820820853455e-5, 3.5920923142131966e-5, 0.0003226659663702083, 0.00015484581765932366, 0.0001545226402896208, 7.852068318219601e-5, 6.674635580332466e-5, 1.2758762469247605e-5, 0.0003489193143297205, 1.3025773586153133e-5, 1.2885413831646546e-5, 1.2996421657299983e-5, 1.307995089540432e-5, 7.455741065724246e-5, 2.2004304032514555e-5, 6.92835538512154e-5, 2.168642281079922e-5, 1.3936851912543852e-5, 0.00036600980906772247, 7.398702065306415e-5, 0.00013037536437217285, 7.433944300255508e-5, 0.0002373077499196364, 1.855683938692942e-5, 2.047302269192182e-5, 1.8716037891692603e-5, 0.00037025964275628185, 2.155902062952222e-5, 2.0992174207137846e-5, 2.3168964349629955e-5, 3.961712150627209e-5, 7.432900732372909e-5, 7.765883516354168e-5, 2.0332637110704405e-5, 0.0003352754974847287, 7.752875372782625e-5, 0.00014814065865731074, 0.000149428490399404, 7.11636037932587e-5, 1.2180076096398397e-5, 1.4941393663445203e-5, 1.2200295298468141e-5, 1.4004128039011908e-5, 0.00037594892837940913, 1.2861788315205815e-5, 7.270595014281173e-5, 2.172539102064846e-5, 2.246155755059975e-5, 7.202950072343933e-5, 1.5430611673346766e-5, 0.003582411455132999, 0.0002825484222289856, 0.00018099728215046192, 0.0003012407141657386, 0.00017876709707790486, 3.9150897234055034e-5, 1.2658245820587652e-5, 1.2771479152914011e-5, 3.279057368090484e-5, 1.297033650205429e-5, 1.0730709603946597e-5, 2.7619123243553578e-5, 2.841895242380056e-5, 2.143839940531424e-5, 2.0002203674570274e-5, 1.6179065936505336e-5
    ], std_error(measured[:PC])[:], atol)
end

@time @testset "DQMC: repulsive HubbardModel Simulation" begin
    Random.seed!(123)
    m = HubbardModelRepulsive(2, 2);
    mc = DQMC(m, beta=2.0, sweeps=1000, thermalization=1000, measure_rate=1);
    mc[:G]    = greens_measurement(mc, m)
    mc[:CDC]  = charge_density_correlation(mc, m)
    mc[:Mx]   = magnetization(mc, m, :x)
    mc[:My]   = magnetization(mc, m, :y)
    mc[:Mz]   = magnetization(mc, m, :z)
    mc[:SDCx] = spin_density_correlation(mc, m, :x)
    mc[:SDCy] = spin_density_correlation(mc, m, :y)
    mc[:SDCz] = spin_density_correlation(mc, m, :z)
    mc[:PC]   = pairing_correlation(mc, m, K=3)
    run!(mc, verbose=false);

    # Check measurements
    measured = measurements(mc)
    atol = 0.05

    # Greens
    @test [
        0.499896, -0.249412, -0.248547, 0.001101, 0.000000, 0.000000, 0.000000, 0.000000, -0.249144, 0.501226, -0.000068, -0.249244, 0.000000, 0.000000, 0.000000, 0.000000, -0.249515, -0.000893, 0.499410, -0.249425, 0.000000, 0.000000, 0.000000, 0.000000, 0.000720, -0.249229, -0.247895, 0.498638, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.500104, -0.249144, -0.249515, -0.000720, 0.000000, 0.000000, 0.000000, 0.000000, -0.249412, 0.498774, 0.000893, -0.249229, 0.000000, 0.000000, 0.000000, 0.000000, -0.248547, 0.000068, 0.500590, -0.247895, 0.000000, 0.000000, 0.000000, 0.000000, -0.001101, -0.249244, -0.249425, 0.501362
    ] ≈ (measured[:G] |> mean)[:]  atol = 0.1
    @test [
        0.015558677710177986 0.008330899904555777 0.008184425352547245 0.009154682712334604 0.0 0.0 0.0 0.0; 0.007410831061669652 0.01648129284257949 0.008625434463440258 0.007294508374789638 0.0 0.0 0.0 0.0; 0.006843413476789877 0.008445936622654146 0.013787211108014254 0.008058662653954283 0.0 0.0 0.0 0.0; 0.009844147766786396 0.008848693062785571 0.008838191580658811 0.020082028056368488 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.01555867771017795 0.007410831061669576 0.0068434134767897965 0.009844147766786278; 0.0 0.0 0.0 0.0 0.008330899904555801 0.01648129284257949 0.008445936622654148 0.008848693062785485; 0.0 0.0 0.0 0.0 0.008184425352547372 0.008625434463440267 0.013787211108014214 0.008838191580658834; 0.0 0.0 0.0 0.0 0.009154682712334624 0.007294508374789657 0.008058662653954276 0.020082028056368446
    ] ≈ measured[:G] |> std_error      atol = atol

    # Charge Density Correlation
    @test [
        1.442123, 0.884519, 0.884452, 0.976506
    ] ≈ mean(measured[:CDC]) atol = atol
    @test [
        0.002734862709277855, 0.0011120013039444745, 0.0010697404888432946, 0.0011271551926820696
    ][:] ≈ std_error(measured[:CDC]) atol = atol

    # Magnetization
    @test [
        0.0, 0.0, 0.0, 0.0
    ] ≈ mean(measured[:Mx])        atol = atol
    @test [
        0.0, 0.0, 0.0, 0.0
    ] ≈ std_error(measured[:Mx])   atol = atol
    @test [
        0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im
    ] ≈ mean(measured[:My])        atol = atol
    @test [
        0.0, 0.0, 0.0, 0.0
    ] ≈ std_error(measured[:My])   atol = atol
    @test [
        0.02, -0.02, -0.04, 0.00 # should this be 0 with error?
    ] ≈ mean(measured[:Mz]) atol = atol
    @test [
        0.0053949097720738, 0.007015652647519458, 0.007894887156482287, 0.007525414225664361
    ] ≈ std_error(measured[:Mz]) atol = atol

    # Spin density correlations (x, y, z)
    @test [
        0.557877, -0.139454, -0.139804, 0.030304
    ] ≈ mean(measured[:SDCx]) atol = atol
    @test [
        0.0027348627092750135, 0.002691219601002444, 0.002555831150492111, 0.0013440578305848455
    ][:] ≈ std_error(measured[:SDCx]) atol = atol
    @test [
        0.557877, -0.139454, -0.139804, 0.030304
    ] ≈ mean(measured[:SDCy]) atol = atol
    @test [
        0.0027348627092750135, 0.002691219601002444, 0.002555831150492111, 0.0013440578305848455
    ][:] ≈ std_error(measured[:SDCy]) atol = atol
    @test [
        0.557877, -0.139273, -0.139662, 0.030327
    ] ≈ mean(measured[:SDCz]) atol = atol
    @test [
        0.0027348627092748106, 0.00791852904864558, 0.007824969460855756, 0.006651107663992822
    ][:] ≈ std_error(measured[:SDCz]) atol = atol

    # Pairing Correlations
    @test [
        0.221061, 0.057740, 0.057774, -0.011747, -0.124435, -0.124531, 0.000064, -0.000029, -0.124539, 0.000081, -0.124630, -0.000032, -0.124437, -0.124534, 0.000029, -0.000064, 0.255948, 0.069727, 0.059859, -0.001775, 0.001648, 0.063452, 0.063422, 0.001848, -0.124504, 0.000032, -0.124595, -0.000081, 0.001632, 0.063422, 0.063452, 0.001863, 0.256029, 0.059847, 0.069902, -0.001776
    ] ≈ mean(measured[:PC])[:] atol = atol
    @test [
        0.001367431354637659, 0.0005560006519724946, 0.0005348702444204311, 0.0005635775963397751, 0.002483170888063298, 0.0019995138471234027, 0.0016765856228872162, 0.0017201184012394855, 0.002131747013786726, 0.0017425693779002111, 0.0022456690936882052, 0.0017307653829596891, 0.002243883064771306, 0.002225762938342551, 0.001720118401239482, 0.0016765856228872218, 0.001947254588881961, 0.001345609800501222, 0.0007977932487371304, 0.0007110834841664547, 0.003803528962932711, 0.0007363554205369239, 0.0007285186026716149, 0.003882312744872355, 0.002142487249272225, 0.0017307653829596863, 0.002077188962229575, 0.0017425693779002172, 0.0038833853421802026, 0.000728518602671597, 0.0007363554205368886, 0.00393584111062321, 0.001906193690791037, 0.0008143580947866242, 0.0012779155752460556, 0.0007228031157513561
    ] ≈ std_error(measured[:PC])[:] atol = atol
end


text:
#
# Abstract syntax tree.
#

# Syntax node.

immutable Syntax
    src::Any
    val::Any
    label::Symbol
    args::Vector{Syntax}
end

Syntax(src, val) =
    Syntax(src, val, NO_LABEL, NO_ARGS)

Syntax(src, label::Symbol, ::Void) =
    Syntax(src, nothing, label, NO_ARGS)

Syntax(src, label::Symbol, args::Vector{Syntax}) =
    Syntax(src, nothing, label, args)

const NO_LABEL = Symbol("")
const NO_ARGS = Syntax[]

haslabel(syn::Syntax) =
    syn.label != NO_LABEL

hasargs(syn::Syntax) =
    syn.args != NO_ARGS

source(syn::Syntax) = syn.src

value(syn::Syntax) = syn.val

label(syn::Syntax) = syn.label

start(syn::Syntax) = 1
next(syn::Syntax, st) = (syn.args[st], st+1)
done(syn::Syntax, st) = st > length(syn.args)
length(syn::Syntax) = length(syn.args)
isempty(syn::Syntax) = isempty(syn.args) : true

function show(io::IO, syn::Syntax)
    if syn.src !== nothing
        print(io, syn.src)
    elseif !haslabel(syn)
        print(io, syn.val)
    else
        print(io, syn.label)
        if hasargs(syn)
            print(io, "(")
            comma = false
            for arg in syn.args
                if comma
                    print(io, ", ")
                end
                comma = true
                print(io, arg)
            end
            print(io, ")")
        end
    end
end

Syntax(syn::Syntax) = syn

Syntax(ex) =
    expr2syntax(parse(string("(", ex, ")")))

# Parses a Julia expression as a Rabbit query.

expr2syntax(ex) =
    Syntax(ex, ex)

expr2syntax(ex::Symbol) =
    ex == :null ?
        Syntax(ex, nothing) :
        Syntax(ex, ex, nothing)

expr2syntax(ex::QuoteNode) =
    Syntax(ex.value, ex.value, nothing)

function expr2syntax(ex::Expr)
    if ex.head == :(:) && length(ex.args) >= 2
        head = length(ex.args) > 2 ?
            Expr(:(:), ex.args[1:end-1]...) :
            ex.args[1]
        syn = expr2syntax(head)
        call, tail = destructcolon(ex.args[end])
        syn =
            if isa(call, Symbol)
                Syntax(Expr(:(:), ex.args[1:end-1]..., call), call, [syn])
            elseif isa(call, Expr) && call.head == :call && isa(call.args[1], Symbol)
                args = Vector{Syntax}(length(call.args))
                args[1] = syn
                for k = 2:endof(args)
                    args[k] = expr2syntax(call.args[k])
                end
                Syntax(Expr(:(:), ex.args[1:end-1]..., call), call.args[1], args)
            else
                error("invalid query postfix notation: $ex")
            end
        if tail !== nothing
            syn = Syntax(ex, :(.), [syn, expr2syntax(tail)])
        end
        return syn
    elseif ex.head == :call && isa(ex.args[1], Symbol)
        args = Vector{Syntax}(length(ex.args)-1)
        for k in eachindex(args)
            args[k] = expr2syntax(ex.args[k+1])
        end
        return Syntax(ex, ex.args[1], args)
    elseif ex.head == :call
        return expr2syntax(destructcall(ex.args[1], ex.args[2:end]))
    elseif ex.head == :(.) || ex.head == :(=>)
        return Syntax(ex, ex.head, [expr2syntax(arg) for arg in ex.args])
    elseif ex.head == :(&&)
        return Syntax(ex, :(&), [expr2syntax(arg) for arg in ex.args])
    elseif ex.head == :(||)
        return Syntax(ex, :(|), [expr2syntax(arg) for arg in ex.args])
    elseif ex.head == :ref
        return Syntax(ex, :get, [expr2syntax(arg) for arg in ex.args])
    elseif ex.head == :vect
        return Syntax(ex, :array, [expr2syntax(arg) for arg in ex.args])
    elseif ex.head == :quote && length(ex.args) == 1
        return expr2syntax(ex.args[1])
    elseif ex.head == :comparison && length(ex.args) >= 3 && isa(ex.args[2], Symbol)
        syn = expr2syntax(Expr(:call, ex.args[2], ex.args[1], ex.args[3]))
        if length(ex.args) > 3
            syn = Syntax(ex, :(&), op, expr2syntax(Expr(:comparison, ex.args[3:end]...)))
        end
        return syn
    elseif ex.head == :macrocall && length(ex.args) == 2 && ex.args[1] == Symbol("@r_str")
        return Syntax(ex, eval(ex))
    else
        error("invalid query notation: $ex")
    end
end

destructcall(ex::QuoteNode, args) =
    Expr(:call, ex.value, args...)

destructcall(ex::Expr, args) =
    if ex.head == :quote && length(ex.args) == 1
        Expr(:call, ex.args[1], args...)
    else
        Expr(ex.head, ex.args[1:end-1]..., destructcall(ex.args[end], args))
    end

destructcall(ex, args) =
    error("invalid query call notation: $ex")

destructcolon(ex) =
    if isa(ex, Expr) && ex.head == :(.) && length(ex.args) == 2
        destructcolon(ex.args[1], ex.args[2])
    elseif (isa(ex, Expr) && ex.head == :call && length(ex.args) == 2 &&
            isa(ex.args[1], Expr) && ex.args[1].head == :(.) && length(ex.args[1].args) == 2)
        destructcolon(ex.args[1].args[1], Expr(:call, ex.args[1].args[2], ex.args[2]))
    else
        ex, nothing
    end

destructcolon(ex, tail) =
    if isa(ex, Expr) && ex.head == :(.) && length(ex.args) == 2
        destructcolon(ex.args[1], Expr(:(.), ex.args[2], tail))
    elseif (isa(ex, Expr) && ex.head == :call && length(ex.args) == 2 &&
            isa(ex.args[1], Expr) && ex.args[1].head == :(.) && length(ex.args[1].args) == 2)
        destructcolon(ex.args[1].args[1], Expr(:call, ex.args[1].args[2], Expr(:(.), ex.args[2], tail)))
    else
        ex, tail
    end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using Documenter
using vOptGeneric

makedocs(
    sitename = "vOptGeneric",
    format = Documenter.HTML(),
    modules = [vOptGeneric]
)

deploydocs(;
branch = "gh-pages",
    devbranch = "master",
    repo="github.com/pnavaro/vOptGeneric.jl",
)


text:
import Base: isempty, ∈, ∩

export Intersection,
       isempty_known,
       set_isempty!,
       swap,
       use_precise_ρ,
       constraints_list,
       linear_map

"""
    IntersectionCache

Container for information cached by a lazy `Intersection` object.

### Fields

- `isempty` -- is the intersection empty? There are three possible states,
               encoded as `Int8` values -1, 0, 1:

    * ``-1`` - it is currently unknown whether the intersection is empty or not
    *  ``0`` - intersection is not empty
    *  ``1`` - intersection is empty
"""
mutable struct IntersectionCache
    isempty::Int8

    # default constructor
    IntersectionCache() = new(Int8(-1))
end

function isempty_known(c::IntersectionCache)
    return c.isempty != Int8(-1)
end

function isempty(c::IntersectionCache)
    @assert isempty_known(c) "'isempty_known' only works if 'isempty' " *
        "returns 'true'"
    return c.isempty == Int8(1)
end

function set_isempty!(c::IntersectionCache, isempty::Bool)
    c.isempty = isempty ? Int8(1) : Int8(0)
end

"""
    Intersection{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}

Type that represents the intersection of two convex sets.

### Fields

- `X`     -- convex set
- `Y`     -- convex set
- `cache` -- internal cache for avoiding recomputation; see
             [`IntersectionCache`](@ref)

### Examples

Create an expression, ``Z``, which lazily represents the intersection of two
squares ``X`` and ``Y``:

```jldoctest lazy_intersection
julia> X, Y = BallInf([0,0.], 0.5), BallInf([1,0.], 0.75);

julia> Z = X ∩ Y;

julia> typeof(Z)
Intersection{Float64,BallInf{Float64,Array{Float64,1}},BallInf{Float64,Array{Float64,1}}}

julia> dim(Z)
2
```

We can check if the intersection is empty with `isempty`:

```jldoctest lazy_intersection
julia> isempty(Z)
false
```

Do not confuse `Intersection` with the concrete operation, which is computed
with the lowercase `intersection` function:

```jldoctest lazy_intersection
julia> W = intersection(X, Y)
Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}([0.375, 0.0], [0.125, 0.5])
```
"""
struct Intersection{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2
    cache::IntersectionCache

    # default constructor with dimension check
    function Intersection(X::S1, Y::S2;
                          cache::IntersectionCache=IntersectionCache()
                         ) where {N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}
        @assert dim(X) == dim(Y) "sets in an intersection must have the same " *
            "dimension"
        return new{N, S1, S2}(X, Y, cache)
    end
end

isoperationtype(::Type{<:Intersection}) = true
isconvextype(::Type{Intersection{N, S1, S2}}) where {N, S1, S2} = isconvextype(S1) && isconvextype(S2)

# Universe is the neutral element for Intersection
@neutral(Intersection, Universe)

# EmptySet is the absorbing element for Intersection
@absorbing(Intersection, EmptySet)

"""
    ∩

Alias for `Intersection`.
"""
∩(X::LazySet, Y::LazySet) = Intersection(X, Y)

# --- cache propagation functions ---

"""
    isempty_known(cap::Intersection)

Ask whether the status of emptiness is known.

### Input

- `cap` -- intersection of two convex sets

### Output

`true` iff the emptiness status is known.
In this case, `isempty(cap)` can be used to obtain the status.
"""
function isempty_known(cap::Intersection)
    return isempty_known(cap.cache)
end

"""
    set_isempty!(cap::Intersection, isempty::Bool)

Set the status of emptiness in the cache.

### Input

- `cap`     -- intersection of two convex sets
- `isempty` -- new status of emptiness
"""
function set_isempty!(cap::Intersection, isempty::Bool)
    return set_isempty!(cap.cache, isempty)
end

"""
    swap(cap::Intersection{N, S1, S2}) where {N<:Real, S1, S2}

Return a new `Intersection` object with the arguments swapped.

### Input

- `cap` -- intersection of two convex sets

### Output

A new `Intersection` object with the arguments swapped.
The old cache is shared between the old and new objects.

### Notes

The advantage of using this function instead of manually swapping the arguments
is that the cache is shared.
"""
function swap(cap::Intersection{N, S1, S2}) where {N<:Real, S1, S2}
    return Intersection(cap.Y, cap.X, cache=cap.cache)
end

# --- LazySet interface functions ---

"""
    dim(cap::Intersection)

Return the dimension of an intersection of two convex sets.

### Input

- `cap` -- intersection of two convex sets

### Output

The ambient dimension of the intersection of two convex sets.
"""
function dim(cap::Intersection)
    return dim(cap.X)
end

"""
    σ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}

Return the support vector of an intersection of two convex sets in a given
direction.

### Input

- `d`   -- direction
- `cap` -- intersection of two convex sets

### Output

The support vector in the given direction.
"""
function σ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}
    error("the exact support vector of an intersection is not implemented")
end

"""
    ρ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}

Return an upper bound on the support function of the intersection of two convex
sets in a given direction.

### Input

- `d`    -- direction
- `cap`  -- intersection of two convex sets

### Output

An uper bound on the support function in the given direction.

### Algorithm

The support function of an intersection of ``X`` and ``Y`` is upper bounded by
the minimum of the support functions of ``X`` and ``Y``.
"""
function ρ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}
    return min(ρ(d, cap.X), ρ(d, cap.Y))
end

function ρ_helper(d::AbstractVector{N},
                  cap::Intersection{N,
                                    <:LazySet{N},
                                    <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}},
                  algorithm::String;
                  kwargs...) where {N<:Real}
    @assert isbounded(cap.X) "the first set in the intersection must be bounded"
    X = cap.X # compact set
    H = cap.Y # halfspace or hyperplane or line

    # if the intersection is empty => stop
    if isempty(cap)
        error("the intersection is empty")
    end

    if !use_precise_ρ(cap) || algorithm == "simple"
        return invoke(ρ, Tuple{typeof(d), Intersection{N}}, d, cap)
    elseif algorithm == "line_search"
        require(:Optim; fun_name="ρ", explanation="(algorithm $algorithm)")
        (s, _) = _line_search(d, X, H; kwargs...)
        return s
    elseif algorithm == "projection"
        @assert H isa Hyperplane "the algorithm $algorithm cannot be used " *
            "with a $(typeof(H)); it only works with hyperplanes"
        return _projection(d, X, H; kwargs...)
    else
        error("algorithm $(algorithm) unknown")
    end
end

"""
    use_precise_ρ(cap::Intersection{N}) where {N<:Real}

Determine whether a precise algorithm for computing ``ρ`` shall be applied.

### Input

- `cap` -- intersection of two convex sets

### Output

`true` if a precise algorithm shall be applied.

### Notes

The default implementation always returns `true`.

If the result is `false`, a coarse approximation of the support function is
returned.

This function can be overwritten by the user to control the policy.
"""
function use_precise_ρ(cap::Intersection{N}) where {N<:Real}
    return true
end

"""
    ρ(d::AbstractVector{N},
      cap::Intersection{N, S1, S2};
      [algorithm]::String="line_search",
      [kwargs...]) where {N<:Real,
                          S1<:LazySet{N},
                          S2<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}

Return the support function of the intersection of a compact set and a
half-space/hyperplane/line in a given direction.

### Input

- `d`         -- direction
- `cap`       -- lazy intersection of a compact set and a half-space/hyperplane/
                 line
- `algorithm` -- (optional, default: `"line_search"`): the algorithm to
                 calculate the support function; valid options are:

    * `"line_search"` -- solve the associated univariate optimization problem
                         using a line search method (either Brent or the
                         Golden Section method)
    * `"projection"`  -- only valid for intersection with a hyperplane;
                         evaluates the support function by reducing the problem
                         to the 2D intersection of a rank 2 linear
                         transformation of the given compact set in the plane
                         generated by the given direction `d` and the
                         hyperplane's normal vector `n`
    * `"simple"`      -- take the ``\\min`` of the support function evaluation
                         of each operand

### Output

The scalar value of the support function of the set `cap` in the given
direction.

### Notes

It is assumed that the set `cap.X` is compact.

Any additional number of arguments to the algorithm backend can be passed as
keyword arguments.

### Algorithm

The algorithms are based on solving the associated optimization problem

```math
\\min_{λ ∈ D_h} ρ(ℓ - λa, X) + λb.
```
where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or
``D_h = \\{ λ : λ ∈ \\mathbb{R} \\}`` if ``H`` is a hyperplane.

For additional information we refer to:

- [G. Frehse, R. Ray. Flowpipe-Guard Intersection for Reachability Computations with
  Support Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809).
- [C. Le Guernic. Reachability Analysis of Hybrid Systems with Linear Continuous
  Dynamics, PhD thesis](https://tel.archives-ouvertes.fr/tel-00422569v2).
- [T. Rockafellar, R. Wets.
  Variational Analysis](https://www.springer.com/us/book/9783540627722).
"""
function ρ(d::AbstractVector{N},
           cap::Intersection{N, S1, S2};
           algorithm::String="line_search",
           kwargs...) where {N<:Real,
                             S1<:LazySet{N},
                             S2<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}
    return ρ_helper(d, cap, algorithm; kwargs...)
end

# symmetric method
function ρ(d::AbstractVector{N},
           cap::Intersection{N, S1, S2};
           algorithm::String="line_search",
           kwargs...) where {N<:Real,
                             S1<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
                             S2<:LazySet{N}}
    return ρ_helper(d, swap(cap), algorithm; kwargs...)
end

"""
    ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2}; kwargs...)
        where {N<:Real, S1<:LazySet{N}, S2<:AbstractPolyhedron{N}}

Return an upper bound on the support function of the intersection between a
compact set and a polyhedron along a given direction.

### Input

- `d`      -- direction
- `cap`    -- intersection of a compact set and a polyhedron
- `kwargs` -- additional arguments that are passed to the support-function
              algorithm

### Output

An upper bound of the support function of the given intersection.

### Algorithm

The idea is to solve the univariate optimization problem `ρ(di, X ∩ Hi)` for each
half-space in the set `P` and then take the minimum. This gives an overapproximation
of the exact support function.

This algorithm is inspired from [G. Frehse, R. Ray. Flowpipe-Guard Intersection
for Reachability Computations with Support
Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809).

### Notes

This method relies on the `constraints_list` of the polyhedron.
"""
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2}; kwargs...
          ) where {N<:Real, S1<:LazySet{N}, S2<:AbstractPolyhedron{N}}
    return ρ_helper(d, cap; kwargs...)
end

# symmetric method
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2}; kwargs...
          ) where {N<:Real, S1<:AbstractPolyhedron{N}, S2<:LazySet{N}}
    return ρ_helper(d, swap(cap); kwargs...)
end

# disambiguation
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2}; kwargs...
          ) where {N<:Real, S1<:AbstractPolytope{N}, S2<:AbstractPolyhedron{N}}
    return ρ_helper(d, cap; kwargs...)
end

function ρ_helper(d::AbstractVector{N}, cap::Intersection{N, S1, S2}; kwargs...
                 ) where {N<:Real, S1<:LazySet{N}, S2<:AbstractPolyhedron{N}}
    if !use_precise_ρ(cap)
        use_simple_method = true
    else
        options = Dict(kwargs)
        use_simple_method = haskey(options, :algorithm) &&
                       options[:algorithm] == "simple"
    end

    if use_simple_method
        # simple algorithm
        return invoke(ρ, Tuple{typeof(d), Intersection{N}}, d, cap)
    end

    # more precise algorithm
    @assert isbounded(cap.X) "the first set in the intersection must be bounded"
    return minimum([ρ(d, cap.X ∩ Hi; kwargs...)
                   for Hi in constraints_list(cap.Y)])
end

"""
    ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2}; kwargs...
     ) where {N<:Real, S1<:AbstractPolyhedron{N}, S2<:AbstractPolyhedron{N}}

Return an upper bound on the support function of the intersection between two
polyhedral sets.

### Input

- `d`      -- direction
- `cap`    -- intersection of two polyhedral sets
- `kwargs` -- additional arguments that are passed to the support-function
              algorithm

### Output

The support function for the given direction.

### Algorithm

We combine the constraints of the two polyhedra to a new `HPolyhedron`, for
which we then evaluate the support function.
"""
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2}; kwargs...
          ) where {N<:Real, S1<:AbstractPolyhedron{N}, S2<:AbstractPolyhedron{N}}
    return ρ(d, HPolyhedron([constraints_list(cap.X); constraints_list(cap.Y)]))
end

# disambiguation
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2};
           algorithm::String="line_search", kwargs...
          ) where {N<:Real, S1<:AbstractPolytope{N},
                   S2<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}
    return ρ_helper(d, cap, algorithm; kwargs...)
end

# symmetric method
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2};
           algorithm::String="line_search", kwargs...
          ) where {N<:Real, S1<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
                   S2<:AbstractPolytope{N}}
    return ρ_helper(d, swap(cap), algorithm; kwargs...)
end

# disambiguation
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2};
           algorithm::String="line_search", kwargs...
          ) where {N<:Real, S1<:AbstractPolyhedron{N},
                   S2<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}
    return ρ(d, HPolyhedron([constraints_list(cap.X); constraints_list(cap.Y)]))
end

# symmetric method
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2};
           algorithm::String="line_search", kwargs...
          ) where {N<:Real, S1<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
                   S2<:AbstractPolyhedron{N}}
    return ρ(d, HPolyhedron([constraints_list(cap.X); constraints_list(cap.Y)]))
end

# disambiguation
function ρ(d::AbstractVector{N}, cap::Intersection{N, S1, S2}
          ) where {N<:Real, S1<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
                   S2<:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}
    return ρ(d, HPolyhedron([constraints_list(cap.X); constraints_list(cap.Y)]))
end

"""
    isbounded(cap::Intersection)

Determine whether an intersection of two convex sets is bounded.

### Input

- `cap` -- intersection of two convex sets

### Output

`true` iff the intersection is bounded.

### Algorithm

We first check if any of the wrapped sets is bounded.
Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref).
"""
function isbounded(cap::Intersection)
    if isbounded(cap.X) || isbounded(cap.Y)
        return true
    end
    return isbounded_unit_dimensions(cap)
end

"""
    ∈(x::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}

Check whether a given point is contained in an intersection of two convex sets.

### Input

- `x`   -- point/vector
- `cap` -- intersection of two convex sets

### Output

`true` iff ``x ∈ cap``.
"""
function ∈(x::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}
    return (x ∈ cap.X) && (x ∈ cap.Y)
end

"""
    constraints_list(cap::Intersection{N}) where {N<:Real}

Return the list of constraints of an intersection of two (polyhedral) sets.

### Input

- `cap` -- intersection of two (polyhedral) sets

### Output

The list of constraints of the intersection.

### Notes

We assume that the underlying sets are polyhedral, i.e., offer a method
`constraints_list`.

### Algorithm

We create the polyhedron by taking the intersection of the `constraints_list`s of
the sets and remove redundant constraints.

This function ignores the boolean output from the in-place `remove_redundant_constraints!`,
which may inform the user that the constraints are infeasible. In that case, the
list of constraints at the moment when the infeasibility was detected is returned.
"""
function constraints_list(cap::Intersection{N}) where {N<:Real}
    constraints = [constraints_list(cap.X); constraints_list(cap.Y)]
    remove_redundant_constraints!(constraints)
    return constraints
end

# --- Intersection functions ---

"""
    isempty(cap::Intersection)

Return if the intersection is empty or not.

### Input

- `cap` -- intersection of two convex sets

### Output

`true` iff the intersection is empty.

### Notes

The result will be cached, so a second query will be fast.
"""
function isempty(cap::Intersection)
    if isempty_known(cap)
        # use cached result
        return isempty(cap.cache)
    end
    # compute result
    empty_intersection = is_intersection_empty(cap.X, cap.Y)
    # update cache
    set_isempty!(cap, empty_intersection)

    return empty_intersection
end

"""
    plot_recipe(cap::Intersection{N}, [ε]::N=-one(N),
                [Nφ]::Int=PLOT_POLAR_DIRECTIONS) where {N<:Real}

Convert a lazy intersection to a pair `(x, y)` of points for plotting.

### Input

- `cap`  -- lazy intersection
- `ε`    -- (optional, default `0`) ignored, used for dispatch
- `Nφ`   -- (optional, default: `PLOT_POLAR_DIRECTIONS`) number of polar
            directions used in the template overapproximation

### Output

A pair `(x, y)` of points that can be plotted.
"""
function plot_recipe(cap::Intersection{N}, ε::N=zero(N),
                     Nφ::Int=PLOT_POLAR_DIRECTIONS) where {N<:Real}
    @assert dim(cap) <= 2 "cannot plot a $(dim(cap))-dimensional intersection"

    if isempty(cap)
        return plot_recipe(EmptySet{N}(dim(cap)), ε)
    elseif dim(cap) == 1
        return plot_recipe(convert(Interval, cap), ε)
    else
        # construct polygon approximation using polar directions
        P = overapproximate(cap, PolarDirections{N}(Nφ))
        return plot_recipe(P, ε)
    end
end

# ==========================================================
# Concrete operations that dispatch on a lazy Intersection
# ==========================================================

"""
    linear_map(M::AbstractMatrix{N}, cap::Intersection{N}) where {N}

Return the concrete linear map of a lazy intersection.

### Input

- `M`   -- matrix
- `cap` -- lazy intersection

### Output

The set obtained by applying the given linear map to the lazy intersection.

### Notes

This function relies on computing `cap` concretely (i.e. as a set representation),
and then applying the linear map.
"""
function linear_map(M::AbstractMatrix{N}, cap::Intersection{N}) where {N}
    return linear_map(M, intersection(cap.X, cap.Y))
end

# ==================================
# Algorithms for lazy intersection
# ==================================

function load_optim_intersection()
return quote

"""
    _line_search(ℓ, X, H::Union{<:HalfSpace, <:Hyperplane, <:Line}; [kwargs...])

Given a compact and convex set ``X`` and a halfspace ``H = \\{x: a^T x ≤ b \\}``
or a hyperplane ``H = \\{x: a^T x = b \\}``, calculate:

```math
\\min_{λ ∈ D_h} ρ(ℓ - λa, X) + λb.
```
where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or
``D_h = \\{ λ : λ ∈ \\mathbb{R} \\}`` if ``H`` is a hyperplane.

### Input

- `ℓ`           -- direction
- `X`           -- set
- `H`           -- halfspace or hyperplane

### Output

The tuple `(fmin, λmin)`, where `fmin` is the minimum value of the function
``f(λ) = ρ(ℓ - λa) + λb`` over the feasible set ``λ ≥ 0``, and ``λmin`` is the
minimizer.

### Notes

This function requires the `Optim` package, and relies on the univariate
optimization interface `Optim.optimize(...)`.

Additional arguments to the `optimize` backend can be passed as keyword arguments.
The default method is `Optim.Brent()`.

### Examples

```jldoctest _line_search
julia> X = Ball1(zeros(2), 1.0);

julia> H = HalfSpace([-1.0, 0.0], -1.0);  # x >= 1

julia> using Optim

julia> using LazySets: _line_search

julia> v = _line_search([1.0, 0.0], X, H);  # uses Brent's method by default

julia> v[1]
1.0
```

We can specify the upper bound in Brent's method:

```julia _line_search
julia> v = _line_search([1.0, 0.0], X, H, upper=1e3);

julia> v[1]
1.0
```

Instead of Brent's method we can use the Golden Section method:

```julia _line_search
julia> v = _line_search([1.0, 0.0], X, H, upper=1e3, method=GoldenSection());

julia> v[1]
1.0
```
"""
function _line_search(ℓ, X, H::Union{<:HalfSpace, <:Hyperplane, <:Line};
                      kwargs...)
    options = Dict(kwargs)

    # Initialization
    a, b = H.a, H.b
    f(λ) = ρ(ℓ - λ[1] * a, X) + λ[1] * b

    if haskey(options, :lower)
        lower = pop!(options, :lower)
    else
        if H isa HalfSpace
            lower = 0.0
        elseif (H isa Hyperplane) || (H isa Line)
            lower = -1e6 # "big": TODO relate with f(λ)
        end
    end

    if haskey(options, :upper)
        upper = pop!(options, :upper)
    else
        upper = 1e6 # "big": TODO relate with f(λ)
    end

    if haskey(options, :method)
        method = pop!(options, :method)
    else
        method = Optim.Brent()
    end

    # Optimization
    sol = Optim.optimize(f, lower, upper, method=method, options...)

    # Recover results
    fmin, λmin = sol.minimum, sol.minimizer
    return (fmin, λmin)
end # _line_search
end # quote
end # load_optim

"""
    _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}};
                [lazy_linear_map]=false,
                [lazy_2d_intersection]=true,
                [algorithm_2d_intersection]=nothing,
                [kwargs...]) where {N}

Given a compact and convex set ``X`` and a hyperplane ``H = \\{x: n ⋅ x = γ \\}``,
calculate the support function of the intersection between the
rank-2 projection ``Π_{nℓ} X`` and the line ``Lγ = \\{(x, y): x = γ \\}``.

### Input

- `ℓ`                    -- direction
- `X`                    -- set
- `H`                    -- hyperplane
- `lazy_linear_map`      -- (optional, default: `false`) to perform the projection
                            lazily or concretely
- `lazy_2d_intersection` -- (optional, default: `true`) to perform the 2D
                            intersection between the projected set and the line
                            lazily or concretely
- `algorithm_2d_intersection` -- (optional, default: `nothing`) if given, fixes the
                                 support function algorithm used for the intersection
                                 in 2D; otherwise the default is implied

### Output

The support function of ``X ∩ H`` along direction ``ℓ``.

### Algorithm

This projection method is based on Prop. 8.2, page 103, [C. Le Guernic.
Reachability Analysis of Hybrid Systems with Linear Continuous Dynamics,
PhD thesis](https://tel.archives-ouvertes.fr/tel-00422569v2).

In the original algorithm, Section 8.2 of Le Guernic's thesis, the linear map
is performed concretely and the intersection is performed lazily (these are the
default options in this algorithm, but here the four combinations are available).
If the set ``X`` is a zonotope, its concrete projection is again a zonotope
(sometimes called "zonogon"). The intersection between this zonogon and the line
can be taken efficiently in a lazy way (see Section 8.2.2 of Le Guernic's thesis),
if one uses dispatch on `ρ(y_dir, Sℓ⋂Lγ; kwargs...)` given that `Sℓ` is itself
a zonotope.

### Notes

This function depends itself on the calculation of the support function of another
set in two dimensions. Obviously one doesn't want to use again `algorithm="projection"`
for this second calculation. The option `algorithm_2d_intersection` is such that,
if it is not given, the default support function algorithm is used (e.g. `"line_search"`).
You can still pass additional arguments to the `"line_search"` backend through the
`kwargs`.
"""
function _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}};
                     lazy_linear_map=false,
                     lazy_2d_intersection=true,
                     algorithm_2d_intersection=nothing,
                     kwargs...) where {N}

    n = H.a                  # normal vector to the hyperplane
    γ = H.b                  # displacement of the hyperplane
    Πnℓ = vcat(n', ℓ')       # projection map

    x_dir = [one(N), zero(N)]
    y_dir = [zero(N), one(N)]

    Lγ = Line(x_dir, γ)

    Xnℓ = lazy_linear_map ? LinearMap(Πnℓ, X) : linear_map(Πnℓ, X)
    Xnℓ⋂Lγ = lazy_2d_intersection ? Intersection(Xnℓ, Lγ) : intersection(Xnℓ, Lγ)

    if algorithm_2d_intersection == nothing
        return ρ(y_dir, Xnℓ⋂Lγ; kwargs...)
    else
        return ρ(y_dir, Xnℓ⋂Lγ, algorithm=algorithm_2d_intersection; kwargs...)
    end
end

"""
    get_constrained_lowdimset(cpa::CartesianProductArray{N, S},
                              P::AbstractPolyhedron{N}
                             ) where {N<:Real, S<:LazySet{N}}

Preprocess step for intersection between Cartesian product array and polyhedron.
Returns low-dimensional a `CartesianProductArray` in the constrained dimensions
of the original cpa,
constrained variables and variables in corresponding blocks, original block
structure of low-dimensional set and list of constrained blocks.

### Input

- `cpa` -- Cartesian product array of convex sets
- `P`   -- polyhedron

### Output

A tuple of low-dimensional set, list of constrained dimensions, original block
structure of low-dimensional set and corresponding blocks indices.
"""
function get_constrained_lowdimset(cpa::CartesianProductArray{N, S},
                                   P::AbstractPolyhedron{N}
                                  ) where {N<:Real, S<:LazySet{N}}

    if isbounded(P)
        blocks, non_empty_length = block_to_dimension_indices(cpa)
    else
        blocks, non_empty_length =
            block_to_dimension_indices(cpa, constrained_dimensions(P))
    end

    array = Vector{S}()
    sizehint!(array, non_empty_length)
    variables = Vector{Int}()
    block_structure = Vector{UnitRange{Int}}()
    sizehint!(block_structure, non_empty_length)

    last_var = 1
    for i in 1:length(blocks)
        start_index, end_index = blocks[i]
        block_end = last_var + end_index - start_index
        if start_index != -1
            push!(array, cpa.array[i])
            append!(variables, start_index : end_index)
            push!(block_structure, last_var : block_end)
            last_var = block_end + 1
        end
    end

    return CartesianProductArray(array), variables, block_structure, blocks
end


text:
# This file contains functions for creating spatial weights from points data using a distane threshold

"""
    dnearneigh(X, Y; threshold)
Build a spatial weights object from a set of coordinates using a distance `threshold`.
"""
function dnearneigh(X::Vector{Float64}, Y::Vector{Float64}; threshold::Real)::SpatialWeights
    n, ny = length(X), length(Y)

    n == ny || throw(DimensionMismatch("dimensions must match: X has length ($(n)), Y has length ($ny)"))

    # Get points within range using a KDTree
    cpoints = vcat(X', Y')

    kdtree = KDTree(cpoints)
    idxs  = inrange(kdtree, cpoints, threshold, true)

    # Check neighbours
    neighs = copy.(fill(Int[], n))
    weights = copy.(fill(Float64[], n))
    nneighs = zeros(Int,n)

    for i in 1:n
        neighsi = idxs[i]
        neighs[i] = sort!(neighsi[neighsi .!= i])
        nneighs[i] = length(neighs[i])
        weights[i] = ones(nneighs[i]) ./ nneighs[i]
    end

    SpatialWeights(n, neighs, weights, nneighs, :row)

end

"""
    dnearneigh(P; threshold)
Build a spatial weights object from a vector of points using a distance `threshold`.
"""
function dnearneigh(P::Vector{T} where T<:Union{Missing,AbstractPoint}; threshold::Real)::SpatialWeights
    cpoints = reduce(hcat, map(a -> coordinates(a), P))

    dnearneigh(cpoints[1,:], cpoints[2,:]; threshold)
end

"""
    dnearneigh(A; threshold)
Build a spatial weights object from a table A that constains a points geometry column using a distance `threshold`.
"""
function dnearneigh(A::Any; threshold::Real)::SpatialWeights
    istable(A) || throw(ArgumentError("Unknown geometry or not points geometry"))

    (:geometry in propertynames(A)) || throw(ArgumentError("table does not have :geometry information"))

    dnearneigh(A.geometry; threshold)
end


text:
"""
Represents the processed spectral data necessary to efficiently filter-fit one or more unknown spectra.
A `FilterFitPacket` contains the data necessary to filter the unknown and to apply pre-filtered references.
If there are duplicate `FilteredReference` for an elemental ROI, the preference is for the first one.  This
allows you to fill in unavailable "FilteredReference" elemental ROIs with more general ones.
"""
struct FilterFitPacket{T<:Detector}
    detector::T
    filter::TopHatFilter
    references::Vector{FilteredReference}

    function FilterFitPacket(det::T, filt::TopHatFilter, refs::Vector{FilteredReference}) where {T <: Detector}
        # Only permit one FilteredReference for each type of label for a set of X-rays
        filtrefs = FilteredReference[]
        for ref in refs
            exists = false
            for (i, fr) in enumerate(filtrefs)
                # If a reference already exists of this type for these X-rays reassign it
                if (typeof(fr.label) == typeof(ref.label))  && (fr.label.xrays==ref.label.xrays)
                    # filtrefs[i] = ref
                    exists = true
                    break
                end
            end
            if !exists
                push!(filtrefs, ref)
            end
        end
        new{T}(det, filt, filtrefs)
    end
end # struct

function Base.show(io::IO, ffp::FilterFitPacket)
    xrays = join(map(r -> "\t$(r.label),", ffp.references), "\n")
    print(io, "References[\n\t$(ffp.detector), \n$xrays\n]")
end

function NeXLUncertainties.asa(::Type{DataFrame}, ffp::FilterFitPacket)
    function p2b(fref)
        croi = first(fref.roi)-first(fref.ffroi):last(fref.roi)-first(fref.ffroi)
        sum(fref.charonly) / (sum(fref.data[croi])-sum(fref.charonly))
    end
    function s2n(fref)
        croi = first(fref.roi)-first(fref.ffroi):last(fref.roi)-first(fref.ffroi)
        sum(fref.charonly) / sqrt(sum(fref.data[croi])-sum(fref.charonly))
    end
    DataFrame(
        :Lines => [ fr.label.xrays for fr in ffp.references],
        :Material => [ get(fr.label.spectrum, :Composition, nothing) for fr in ffp.references],
        :ROI => [ fr.roi for fr in ffp.references],
        # :FullROI => [ fr.ffroi for fr in ffp.references],
        Symbol("P-to-B") => p2b.(ffp.references),
        Symbol("S-to-N") => s2n.(ffp.references)
    )
end

"""
    missingReferences(ffp::FilterFitPacket, elms::Vector{Element}, e0::Float64, ampl=1.0e-5)

Returns a `Vector{Tuple{Vector{CharXRay}, UnitRange{Int64}}}` containing the ROIs for which a
`FilteredReference` is missing from the `FilterFitPacket`.
"""
function missingReferences(ffp::FilterFitPacket, elms::Vector{Element}, e0::Float64, ampl=1.0e-5)
    # Find all collections of X-rays for which we'll need references
    allXrays = mapreduce(elm->NeXLSpectrum.labeledextents(elm, ffp.detector, ampl, e0), append!, elms)
    # Remove those for which there are references
    filter(allXrays) do xr
        isnothing(findfirst(ref->ref.label.xrays == xr[1], ffp.references))
    end
end

"""
Returns a tuple containing the unfiltered spectra associated with the references.
"""
spectra(ffp::FilterFitPacket) = (ref.label.spectrum for ref in ffp.references)

"""
A list of elements for which there are filtered references.
"""
NeXLCore.elms(ffp::FilterFitPacket) = union( [ element(ref.label) for ref in filter(r->r.label isa CharXRayLabel, ffp.references) ] )

struct ReferencePacket
    spectrum::Spectrum
    element::Element
    material::Material
end

Base.show(io::IO, rp::ReferencePacket) = print(io, "ReferencePacket[$(symbol(rp.element)), $(name(rp.material)), $(name(rp.spectrum))]")

"""
    reference( elm::Element, spec::Spectrum, mat::Material=spec[:Composition]; pc = nothing, lt = nothing, e0 = nothing, coating = nothing)::ReferencePacket
    reference(els::AbstractVector{Element}, spec::Spectrum, mat::Material = spec[:Composition]; pc = nothing, lt = nothing, e0 = nothing, coating = nothing)::Vector{ReferencePacket}

Construct a `ReferencePacket` from a `Spectrum` collected from the specified `Material` for the specified `Element`.
Often used with `references(...)` to build `FilterFitPacket`s.

Optional named arguments `pc`, `lt`, `e0`, `coating` allow you to specify the probe current, live time, beam energy and
sample coating.
"""
function reference(
    elm::Element,
    spec::Spectrum,
    mat::Material = spec[:Composition];
    pc = nothing,
    lt = nothing,
    e0 = nothing,
    coating = nothing,
)::ReferencePacket
    if !isnothing(lt)
        @assert lt > 0.0 "The live time must be larger than zero for $(spec[:Name])."
        spec[:LiveTime] = lt
    end
    if !isnothing(pc)
        @assert pc > 0.0 "The probe current must be larger than zero for $(spec[:Name])."
        spec[:ProbeCurrent] = pc
    end
    if !isnothing(e0)
        @assert e0 >= 1.0e3 "The beam energy ($(e0/1000.0)) keV is too low for $(spec[:Name])."
        spec[:BeamEnergy] = e0
    end
    if !isnothing(coating)
        @assert coating isa Film
        spec[:Coating] = coating
    end
    spec[:Composition] = mat
    @assert haskey(spec, :LiveTime) "The :LiveTime property must be defined for $(spec[:Name]).  (Use the `lt` keyword argument)"
    @assert haskey(spec, :ProbeCurrent) "The :ProbeCurrent property must be defined for $(spec[:Name]).  (Use the `pc` keyword argument)"
    @assert haskey(spec, :BeamEnergy) "The :BeamEnergy property must be defined for $(spec[:Name]).  (Use the `e0` keyword argument)"
    @assert haskey(mat, elm) "$(Symbol(elm)) is not present in $mat."
    return ReferencePacket(spec, elm, mat)
end

reference(elm::Element, filename::AbstractString, mat::Material; kwargs...) =
    reference(elm, loadspectrum(filename), mat; kwargs...)

reference(elm::Element, filename::AbstractString; kwargs...) =
    reference(elm, loadspectrum(filename); kwargs...)

function reference(
    els::AbstractVector{Element},
    spec::Spectrum,
    mat::Material = spec[:Composition];
    kwargs...)
    map(el->reference(el, spec, mat; kwargs...), els)
end
reference(elm::AbstractVector{Element}, filename::AbstractString, mat::Material; kwargs...) =
    reference(elm, loadspectrum(filename), mat; kwargs...)

"""
    references(refs::AbstractVector{ReferencePacket}, det::EDSDetector)::FilterFitPacket
    references(refs::AbstractVector{ReferencePacket}, fwhm::Float64)::FilterFitPacket

Constructs a FilterFitPacket from a vector of `ReferencePackets`.  Each `ReferencePacket` represents a
single ROI for an element.  It is possible more than one `ReferencePacket` might be defined for an
elemental ROI.  In this case, the `ReferencePacket` with the lower index will take preference over
later ones.  This allows you to fill in only the missing elemental ROIs using spectra collected from
alternative materials.  For example, a spectrum from F₂Fe is suitable for the Fe K-lines but not the
Fe L-lines. So we might specify F₂Fe first to specify the references for the Fe K-lines and then fill
in the L-lines with a spectrum from pure Fe.
"""
function references(
    refs::AbstractVector{ReferencePacket},
    det::EDSDetector,
)::FilterFitPacket
    chcount = det.channelcount
    @assert all(length(r.spectrum) == chcount for r in refs) "The number of spectrum channels must match the detector for all spectra."
    @assert length(refs) > 0 "Please provide at least one ReferencePacket in references(...)"
    # Build the top-hat filter for det
    ff = buildfilter(det)
    # Apply the top-hat filter to all refs
    frefs = mapreduce(append!, refs) do ref
        frefs = filterreference(ff, ref.spectrum, ref.element, ref.material)
        length(frefs)==0 && @warn "Unable to create any filtered ROI references for $(ref.element) from $(name(ref.material))."
        frefs
    end
    return FilterFitPacket(det, ff, frefs)
end
references(refs::AbstractVector{ReferencePacket}, fwhm::Float64)::FilterFitPacket =
    references(refs, matching(refs[1].spectrum, fwhm))

fit_spectrum(spec::Spectrum, ffp::FilterFitPacket) =
    fit_spectrum(FilteredUnknownW, spec, ffp.filter, ffp.references)

fit_spectrum(specs::AbstractVector{<:Spectrum}, ffp::FilterFitPacket) =
    map(spec -> fit_spectrum(FilteredUnknownW, spec, ffp.filter, ffp.references), specs)

"""
fit_spectrum(hs::HyperSpectrum, ffp::FilterFitPacket; mode::Symbol=:Fast, zero = x -> max(0.0, x))::Array{KRatios}

  * mode = :Fast - Uses precomputed, filtered "vector" fit method.  No uncertainties are available.
  * mode = :Intermediate - Uses an optimized full fit without refits for negative k-ratios.
  * mode = :Full - Uses the full single spectrum fit algorithm including refitting when one or more k-ratio is less than zero.

Performs a filtered fit on a hyperspectrum returning an `Array{KRatios}`.

Selecting a mode:
  :Fast is good for generating k-ratio maps or exploratory analysis of a k-ratio map. :Full is best when a
  quantitative map of a high count hyperspectrum is desired.  Fit frefsults for major elements are similar for
  all three but differ for minor and trace elements.  Particularly when a k-ratio is slightly negative. This
  negative k-ratio can effect the other k-ratios.  :Fast also works less well when many elements (>>10) (particularly
  interfering elements) are included in the fit. Unfortunately, :Intermediate and :Full slow down when many elements
  are fit - O(n(elements)²).

The following timing on a 512 x 512 x 2048 hyperspectrum fitting 15 elements with 25 distinct ROIs on a fast laptop
with 64 GiB memory give a relative feel for the speed of each algorithm.  Yes, :Fast is approximately 20x faster than
:Intermediate and used almost 100x less memory.  (Single thread timings)

|---------------|----------------|--------------|--------------|----------|
| mode          | Run time (s)   | Allocations  | Memory (GiB) | GC time  |
|---------------|----------------|--------------|--------------|----------|
| :Fast         | 44.8           | 3.95 M       | 5.72         | 4.2%     |
| :Intermediate | 1305.1         | 24.83 M      | 523.7        | 2.7%     |
| :Full         | 2056           | 2.7 G        | 786.1        | 4.2%     |
|---------------|----------------|--------------|--------------|----------|
"""
function fit_spectrum(
    hs::HyperSpectrum,
    ffp::FilterFitPacket;
    mode::Symbol = :Fast,
    zero = x -> max(0.0, x),
    sigma = 0.0
)::Array{KRatios}
    unklabel = UnknownLabel(hs)
     function _tophatfilterhs(unklabel, data, thf, scale)
        @assert length(data) <= length(thf) "The reference spectra have fewer channels than the hyperspectrum data."
        filtered = Float64[filtereddatum(thf, data, i) for i in eachindex(data)]
        dp = Float64[max(x, 1.0) for x in data] # To ensure covariance isn't zero or infinite precision
        covar = Float64[filteredcovar(thf, dp, i, i) for i in eachindex(data)]
        FilteredUnknownW(unklabel, scale, eachindex(data), data, filtered, covar)
    end
    function fitcontiguousx(unk, ffs, chs)
        _buildscale(unk, ffs) * pinv(_buildmodel(ffs, chs), rtol = 1.0e-6) * extract(unk, chs)
    end
    function _filterfitx(unk, ffs, fitrois)
        mapreduce(append!, fitrois) do fr
            fitcontiguousx(
                unk,
                filter(ff -> length(intersect(fr, ff.ffroi)) > 0, ffs),
                fr,
            )
        end
    end
    @assert matches(hs[CartesianIndices(hs)[1]], ffp.detector) "The detector for the hyper-spectrum must match the detector for the filtered references."
    if mode == :Fast
        vq = VectorQuant(ffp.references, ffp.filter)
        return fit_spectrum(vq, hs, zero)
    elseif mode == :Intermediate
        krs = zeros(Float32, length(ffp.references), size(hs)...)
        len = 1:depth(hs)
        data = zeros(Float64, length(ffp.filter))
        fitrois = ascontiguous(map(fd -> fd.ffroi, ffp.references))
        ThreadsX.foreach(CartesianIndices(hs)) do ci
            data[len] .= view(counts(hs), len, ci)
            unk = _tophatfilterhs(unklabel, data, ffp.filter, 1.0/dose(hs,ci))
            krs[:, ci] .= zero.(_filterfitx(unk, ffp.references, fitrois))
        end
        return ThreadsX.map(filter(ii -> ffp.references[ii].label isa CharXRayLabel, eachindex(ffp.references))) do i
            k, lbl = krs[i], ffp.references[i].label
            rprops = properties(spectrum(lbl))
            KRatios(xrays(lbl), properties(hs), rprops, rprops[:Composition], krs[i, :, :])
        end
    elseif mode == :Full
        krs = zeros(Float32, length(ffp.references), size(hs)...)
        len = 1:depth(hs)
        data = zeros(Float64, length(ffp.filter))
        ThreadsX.foreach(CartesianIndices(hs)) do ci
            data[len] .= view(counts(hs), len, ci)
            unk = _tophatfilterhs(unklabel, data, ffp.filter, 1.0/dose(hs,ci))
            uvs = _filterfit(unk, ffp.references, true)
            krs[:, ci] = map(ref.label for ref in ffp.references) do id
                isnan(uvs, id) || (value(uvs, id) < sigma*σ(uvs, id)) ? zero(eltype(krs)) : convert(eltype(krs), value(uvs, id))
            end
        end
        return ThreadsX.map(filter(ii -> ffp.references[ii].label isa CharXRayLabel, eachindex(ffp.references))) do i
            k, lbl = krs[i], ffp.references[i].label
            rprops = properties(spectrum(lbl))
            KRatios( xrays(lbl), properties(hs), rprops, rprops[:Composition], krs[i, :, :])
        end
    else
        @assert false "The mode argument must be :Fast, :Intermediate or :Full."
    end
end


text:
# functions used for curvilinear calculations

"""
  Calculates the S matrices in the cartesian directions for a curvilinear
  element.

  Inputs:
    sbp: an AbstractOperator
    dxidx: the scaled mapping jacobian, dxidx/|J| at every node, Tdim x Txim
           x numNodesPerElement

  Inputs/Outputs:
    S: array to be populated with S, numNodesPerElement x numNodesPerElement x
       x Tdim

"""
function calcSCurvilinear(sbp::AbstractOperator, dxidx::AbstractArray{Tmsh, 3}, S::AbstractArray{T, 3}) where {T, Tmsh}
# calculate S in x, y, z directions

  Tdim = size(S, 3)
  nnodes = size(S, 1)
  fill!(S, 0.0)

  #TODO: test performance
  #      maybe doing 2 passed over the data in order would be better?
  for d1=1:Tdim  # cartesian direction
    for d2=1:Tdim  # parametric direction
      for i=1:nnodes
        for j=1:(i - 1)
          S[j, i, d1] += 0.5*dxidx[d2, d1, j ]*sbp.Q[j, i, d2] - 0.5*dxidx[d2, d1, i]*sbp.Q[i, j, d2]
          S[i, j, d1] = -S[j, i, d1]

        end
      end
    end
  end

  return nothing
end

"""
  Reverse mode of [`calcSCurvilinar`](@ref)

  **Inputs**

   * sbp
   * dxidx
   * S
   * S_bar: seed values for reverse mode

  **Outputs**

   * dxidx_bar: updated with reverse propigation of S_bar
"""
function calcSCurvilinear_rev(sbp::AbstractOperator, dxidx::AbstractArray{Tmsh, 3},
                               dxidx_bar::AbstractArray{Tmsh, 3},
                               S::AbstractArray{T, 3}, S_bar::AbstractArray{T, 3}) where {T, Tmsh}

  Tdim = size(S, 3)
  nnodes = size(S, 1)
  fill!(S, 0.0)

  for d1=1:Tdim  # cartesian direction
    for d2=1:Tdim  # parametric direction
      for i=1:nnodes
        for j=1:(i - 1)
          S[j, i, d1] += 0.5*dxidx[d2, d1, j ]*sbp.Q[j, i, d2] - 0.5*dxidx[d2, d1, i]*sbp.Q[i, j, d2]
          S[i, j, d1] = -S[j, i, d1]

          dxidx_bar[d2, d1, j] += 0.5*sbp.Q[j, i, d2]*S_bar[j, i, d1]
          dxidx_bar[d2, d1, i] -= 0.5*sbp.Q[i, j, d2]*S_bar[j, i, d1]

          dxidx_bar[d2, d1, j] -= 0.5*sbp.Q[j, i, d2]*S_bar[i, j, d1]
          dxidx_bar[d2, d1, i] += 0.5*sbp.Q[i, j, d2]*S_bar[i, j, d1]

        end
      end
    end
  end

  return nothing
end

"""
  Calculates the E matrices in the cartesian directions using the metric
  terms.  This is *not* the E used in the final curvilinear discretization.
  Note that this E is not symmetric

  Inputs:
    sbp: an AbstractOperator
    dxidx: the scaled mapping jacobian, dxidx/|J| at every node, Tdim x Txim
           x numNodesPerElement

  Inputs/Outputs:
    E: array to be populated with E, numNodesPerElement x numNodesPerElement x
       x Tdim

"""
function calcECurvilinear(sbp::AbstractOperator, dxidx::AbstractArray{Tmsh, 3}, E::AbstractArray{T, 3}) where {T, Tmsh}
# calculate E in x, y, z directions
# this isn't the E actually used by the code, it is Ex =  E_xi*metrics_xi_x +
# E_eta*metrics_eta_x

  Tdim = size(E, 3)
  nnodes = size(E, 1)
  fill!(E, 0.0)
  for d1=1:Tdim  # cartesian direction
    for d2=1:Tdim  # parametric direction
      for i=1:nnodes
        for j=1:nnodes
          E_d2 = sbp.Q[j, i, d2] + sbp.Q[i, j, d2]
          E[j, i, d1] += E_d2*dxidx[d2, d1, i]
#          E[i, j, d1] = E[j, i, d1]
        end
      end
    end
  end

  return nothing
end

"""
  TODO: this function appears to compute D in the xi directions, not the
  cartesian directions.  Fix this

  Calculates the D matrices in the cartesian directions using the metric
  terms.  This is *not* the D used in the final curvilinear discretization.

  Inputs:
    sbp: an AbstractOperator
    dxidx: the scaled mapping jacobian, dxidx/|J| at every node, Tdim x Txim
           x numNodesPerElement

  Inputs/Outputs:
    D: array to be populated with D, numNodesPerElement x numNodesPerElement x
       x Tdim

"""

function calcDCurvilinear(sbp::AbstractOperator, dxidx::AbstractArray{Tmsh, 3}, D::AbstractArray{T, 3}) where {T, Tmsh}
# calculate D in x, y, z directions
# this isn't the D actually used by the code, it is Dx =  D_xi*metrics_xi_x +
# D_eta*metrics_eta_x

  Tdim = size(D, 3)
  nnodes = size(D, 1)
  fill!(D, 0.0)
  for d1=1:Tdim  # cartesian direction
    for d2=1:Tdim  # parametric direction
      for i=1:nnodes
        for j=1:nnodes
          h_j = 1./sbp.w[j]
          D[j, i, d1] += h_j*sbp.Q[j, i, d2]
        end
      end
    end
  end

  return nothing
end