text:
type RelativePriceError <: CalibrationErrorType end
type PriceError <: CalibrationErrorType end
type ImpliedVolError <: CalibrationErrorType end

type NaiveBasketType <: CalibrationBasketType end
type MaturityStrikeByDeltaGammaBasketType <: CalibrationBasketType end

type CalibrationHelperCommon
  marketValue::Float64
  calibrationErrorType::CalibrationErrorType

  CalibrationHelperCommon() = new(0.0, RelativePriceError())
end

type ImpliedVolatilityHelper{C <: CalibrationHelper}
  helper::C
  value::Float64
end

function operator(iv::ImpliedVolatilityHelper)
  function _inner(x::Float64)
    return iv.value - black_price!(iv.helper, x)
  end

  return _inner
end

calibration_error{C <: CalibrationHelper}(::RelativePriceError, helper::C) =
  abs(helper.calibCommon.marketValue - model_value!(helper)) / helper.calibCommon.marketValue

function implied_volatility!(ch::CalibrationHelper, targetValue::Float64, accuracy::Float64, maxEvals::Int, minVol::Float64, maxVol::Float64)
  ivh = ImpliedVolatilityHelper(ch, targetValue)

  solv = BrentSolver(maxEvals)

  return solve(solv, operator(ivh), accuracy, ch.volatility.value, minVol, maxVol)
end

function update_pricing_engine{C <: CalibrationHelper, P <: PricingEngine}(ch::C, pe::P)
  newCh = clone(ch, pe)

  return newCh
end

type SwaptionHelper{Dm <: Dates.Period, Dl <: Dates.Period, TP <: TenorPeriod, DC_fix <: DayCount, DC_float <: DayCount, T <: YieldTermStructure, P <: PricingEngine} <: CalibrationHelper
  lazyMixin::LazyMixin
  exerciseDate::Date
  endDate::Date
  maturity::Dm
  swapLength::Dl
  volatility::Quote
  iborIndex::IborIndex
  fixedLegTenor::TP
  fixedLegDayCount::DC_fix
  floatingLegDayCount::DC_float
  strike::Float64
  nominal::Float64
  shift::Float64
  exerciseRate::Float64
  calibCommon::CalibrationHelperCommon
  yts::T
  pricingEngine::P
  swaption::Swaption

  SwaptionHelper(lazyMixin::LazyMixin, exerciseDate::Date, endDate::Date, maturity::Dm, swapLength::Dl, volatility::Quote, iborIndex::IborIndex, fixedLegTenor::TP,
                fixedLegDayCount::DC_fix, floatingLegDayCount::DC_float, strike::Float64, nominal::Float64, shift::Float64, exerciseRate::Float64,
                calibCommon::CalibrationHelperCommon, yts::T, pe::P) =
                new{Dm, Dl, TP, DC_fix, DC_float, T, P}(lazyMixin, exerciseDate, endDate, maturity, swapLength, volatility, iborIndex, fixedLegTenor,
                fixedLegDayCount, floatingLegDayCount, strike, nominal, shift, exerciseRate, calibCommon, yts, pe)

  SwaptionHelper(lazyMixin::LazyMixin, exerciseDate::Date, endDate::Date, maturity::Dm, swapLength::Dl, volatility::Quote, iborIndex::IborIndex, fixedLegTenor::TP,
                fixedLegDayCount::DC_fix, floatingLegDayCount::DC_float, strike::Float64, nominal::Float64, shift::Float64, exerciseRate::Float64,
                calibCommon::CalibrationHelperCommon, yts::T, pe::P, swaption::Swaption) =
                new{Dm, Dl, TP, DC_fix, DC_float, T, P}(lazyMixin, exerciseDate, endDate, maturity, swapLength, volatility, iborIndex, fixedLegTenor,
                fixedLegDayCount, floatingLegDayCount, strike, nominal, shift, exerciseRate, calibCommon, yts, pe, swaption)
end

SwaptionHelper{Dm <: Dates.Period, Dl <: Dates.Period, TP <: TenorPeriod, DC_fix <: DayCount, DC_float <: DayCount, T <: YieldTermStructure, P <: PricingEngine}(maturity::Dm,
              swapLength::Dl, volatility::Quote, iborIndex::IborIndex, fixedLegTenor::TP, fixedLegDayCount::DC_fix, floatingLegDayCount::DC_float,
              yts::T, pe::P, strike::Float64 = -1.0, nominal::Float64 = 1.0, shift::Float64 = 0.0, exerciseRate::Float64 = 0.0) =
              SwaptionHelper{Dm, Dl, TP, DC_fix, DC_float, T, P}(LazyMixin(), Date(), Date(), maturity, swapLength, volatility, iborIndex, fixedLegTenor, fixedLegDayCount, floatingLegDayCount, strike,
              nominal, shift, exerciseRate, CalibrationHelperCommon(), yts, pe)

SwaptionHelper{TP <: TenorPeriod, DC_fix <: DayCount, DC_float <: DayCount, T <: YieldTermStructure, P <: PricingEngine}(expiryDate::Date,
              endDate::Date, volatility::Quote, iborIndex::IborIndex, fixedLegTenor::TP, fixedLegDayCount::DC_fix, floatingLegDayCount::DC_float,
              yts::T, pe::P = NullSwaptionEngine(), strike::Float64 = -1.0, nominal::Float64 = 1.0, shift::Float64 = 0.0, exerciseRate::Float64 = 0.0) =
              SwaptionHelper{Dates.Day, Dates.Day, TP, DC_fix, DC_float, T, P}(LazyMixin(), expiryDate, endDate, Dates.Day(0), Dates.Day(0), volatility, iborIndex, fixedLegTenor,
              fixedLegDayCount, floatingLegDayCount, strike, nominal, shift, exerciseRate, CalibrationHelperCommon(), yts, pe)

SwaptionHelper{Dm <: Dates.Period, Dl <: Dates.Period, TP <: TenorPeriod, DC_fix <: DayCount, DC_float <: DayCount, T <: YieldTermStructure, P <: PricingEngine}(lazyMixin::LazyMixin,
              exerciseDate::Date, endDate::Date, maturity::Dm, swapLength::Dl, volatility::Quote, iborIndex::IborIndex, fixedLegTenor::TP,
              fixedLegDayCount::DC_fix, floatingLegDayCount::DC_float, strike::Float64, nominal::Float64, shift::Float64, exerciseRate::Float64,
              calibCommon::CalibrationHelperCommon, yts::T, pe::P, swaption::Swaption) =
              SwaptionHelper{Dm, Dl, TP, DC_fix, DC_float, T, P}(lazyMixin, exerciseDate, endDate, maturity, swapLength, volatility, iborIndex, fixedLegTenor,
              fixedLegDayCount, floatingLegDayCount, strike, nominal, shift, exerciseRate, calibCommon, yts, pe, swaption)

function perform_calculations!(swaptionHelper::SwaptionHelper)
  calendar = swaptionHelper.iborIndex.fixingCalendar
  fixingDays = swaptionHelper.iborIndex.fixingDays
  convention = swaptionHelper.iborIndex.convention

  exerciseDate = swaptionHelper.exerciseDate == Date() ? advance(swaptionHelper.maturity, calendar, swaptionHelper.yts.referenceDate, convention) : swaptionHelper.exerciseDate

  startDate = advance(Dates.Day(fixingDays), calendar, exerciseDate, convention)

  endDate = swaptionHelper.endDate == Date() ? advance(swaptionHelper.swapLength, calendar, startDate, convention) : swaptionHelper.endDate

  fixedSchedule = QuantLib.Time.Schedule(startDate, endDate, swaptionHelper.fixedLegTenor, convention, convention, QuantLib.Time.DateGenerationForwards(), false, calendar)
  floatSchedule = QuantLib.Time.Schedule(startDate, endDate, swaptionHelper.iborIndex.tenor, convention, convention, QuantLib.Time.DateGenerationForwards(), false, calendar)

  swapEngine = DiscountingSwapEngine(swaptionHelper.yts, false)

  swapT = Receiver()

  tempSwap = VanillaSwap(swapT, swaptionHelper.nominal, fixedSchedule, 0.0, swaptionHelper.fixedLegDayCount, swaptionHelper.iborIndex, 0.0, floatSchedule, swaptionHelper.floatingLegDayCount, swapEngine)

  forward = fair_rate(tempSwap)

  if swaptionHelper.strike == -1.0
    swaptionHelper.exerciseRate = forward
  else
    swaptionHelper.exerciseRate = strike
    swapT = strike <= forward ? Receiver() : Payer()
  end

  swap = VanillaSwap(swapT, swaptionHelper.nominal, fixedSchedule, swaptionHelper.exerciseRate, swaptionHelper.fixedLegDayCount, swaptionHelper.iborIndex, 0.0, floatSchedule, swaptionHelper.floatingLegDayCount, swapEngine)
  exercise = EuropeanExercise(exerciseDate)

  swaptionHelper.swaption = Swaption(swap, exercise)

  # calibration calc
  _calibration_calculate!(swaptionHelper)

  return swaptionHelper
end

function _calibration_calculate!(swaptionHelper::SwaptionHelper)
  swaptionHelper.calibCommon.marketValue = black_price!(swaptionHelper, swaptionHelper.volatility.value)

  return swaptionHelper
end

function add_times_to!(swaptionHelper::SwaptionHelper, times::Vector{Float64})
  calculate!(swaptionHelper)
  discretizedSwap = DiscretizedSwaption(swaptionHelper.swaption, reference_date(swaptionHelper.yts), swaptionHelper.yts.dc)
  times = vcat(times, mandatory_times(discretizedSwap))
  return times
end

function model_value!(sh::SwaptionHelper)
  calculate!(sh)
  sh.swaption = update_pricing_engine(sh.swaption, sh.pricingEngine) # this might clone swaption
  return npv(sh.swaption)
end

function black_price!(swaptionHelper::SwaptionHelper, sigma::Float64)
  calculate!(swaptionHelper)
  # stuff
  black = BlackSwaptionEngine(swaptionHelper.yts, Quote(sigma), Actual365(), swaptionHelper.shift)
  swaptionHelper.swaption = update_pricing_engine(swaptionHelper.swaption, black)
  # swaptionHelper.swaption.pricingEngine = black
  value = npv(swaptionHelper.swaption)
  return value
end

function underlying_swap!(swaptionHelper::SwaptionHelper)
  calculate!(swaptionHelper)
  return swaptionHelper.swaption.swap
end

# clone functions #
function clone(swaptionHelper::SwaptionHelper, pe::PricingEngine = swaptionHelper.pricingEngine)
  lazyMixin, calibCommon = pe == swaptionHelper.pricingEngine ? (swaptionHelper.lazyMixin, swaptionHelper.calibCommon) : (LazyMixin(), CalibrationHelperCommon())
  if isdefined(swaptionHelper, :swaption)
    newSwaptionHelper = SwaptionHelper(lazyMixin, swaptionHelper.exerciseDate, swaptionHelper.endDate, swaptionHelper.maturity, swaptionHelper.swapLength,
                        swaptionHelper.volatility, swaptionHelper.iborIndex, swaptionHelper.fixedLegTenor, swaptionHelper.fixedLegDayCount, swaptionHelper.floatingLegDayCount,
                        swaptionHelper.strike, swaptionHelper.nominal, swaptionHelper.shift, swaptionHelper.exerciseRate, calibCommon, swaptionHelper.yts,
                        pe, swaptionHelper.swaption)
  else
    newSwaptionHelper = SwaptionHelper(lazyMixin, swaptionHelper.exerciseDate, swaptionHelper.endDate, swaptionHelper.maturity, swaptionHelper.swapLength,
                        swaptionHelper.volatility, swaptionHelper.iborIndex, swaptionHelper.fixedLegTenor, swaptionHelper.fixedLegDayCount, swaptionHelper.floatingLegDayCount,
                        swaptionHelper.strike, swaptionHelper.nominal, swaptionHelper.shift, swaptionHelper.exerciseRate, calibCommon, swaptionHelper.yts,
                        pe)
  end

  return newSwaptionHelper
end


text:
@testset "kl_divergence" begin
    p = MvNormal(zeros(2), ones(2))
    q = MvNormal(zeros(2), ones(2))
    @test all(kl_divergence(p, q) .≈ 0.0)

    cp = ConditionalMvNormal(identity)

    kld = kl_divergence(cp, q, zeros(2,10))
    @test size(kld) == (1,10)
    @test all(kld .≈ 0.0)
end

ext: jl


text:
relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 3, 4, 2, true, false, false), nothing, relaxed_tols),
    ((2, 3, 4, 2, true, false, false), :SOCExpPSD, relaxed_tols),
    ((2, 3, 4, 2, false, true, false),),
    ((2, 3, 4, 2, false, true, false), :ExpPSD),
    ((2, 3, 4, 2, false, false, true),),
    ((2, 3, 4, 2, false, false, true), :ExpPSD),
    ]
insts["fast"] = [
    ((3, 5, 7, 2, true, false, false),),
    ((3, 5, 7, 2, true, false, false), :SOCExpPSD),
    ((3, 5, 7, 2, false, true, false),),
    ((3, 5, 7, 2, false, true, false), :ExpPSD),
    ((3, 5, 7, 2, false, false, true),),
    ((3, 5, 7, 2, false, false, true), :ExpPSD),
    ((5, 15, 25, 5, true, false, false),),
    ((5, 15, 25, 5, false, true, false),),
    ((5, 15, 25, 5, false, false, true),),
    ((10, 30, 50, 5, true, false, false),),
    ((10, 30, 50, 5, false, true, false),),
    ((10, 30, 50, 5, false, false, true),),
    ((25, 75, 125, 10, true, false, false),),
    ((25, 75, 125, 10, false, true, false),),
    ((25, 75, 125, 10, false, false, true),),
    ]
insts["various"] = [
    ((25, 45, 75, 10, true, false, false),),
    ((25, 45, 75, 10, true, false, false), :SOCExpPSD),
    ((25, 45, 75, 10, false, true, false),),
    ((25, 45, 75, 10, false, true, false), :ExpPSD),
    ((25, 45, 75, 10, false, true, false), :SOCExpPSD),
    ((25, 45, 75, 10, false, false, true),),
    ((25, 45, 75, 10, false, false, true), :ExpPSD),
    ((25, 45, 75, 10, false, false, true), :SOCExpPSD),
    ((50, 90, 150, 10, true, false, false),),
    ((50, 90, 150, 10, true, false, false), :SOCExpPSD),
    ((50, 90, 150, 10, false, true, false),),
    ((50, 90, 150, 10, false, true, false), :ExpPSD),
    ((50, 90, 150, 10, false, true, false), :SOCExpPSD),
    ((50, 90, 150, 10, false, false, true),),
    ((100, 180, 300, 10, true, false, false),),
    ((100, 180, 300, 10, false, true, false),),
    ]
return (DOptimalDesignJuMP, insts)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct FilesCommentsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_filesCommentsDelete(_api::FilesCommentsApi; file=nothing, id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/files.comments.delete", ["slackAuth"])
    Swagger.set_param(_ctx.form, "file", file)  # type String
    Swagger.set_param(_ctx.form, "id", id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Deletes an existing comment on a file.
Param: file::String
Param: id::String
Return: Dict{String, Any}
"""
function filesCommentsDelete(_api::FilesCommentsApi; file=nothing, id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_filesCommentsDelete(_api; file=file, id=id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function filesCommentsDelete(_api::FilesCommentsApi, response_stream::Channel; file=nothing, id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_filesCommentsDelete(_api; file=file, id=id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export filesCommentsDelete


text:
# Import modules to all processes
@everywhere using Pkg
@everywhere Pkg.activate("../..")
@everywhere using Distributions
@everywhere using StatsBase
@everywhere using LinearAlgebra
# Import EKP modules
@everywhere using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
@everywhere using EnsembleKalmanProcesses.Observations
@everywhere using EnsembleKalmanProcesses.ParameterDistributionStorage
@everywhere include(joinpath(@__DIR__, "../../src/helper_funcs.jl"))
using JLD2
using NPZ

#########
#########  Define the parameters and their priors
#########

# Define the parameters that we want to learn
param_names = ["entrainment_factor", "detrainment_factor", "sorting_power",
    "tke_ed_coeff", "tke_diss_coeff", "pressure_normalmode_adv_coeff",
         "pressure_normalmode_buoy_coeff1", "pressure_normalmode_drag_coeff", "static_stab_coeff"]
n_param = length(param_names)

# Prior information: Define transform to unconstrained gaussian space
constraints = [ [bounded(0.01, 0.3)],
                [bounded(0.01, 0.9)],
                [bounded(0.25, 4.0)],
                [bounded(0.01, 0.5)],
                [bounded(0.01, 0.5)],
                [bounded(0.0, 0.5)],
                [bounded(0.0, 0.5)],
                [bounded(5.0, 15.0)],
                [bounded(0.1, 0.8)]]
# All vars are standard Gaussians in unconstrained space
prior_dist = [Parameterized(Normal(0.0, 0.5))
                for x in range(1, n_param, length=n_param) ]
priors = ParameterDistribution(prior_dist, constraints, param_names)

#########
#########  Retrieve true LES samples from PyCLES data and transform
#########

# Define observation window per flow configuration
ti = [7200.0, 25200.0, 7200.0, 14400.0]
tf = [14400.0, 32400.0, 14400.0, 21600.0]

#ti = [ [7200.0, 14400.0], [14400.0, 25200.0], [7200.0, 14400.0], [7200.0, 14400.0] ]
#tf = [ [14400.0, 21600.0], [21600.0, 32400.0], [14400.0, 21600.0], [14400.0, 21600.0] ]

# Define variables per flow configuration
y_names = Array{String, 1}[]
push!(y_names, ["thetal_mean", "ql_mean", "qt_mean", "total_flux_h", "total_flux_qt"]) #DYCOMS_RF01
push!(y_names, ["thetal_mean", "u_mean", "v_mean", "tke_mean"]) #GABLS
push!(y_names, ["thetal_mean", "total_flux_h"]) #Nieuwstadt
push!(y_names, ["thetal_mean", "ql_mean", "qt_mean", "total_flux_h", "total_flux_qt"]) #Bomex

# Define preconditioning and regularization of inverse problem
normalized = true # Variable normalization
perform_PCA = true # PCA on config covariance
cutoff_reg = true # Regularize above PCA cutoff
beta = 10.0 # Regularization hyperparameter
variance_loss = 1.0e-3 # PCA variance loss
noisy_obs = true # Choice of covariance in evaluation of y_{j+1} in EKI. True -> Γy, False -> 0

sim_names = ["DYCOMS_RF01", "GABLS", "Nieuwstadt", "Bomex"]
les_suffixes = ["may20", "iles128wCov", "dry11", "may18"]
les_root = "/groups/esm/ilopezgo"
scm_data_root = pwd()
scm_dir = "/home/ilopezgo/SCAMPy"
save_full_EDMF_data = false  # if true, save each ensemble output file
outdir_root = pwd()

C = sum([length(ti_) for ti_ in ti])
sim_num = length(sim_names)

# Init arrays
yt = zeros(0)
yt_var_list = Array{Float64, 2}[]
yt_big = zeros(0)
yt_var_list_big = Array{Float64, 2}[]
P_pca_list = Array{Float64, 2}[]
pool_var_list = []

@assert (  # Each entry in these lists correspond to one simulation case
     length(sim_names) == length(les_suffixes)
     == length(y_names) == length(ti) == length(tf)
 )

for (i, sim_name) in enumerate(sim_names) # Loop on simulations
    for (j, ti_j) in enumerate(ti[i]) # Loop on time intervals
        tf_j = tf[i][j]
        les_dir = occursin("Nieuwstadt", sim_name) ?
            joinpath(les_root, "Output.Soares.$(les_suffixes[i])") : joinpath(les_root, "Output.$sim_name.$(les_suffixes[i])")
        # Get SCM vertical levels for interpolation
        z_scm = get_profile(joinpath(scm_data_root, "Output.$sim_name.00000"), ["z_half"])
        # Get (interpolated and pool-normalized) observations, get pool variance vector
        yt_, yt_var_, pool_var = obs_LES(y_names[i], les_dir, ti_j, tf_j, z_scm = z_scm, normalize=normalized)
        push!(pool_var_list, pool_var)
        if perform_PCA
            yt_pca, yt_var_pca, P_pca = obs_PCA(yt_, yt_var_, variance_loss)
            append!(yt, yt_pca)
            push!(yt_var_list, yt_var_pca)
            push!(P_pca_list, P_pca)
        else
            append!(yt, yt_)
            push!(yt_var_list, yt_var_)
            global P_pca_list = nothing
        end
        # Save full dimensionality (normalized) output for error computation
        append!(yt_big, yt_)
        push!(yt_var_list_big, yt_var_)
    end
end

d = length(yt) # Length of data array
# Construct global observational covariance matrix for error diagnostic, no TSVD
yt_var_big = cov_from_cov_list(yt_var_list_big)

# Construct global observational covariance matrix for EKP, TSVD
yt_var = zeros(d, d)
vars_num = 1
for (k,config_cov) in enumerate(yt_var_list)
    vars = length(config_cov[1,:])
    yt_var[vars_num:vars_num+vars-1, vars_num:vars_num+vars-1] = cutoff_reg ?
         (config_cov + (beta*minimum(diag(config_cov))) .* Matrix(1.0I, size(config_cov)) ) .* vars : config_cov .* vars
    global vars_num = vars_num+vars
end
println("NUMBER OF OUTPUTS CONSIDERED ", d, " CAPTURING FRACTION OF VARIANCE: ", 1.0-variance_loss)

Γy = yt_var
println("DETERMINANT OF FULL Γy, ", det(Γy))

#########
#########  Calibrate: Ensemble Kalman inversion
#########

algo = Unscented(vcat(get_mean(priors)...), get_cov(priors), 1.0, 0 ) # Sampler(vcat(get_mean(priors)...), get_cov(priors)) # Inversion()
N_ens = typeof(algo) == Unscented{Float64,Int64} ? 2*n_param + 1 : 50 # number of ensemble members
N_iter = 10 # number of EKP iterations.
Δt = 1.0 # follows scaling by batch size

println("NUMBER OF ENSEMBLE MEMBERS: ", N_ens)
println("NUMBER OF ITERATIONS: ", N_iter)
deterministic_forward_map = noisy_obs ? true : false

if typeof(algo) != Unscented{Float64,Int64}
    initial_params = construct_initial_ensemble(priors, N_ens, rng_seed=rand(1:1000))
    # Discard unstable parameter combinations, parallel
    #precondition_ensemble!(initial_params, priors, param_names, y_names, ti, tf=tf)
end

# UKI does not require sampling from the prior
ekobj = typeof(algo) == Unscented{Float64,Int64} ?  EnsembleKalmanProcess(yt, Γy, algo ) : EnsembleKalmanProcess(initial_params, yt, Γy, algo )
@everywhere g_(x::Array{Float64,1}) = run_SCAMPy(
     x, $param_names, $y_names, $scm_dir,
     $scm_data_root, $sim_names, $ti, $tf,
     P_pca_list = $P_pca_list, norm_var_list = $pool_var_list,
 )
# Create output dir
prefix = perform_PCA ? "results_pycles_PCA_" : "results_pycles_"
prefix = cutoff_reg ? string(prefix, "creg", beta, "_") : prefix
prefix = typeof(algo) == Sampler{Float64} ? string(prefix, "eks_") : prefix
prefix = typeof(algo) == Unscented{Float64,Int64} ? string(prefix, "uki_") : prefix
prefix = noisy_obs ? prefix : string(prefix, "nfo_")
prefix = Δt ≈ 1 ? prefix : string(prefix, "dt", Δt, "_")
outdir_path = string(prefix, "p", n_param, "_e", N_ens, "_i", N_iter, "_d", d)
println("Name of outdir path for this EKP, ", outdir_path)
command = `mkdir $outdir_path`
try
    run(command)
catch e
    println("Output directory already exists. Output may be overwritten.")
end

# EKP iterations
g_ens = zeros(N_ens, d)
norm_err_list = []
g_big_list = []
Δt_scaled = Δt / C # Scale artificial timestep by batch size
for i in 1:N_iter
    # Note that the parameters are transformed when used as input to SCAMPy
    params_cons_i = deepcopy(transform_unconstrained_to_constrained(priors,
        get_u_final(ekobj)) )
    params = [row[:] for row in eachrow(params_cons_i')]
    @everywhere params = $params
    array_of_tuples = pmap(g_, params) # Outer dim is params iterator
    (sim_dirs_arr, g_ens_arr, g_ens_arr_pca) = ntuple(l->getindex.(array_of_tuples,l),3) # Outer dim is G̃, G
    println("LENGTH OF G_ENS_ARR", length(g_ens_arr))
    println("LENGTH OF G_ENS_ARR_PCA", length(g_ens_arr_pca))
    println(string("\n\nEKP evaluation $i finished. Updating ensemble ...\n"))
    for j in 1:N_ens
      g_ens[j, :] = g_ens_arr_pca[j]
    end
    # Get normalized error
    push!(norm_err_list, compute_errors(g_ens_arr, yt_big))
    push!(g_big_list, g_ens_arr)
    if typeof(algo) == Inversion
        update_ensemble!(ekobj, Array(g_ens'), Δt_new=Δt_scaled, deterministic_forward_map = deterministic_forward_map)
    elseif typeof(algo) == Unscented{Float64,Int64}
        update_ensemble!(ekobj, Array(g_ens'), Δt_new=Δt_scaled )
    else
        update_ensemble!(ekobj, Array(g_ens') )
    end
    println("\nEnsemble updated. Saving results to file...\n")
    # Save EKP information to file
    save(joinpath(outdir_path, "ekp.jld2"),
        "ekp_u", transform_unconstrained_to_constrained(priors, get_u(ekobj)),
        "ekp_g", get_g(ekobj),
        "truth_mean", ekobj.obs_mean,
        "truth_cov", ekobj.obs_noise_cov,
        "ekp_err", ekobj.err,
        "norm_err", norm_err_list,
        "truth_mean_big", yt_big,
        "g_big", g_big_list,
        "truth_cov_big", yt_var_big,
        "P_pca", P_pca_list,
        "pool_var", pool_var_list,
        )
    # Convert to arrays
    phi_params = Array{Array{Float64,2},1}(transform_unconstrained_to_constrained(priors, get_u(ekobj)))
    phi_params_arr = zeros(i+1, n_param, N_ens)
    g_big_arr = zeros(i, N_ens, length(yt_big))
    for (k,elem) in enumerate(phi_params)
      phi_params_arr[k,:,:] = elem
      if k < i + 1
        g_big_arr[k,:,:] = hcat(g_big_list[k]...)'
      end
    end
    norm_err_arr = hcat(norm_err_list...)' # N_iter, N_ens
    npzwrite(string(outdir_path,"/y_mean.npy"), ekobj.obs_mean)
    npzwrite(string(outdir_path,"/Gamma_y.npy"), ekobj.obs_noise_cov)
    npzwrite(string(outdir_path,"/y_mean_big.npy"), yt_big)
    npzwrite(string(outdir_path,"/Gamma_y_big.npy"), yt_var_big)
    npzwrite(string(outdir_path,"/phi_params.npy"), phi_params_arr)
    npzwrite(string(outdir_path,"/norm_err.npy"), norm_err_arr)
    npzwrite(string(outdir_path,"/g_big.npy"), g_big_arr)
    for (l, P_pca) in enumerate(P_pca_list)
      if C ≈ sim_num
         npzwrite(string(outdir_path,"/P_pca_",sim_names[l],".npy"), P_pca)
         npzwrite(string(outdir_path,"/pool_var_",sim_names[l],".npy"), pool_var_list[l])
      else
         npzwrite(string(outdir_path,"/P_pca_",l,".npy"), P_pca)
         npzwrite(string(outdir_path,"/pool_var_",l,".npy"), pool_var_list[l])
      end
    end

    if save_full_EDMF_data
         # Save full EDMF data from every ensemble
         eki_iter_path = joinpath(outdir_path, "EKI_iter_$i")
         mkpath(eki_iter_path)
         # get a simulation directory `.../Output.SimName.UUID`, and corresponding parameter name
         for (ens_i, sim_dirs) in enumerate(sim_dirs_arr)  # each ensemble returns a list of simulation directories
             ens_i_path = joinpath(eki_iter_path, "ens_$ens_i")
             mkpath(ens_i_path)
             for (sim_name, sim_dir) in zip(sim_names, sim_dirs)
                 # Copy simulation data to output directory
                 dirname = splitpath(sim_dir)[end]
                 @assert dirname[1:7] == "Output."  # sanity check
                 # Stats file
                 tmp_data_path = joinpath(sim_dir, "stats/Stats.$sim_name.nc")
                 save_data_path = joinpath(ens_i_path, "Stats.$sim_name.$ens_i.nc")
                 run(`cp $tmp_data_path $save_data_path`)
                 # namefile and paramfile
                 tmp_namefile_path = joinpath(sim_dir, "$sim_name.in")
                 save_namefile_path = joinpath(ens_i_path, "$sim_name.in")
                 run(`cp $tmp_namefile_path $save_namefile_path`)
                 tmp_paramfile_path = joinpath(sim_dir, "paramlist_$sim_name.in")
                 save_paramfile_path = joinpath(ens_i_path, "paramlist_$sim_name.in")
                 run(`cp $tmp_paramfile_path $save_paramfile_path`)
             end
         end
     end

end

# EKP results: Has the ensemble collapsed toward the truth?
println("\nEKP ensemble mean at last stage (original space):")
println(mean(transform_unconstrained_to_constrained(priors, get_u(ekobj)), dims=1))

println("\nEnsemble covariance det. 1st iteration, transformed space.")
println(det( cov((get_u(ekobj, 1)), dims=1) ))
println("\nEnsemble covariance det. last iteration, transformed space.")
println(det( cov(get_u_final(ekobj), dims=2) ))


text:
module TestCodim2
include("preamble_plots.jl")

using Bifurcations.Codim2: NormalizingAS, BackReferencingAS
using Bifurcations.Examples.Calcium

@testset "smoke Calcium codim-2 ($(nameof(ASType)))" for
        ASType in [NormalizingAS, BackReferencingAS]

    codim1 = init(Calcium.prob)
    solve!(codim1)

    for point in special_intervals(codim1)
        codim2_prob = BifurcationProblem(
            point,
            codim1,
            (@lens _.gca),
            (0.0, 8.0),
            augmented_system = ASType(),
        )
        codim2_solver = init(codim2_prob)
        @test_nothrow solve!(codim2_solver)

        @testset "show" begin
            smoke_test_solver_show(codim2_solver)
        end

        @testset "plot" begin
            smoke_test_solver_plot(codim2_solver)
        end
    end
end

end  # module


text:
module BDPlusExceptions
    function raise_discoverydoc_err()
        printstyled("ERROR: ", color=:red, bold=true)
        printstyled("Unable to discover BigQuery REST API documentation. \n",
        "Please contact the developer reporting this error.", bold=true)
    end
end

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
#### Caching
function CachedOperator{T}(io::InterlaceOperator{T,1};padding::Bool=false)
    ds=domainspace(io)
    rs=rangespace(io)

    ind=find(op->isinf(size(op,1)),io.ops)
    if length(ind) ≠ 1  || !isbanded(io.ops[ind[1]])  # is almost banded
        return CachedOperator(io,Array(eltype(op),0,0))
    end
    i=ind[1]
    bo=io.ops[i]
    lin,uin=bandwidths(bo)

    # calculate number of rows interlaced
    # each each row in the rangespace increases the lower bandwidth by 1
    nds=0
    md=0
    for k=1:length(io.ops)
        if k ≠ i
            d=dimension(rs[k])
            nds+=d
            md=max(md,d)
        end
    end

    isend=true
    for k=i+1:length(io.ops)
        if dimension(rs[k]) == md
            isend=false
        end
    end

    numoprows=isend?md-1:md
    n=nds+numoprows

    (l,u) = (max(lin+nds,n-1),max(0,uin+1-ind[1]))

    # add extra rows for QR
    if padding
        u+=l
    end

    ret=abzeros(T,n,n+u,l,u,nds)

    # populate the finite rows
    jr=1:n+u
    bcrow=1
    oprow=0
    for k=1:n
        K,J=io.rangeinterlacer[k]

        if K ≠ i
            # fill the fill matrix
            ret.fill.V[:,bcrow] = Matrix(view(io.ops[K],J:J,jr))
            ret.fill.U[k,bcrow] = 1
            bcrow += 1
        else
            oprow+=1
        end

        for j=rowrange(ret.bands,k)
            ret[k,j] = io[k,j]
        end
    end

    CachedOperator(io,ret,(n,n+u),ds,rs,(l,∞))
end

function CachedOperator{T}(io::InterlaceOperator{T,2};padding::Bool=false)
    ds=domainspace(io)
    rs=rangespace(io)
    di=io.domaininterlacer
    ri=io.rangeinterlacer
    ddims=di.iterator.dimensions
    rdims=ri.iterator.dimensions

    # we are only almost banded if every operator is either finite
    # range or banded, and if the # of  ∞ spaces is the same
    # for the domain and range
    isab=all(op->isfinite(size(op,1)) || isbanded(op),io.ops) &&
        count(isinf,ddims) == count(isinf,rdims)

    if !isab
        return CachedOperator(io,Array(T,0,0))
    end

    # these are the bandwidths if we only had ∞-dimensional operators
    d∞=find(isinf,[ddims...])
    r∞=find(isinf,[rdims...])
    p=length(d∞)

    l∞,u∞ = 0,0
    for k=1:p,j=1:p
        l∞=min(l∞,p*bandinds(io.ops[r∞[k],d∞[j]],1)+j-k)
    end
    for k=1:p,j=1:p
        u∞=max(u∞,p*bandinds(io.ops[r∞[k],d∞[j]],2)+j-k)
    end

    # now we move everything by the finite rank
    ncols=mapreduce(d->isfinite(d)?d:0,+,ddims)
    nbcs=mapreduce(d->isfinite(d)?d:0,+,rdims)
    shft=ncols-nbcs
    l∞,u∞=l∞+shft,u∞+shft

    # iterate through finite rows to find worst case bandinds
    l,u=l∞,u∞
        for k=1+nbcs+p,j=1:ncols+p
            N,n=ri[k]
            M,m=di[j]
            l=min(l,bandinds(io.ops[N,M],1)+n-m-k+j)
            u=max(u,bandinds(io.ops[N,M],2)+n-m-k+j)
        end

    # add extra rows for QR
    if padding
        u-=l
    end

    n=1+nbcs+p
    ret=abzeros(T,n,n+u,-l,u,nbcs)

    # populate entries and fill functionals
    bcrow=1
    oprow=0
    jr=1:n+u
    for k=1:n
        K,J=io.rangeinterlacer[k]

        if isfinite(rdims[K] )
            # fill the fill matrix
            ret.fill.V[:,bcrow] = Matrix(view(io,k:k,jr))
            ret.fill.U[k,bcrow] = 1
            bcrow += 1
        else
            oprow+=1
        end

        for j=rowrange(ret.bands,k)
            ret[k,j] = io[k,j]
        end
    end

    CachedOperator(io,ret,(n,n+u),ds,rs,(l,∞))
end

function resizedata!{T<:Number,DS,RS,DI,RI,BI}(co::CachedOperator{T,AlmostBandedMatrix{T},
                                                        InterlaceOperator{T,1,DS,RS,DI,RI,BI}},
                                     n::Integer,::Colon)
    if n ≤ co.datasize[1]
        return co
    end

    (l,u)=bandwidths(co.data.bands)
    pad!(co.data,n,n+u)

    r=rank(co.data.fill)
    ind=findfirst(op->isinf(size(op,1)),co.op.ops)

    k=1
    for (K,J) in co.op.rangeinterlacer
        if K ≠ ind
            co.data.fill.V[co.datasize[2]:end,k] = co.op.ops[K][J,co.datasize[2]:n+u]
            k += 1
            if k > r
                break
            end
        end
    end

    kr=co.datasize[1]+1:n
    jr=max(1,kr[1]-l):n+u
    BLAS.axpy!(1.0,view(co.op.ops[ind],kr-r,jr),
                    view(co.data.bands,kr,jr))

    co.datasize=(n,n+u)
    co
end

function resizedata!{T<:Number,DS,RS,DI,RI,BI}(co::CachedOperator{T,AlmostBandedMatrix{T},
                                                        InterlaceOperator{T,2,DS,RS,DI,RI,BI}},
                                     n::Integer,::Colon)
    if n ≤ co.datasize[1]
        return co
    end

    io=co.op
    ds=domainspace(io)
    rs=rangespace(io)
    di=io.domaininterlacer
    ri=io.rangeinterlacer
    ddims=di.iterator.dimensions
    rdims=ri.iterator.dimensions

    d∞=find(isinf,[ddims...])
    r∞=find(isinf,[rdims...])
    p=length(d∞)

    (l,u)=bandwidths(co.data.bands)
    pad!(co.data,n,n+u)
    co.data
    # r is number of extra rows, ncols is number of extra columns
    r=rank(co.data.fill)
    ncols=mapreduce(d->isfinite(d)?d:0,+,ddims)

    # fill rows
    K=k=1
    while k ≤ r
        if isfinite(dimension(rs[ri[K][1]]))
            co.data.fill.V[co.datasize[2]:end,k] = co.op[K,co.datasize[2]:n+u]
            k += 1
        end
        K += 1
    end

    kr=co.datasize[1]+1:n
    jr=max(1,kr[1]-l):n+u
    io∞=InterlaceOperator(io.ops[r∞,d∞])

    BLAS.axpy!(1.0,view(io∞,kr-r,jr-ncols),view(co.data.bands,kr,jr))

    co.datasize=(n,n+u)
    co
end

resizedata!{T<:Number,DS,RS,DI,RI,BI}(co::CachedOperator{T,AlmostBandedMatrix{T},
                                                        InterlaceOperator{T,1,DS,RS,DI,RI,BI}},
                     n::Integer,m::Integer) = resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:)

resizedata!{T<:Number,DS,RS,DI,RI,BI}(co::CachedOperator{T,AlmostBandedMatrix{T},
                                                         InterlaceOperator{T,2,DS,RS,DI,RI,BI}},
                      n::Integer,m::Integer) = resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:)

##
# These give highly optimized routines for delaying with Cached
#

## populate data

function resizedata!{T,MM,DS,RS,BI}(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T},
                                                                 MM,DS,RS,BI}},
                        ::Colon,col)
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R
    W=QR.H

    R=MO.data.bands
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    F=MO.data.fill.U

    for k=QR.ncols+1:col
        W[:,k] = view(R.data,R.u+1:R.u+R.l+1,k) # diagonal and below
        wp=view(W,:,k)
        W[1,k]+= flipsign(norm(wp),W[1,k])
        normalize!(wp)

        # scale banded entries
        for j=k:k+R.u
            dind=R.u+1+k-j
            v=view(R.data,dind:dind+M-1,j)
            dt=dot(wp,v)
            Base.axpy!(-2*dt,wp,v)
        end

        # scale banded/filled entries
        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=view(R.data,1:M-p,j)  # shift down each time
            wp2=view(wp,p+1:M)
            dt=dot(wp2,v)
            for ℓ=k:k+p-1
                @inbounds dt=muladd(conj(W[ℓ-k+1,k]),
                                    unsafe_getindex(MO.data.fill,ℓ,j),dt)
            end
            Base.axpy!(-2*dt,wp2,v)
        end

        # scale filled entries

        for j=1:size(F,2)
            v=view(F,k:k+M-1,j) # the k,jth entry of F
            dt=dot(wp,v)
            Base.axpy!(-2*dt,wp,v)
        end
    end
    QR.ncols=col
    QR
end

function resizedata!{T,MM,DS,RS,BI}(QR::QROperator{CachedOperator{T,BandedMatrix{T},
                                                                 MM,DS,RS,BI}},
                        ::Colon,col)
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R
    W=QR.H

    R=MO.data
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    for k=QR.ncols+1:col
        W[:,k] = view(R.data,R.u+1:R.u+R.l+1,k) # diagonal and below
        wp=view(W,:,k)
        W[1,k]+= flipsign(norm(wp),W[1,k])
        normalize!(wp)

        # scale banded entries
        for j=k:k+R.u
            dind=R.u+1+k-j
            v=view(R.data,dind:dind+M-1,j)
            dt=dot(wp,v)
            Base.axpy!(-2*dt,wp,v)
        end

        # scale banded/filled entries
        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=view(R.data,1:M-p,j)  # shift down each time
            wp2=view(wp,p+1:M)
            dt=dot(wp2,v)
            Base.axpy!(-2*dt,wp2,v)
        end
    end
    QR.ncols=col
    QR
end

function resizedata!{T,MM,DS,RS,BI}(QR::QROperator{CachedOperator{T,RaggedMatrix{T},
                                                                 MM,DS,RS,BI}},
                        ::Colon,col)
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R
    W=QR.H

    if col ≥ MO.datasize[2]
        m = MO.datasize[2]
        resizedata!(MO,:,col+100)  # double the last rows

        # apply previous Householders to new columns of R
        for J=1:size(W,2)
            wp=view(W,:,J)
            for j=m+1:MO.datasize[2]
                kr=j:j+length(wp)-1
                v=view(MO,kr,j)
                dt=dot(wp,v)
                Base.axpy!(-2*dt,wp,v)
            end
        end
    end

    if col > size(W,2)
        m=size(W,2)
        resize!(W.cols,2col+1)

        for j=m+1:2col
            cs=colstop(MO,j)
            W.cols[j+1]=W.cols[j] + cs-j+1
            W.m=max(W.m,cs-j+1)
        end

        resize!(W.data,W.cols[end]-1)
    end

    for k=QR.ncols+1:col
        cs = colstop(QR.R,k)
        W[1:cs-k+1,k] = view(MO.data,k:cs,k) # diagonal and below
        wp=view(W,1:cs-k+1,k)
        W[1,k]+= flipsign(norm(wp),W[1,k])
        normalize!(wp)

        # scale rows entries
        kr=k:k+length(wp)-1
        for j=k:MO.datasize[2]
            v=view(MO.data,kr,j)
            dt=dot(wp,v)
            Base.axpy!(-2*dt,wp,v)
        end
    end
    QR.ncols=col
    QR
end

# BLAS versions, requires BlasFloat

function resizedata!{T<:BlasFloat,MM,DS,RS,BI}(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T},
                                                                 MM,DS,RS,BI}},
                        ::Colon,col)
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R
    W=QR.H

    R=MO.data.bands
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    F=MO.data.fill.U

    f=pointer(F)
    m,n=size(R)
    w=pointer(W)
    r=pointer(R.data)
    sz=sizeof(T)
    st=stride(R.data,2)
    stw=stride(W,2)

    for k=QR.ncols+1:col
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M,wp,1),W[1,k])
        normalize!(M,wp)

        for j=k:k+R.u
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt=dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end

        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=r+sz*((j-1)*st)  # shift down each time
            dt=dot(M-p,wp+p*sz,1,v,1)
            for ℓ=k:k+p-1
                @inbounds dt=muladd(conj(W[ℓ-k+1,k]),
                                    unsafe_getindex(MO.data.fill,ℓ,j),dt)
            end
            BLAS.axpy!(M-p,-2*dt,wp+p*sz,1,v,1)
        end

        fp=f+(k-1)*sz
        fst=stride(F,2)
        for j=1:size(F,2)
            v=fp+fst*(j-1)*sz   # the k,jth entry of F
            dt=dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end
    end
    QR.ncols=col
    QR
end

function resizedata!{T<:BlasFloat,MM,DS,RS,BI}(QR::QROperator{CachedOperator{T,BandedMatrix{T},
                                                                 MM,DS,RS,BI}},
                        ::Colon,col)
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R
    W=QR.H

    R=MO.data
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    m,n=size(R)
    w=pointer(W)
    r=pointer(R.data)
    sz=sizeof(T)
    st=stride(R.data,2)
    stw=stride(W,2)

    for k=QR.ncols+1:col
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M,wp,1),W[1,k])
        normalize!(M,wp)

        for j=k:k+R.u
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt=dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end

        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=r+sz*((j-1)*st)  # shift down each time
            dt=dot(M-p,wp+p*sz,1,v,1)
            BLAS.axpy!(M-p,-2*dt,wp+p*sz,1,v,1)
        end
    end
    QR.ncols=col
    QR
end

# back substitution
trtrs!(::Type{Val{'U'}},co::CachedOperator,u::Array) =
                trtrs!(Val{'U'},resizedata!(co,size(u,1),size(u,1)).data,u)

function trtrs!(::Type{Val{'U'}},B::AlmostBandedMatrix,u::Array)
    n=size(u,1)
    A=B.bands
    F=B.fill
    b=bandwidth(A,2)
    nbc = rank(B.fill)
    T=eltype(u)
    pk = zeros(T,nbc)

    for c=1:size(u,2)
        fill!(pk,zero(T))

        # before we get to filled rows
        for k=n:-1:max(1,n-b)
            @simd for j=k+1:n
                @inbounds u[k,c]=muladd(-A.data[k-j+A.u+1,j],u[j,c],u[k,c])
            end

            @inbounds u[k,c] /= A.data[A.u+1,k]
        end

       #filled rows
        for k=n-b-1:-1:1
            @simd for j=1:nbc
                @inbounds pk[j] = muladd(u[k+b+1,c],F.V[k+b+1,j],pk[j])
            end

            @simd for j=k+1:k+b
                @inbounds u[k,c]=muladd(-A.data[k-j+A.u+1,j],u[j,c],u[k,c])
            end

            @simd for j=1:nbc
                @inbounds u[k,c] = muladd(-F.U[k,j],pk[j],u[k,c])
            end

            @inbounds u[k,c] /= A.data[A.u+1,k]
        end
    end
    u
end

function trtrs!(::Type{Val{'U'}},A::BandedMatrix,u::Array)
    n=size(u,1)
    b=bandwidth(A,2)
    T=eltype(u)

    for c=1:size(u,2)
        for k=n:-1:1
            @simd for j=k+1:min(n,k+b)
                @inbounds u[k,c]=muladd(-A.data[k-j+A.u+1,j],u[j,c],u[k,c])
            end

            @inbounds u[k,c] /= A.data[A.u+1,k]
        end
    end
    u
end

function trtrs!(::Type{Val{'U'}},A::RaggedMatrix,u::Array)
    if size(A,1) < size(u,1)
        throw(BoundsError())
    end

    n=size(u,1)
    b=bandwidth(A,2)
    T=eltype(u)

    for c=1:size(u,2)
        for k=n:-1:1
            @inbounds ck = A.cols[k]
            @inbounds u[k,c] /= A.data[ck+k-1]
            BLAS.axpy!(-u[k,c],view(A.data,ck:ck+k-2),view(u,1:k-1,c))
        end
    end
    u
end

## Ac_mul_B! for QROperatorQ

function Ac_mul_Bpars(A::QROperatorQ,B::Vector,tolerance,maxlength)
    T = promote_type(eltype(A),eltype(B))
    Ac_mul_Bpars(Operator{T}(A),Vector{T}(B),tolerance,maxlength)
end

function Ac_mul_Bpars{QR,T}(A::QROperatorQ{QR,T},B::Vector{T},tolerance,maxlength)
    if length(B) > A.QR.ncols
        # upper triangularize extra columns to prepare for \
        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
    end

    H=A.QR.H
    M=size(H,1)
    m=length(B)
    Y=pad(B,m+M+10)

    k=1
    yp=view(Y,1:M)
    while (k ≤ m+M || norm(yp) > tolerance )
        if k > maxlength
            warn("Maximum length $maxlength reached.")
            break
        end

        if k+M-1>length(Y)
            pad!(Y,2*(k+M))
        end
        if k > A.QR.ncols
            # upper triangularize extra columns to prepare for \
            resizedata!(A.QR,:,2*(k+M))
            H=A.QR.H
        end

        wp=view(H,:,k)
        yp=view(Y,k:k+M-1)

        dt=dot(wp,yp)
        Base.axpy!(-2*dt,wp,yp)
        k+=1
    end
    Fun(resize!(Y,k),domainspace(A))  # chop off zeros
end

# BLAS

function Ac_mul_Bpars{RR,T<:BlasFloat}(A::QROperatorQ{QROperator{RR,Matrix{T},T},T},
                                        B::Vector{T},
                                        tolerance,maxlength)
    if length(B) > A.QR.ncols
        # upper triangularize extra columns to prepare for \
        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
    end

    H=A.QR.H
    h=pointer(H)

    M=size(H,1)

    b=pointer(B)
    st=stride(H,2)

    sz=sizeof(T)

    m=length(B)
    Y=pad(B,m+M+10)
    y=pointer(Y)

    k=1
    yp=y
    while (k ≤ m+M || BLAS.nrm2(M,yp,1) > tolerance )
        if k > maxlength
            warn("Maximum length $maxlength reached.")
            break
        end

        if k+M-1>length(Y)
            pad!(Y,2*(k+M))
            y=pointer(Y)
        end
        if k > A.QR.ncols
            # upper triangularize extra columns to prepare for \
            resizedata!(A.QR,:,2*(k+M))
            H=A.QR.H
            h=pointer(H)
        end

        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M,wp,1,yp,1)
        BLAS.axpy!(M,-2*dt,wp,1,yp,1)
        k+=1
    end
    Fun(resize!(Y,k),domainspace(A))  # chop off zeros
end


text:
using Jive
runtests(@__DIR__, skip=["loltools/helpers.jl"])


text:
function template(t::String)
    """
    <!DOCTYPE html>

            <title>BAWJ</title>
            <meta charset="utf-8">
            <meta name="viewport" content="width=device-width, initial-scale=1">
            <link rel="shortcut icon" href="/favicon.ico" type="image/x-icon" />
            <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/4.5.2/css/bootstrap.min.css">

            AppliGate

            <button class="navbar-toggler" type="button" data-toggle="collapse" data-target="#collapsibleNavbar">

            </button>

                            Home

                            Aging Report

            $(t)

    """
end

function index(c::WebController)
    render(HTML, template("
Hello World!
"))
end

function aging_report(c::WebController)
    r = @fetchfrom ar_pid report()
    result =
  """

Aging Report

    Invoice	Customer	Date	Amount	Age
  """
    for n = 1:length(r)
      result = result * """

          $(r[n].id_inv)
          $(r[n].csm)	$(r[n].inv_date)
          $(r[n].amount)
          $(r[n].days)

"""
    end
    result * "
"
    render(HTML, template("$(result)"))
end


text:
using JLD2, FileIO
using ClimateMachine
using ClimateMachine.SingleStackUtils
using ClimateMachine.Checkpoint
using ClimateMachine.BalanceLaws: vars_state
import ClimateMachine.BalanceLaws: projection
import ClimateMachine.DGMethods
using ClimateMachine.Atmos
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
import CLIMAParameters

include(joinpath(clima_dir, "experiments", "AtmosLES", "stable_bl_model.jl"))
include("edmf_model.jl")
include("edmf_kernels.jl")

# CLIMAParameters.Planet.T_surf_ref(::EarthParameterSet) = 290.0 # default
CLIMAParameters.Planet.T_surf_ref(::EarthParameterSet) = 265
CLIMAParameters.Atmos.EDMF.a_surf(::EarthParameterSet) = 0.0
function set_clima_parameters(filename)
    eval(:(include($filename)))
end

"""
    init_state_prognostic!(
            turbconv::EDMF{FT},
            m::AtmosModel{FT},
            state::Vars,
            aux::Vars,
            localgeo,
            t::Real,
        ) where {FT}

Initialize EDMF state variables.
This method is only called at `t=0`.
"""
function init_state_prognostic!(
    turbconv::EDMF{FT},
    m::AtmosModel{FT},
    state::Vars,
    aux::Vars,
    localgeo,
    t::Real,
) where {FT}
    # Aliases:
    gm = state
    en = state.turbconv.environment
    up = state.turbconv.updraft
    N_up = n_updrafts(turbconv)
    # GCM setting - Initialize the grid mean profiles of prognostic variables (ρ,e_int,q_tot,u,v,w)
    z = altitude(m, aux)

    # SCM setting - need to have separate cases coded and called from a folder - see what LES does
    # a thermo state is used here to convert the input θ to e_int profile
    e_int = internal_energy(m, state, aux)
    param_set = parameter_set(m)
    ts = PhaseDry(param_set, e_int, state.ρ)
    T = air_temperature(ts)
    p = air_pressure(ts)
    q = PhasePartition(ts)
    θ_liq = liquid_ice_pottemp(ts)

    a_min = turbconv.subdomains.a_min
    @unroll_map(N_up) do i
        up[i].ρa = gm.ρ * a_min
        up[i].ρaw = gm.ρu[3] * a_min
        up[i].ρaθ_liq = gm.ρ * a_min * θ_liq
        up[i].ρaq_tot = FT(0)
    end

    # initialize environment covariance with zero for now
    if z <= FT(250)
        en.ρatke =
            gm.ρ *
            FT(0.4) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0)
        en.ρaθ_liq_cv =
            gm.ρ *
            FT(0.4) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0)
    else
        en.ρatke = FT(0)
        en.ρaθ_liq_cv = FT(0)
    end
    en.ρaq_tot_cv = FT(0)
    en.ρaθ_liq_q_tot_cv = FT(0)
    return nothing
end;

function main(::Type{FT}, cl_args) where {FT}

    surface_flux = cl_args["surface_flux"]

    # Choice of compressibility and CFL
    # compressibility = Compressible()
    compressibility = Anelastic1D()
    str_comp = compressibility == Compressible() ? "COMPRESS" : "ANELASTIC"

    # DG polynomial order
    N = 4
    nelem_vert = 20

    # Prescribe domain parameters
    zmax = FT(400)
    # Simulation time
    t0 = FT(0)
    timeend = FT(1800 * 1)
    CFLmax = compressibility == Compressible() ? FT(1) : FT(100)

    config_type = SingleStackConfigType
    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    # Choice of SGS model
    N_updrafts = 1
    N_quad = 3
    turbconv = NoTurbConv()
    # turbconv = EDMF(
    #     FT,
    #     N_updrafts,
    #     N_quad,
    #     param_set,
    #     surface = NeutralDrySurfaceModel{FT}(param_set),
    # )

    C_smag_ = C_smag(param_set)
    # turbulence = ConstantKinematicViscosity(FT(0.1))
    turbulence = SmagorinskyLilly{FT}(C_smag_)
    model = stable_bl_model(
        FT,
        config_type,
        zmax,
        surface_flux;
        turbulence = turbulence,
        turbconv = turbconv,
        compressibility = compressibility,
    )

    # Assemble configuration
    driver_config = ClimateMachine.SingleStackConfiguration(
        string("SBL_", str_comp, "_1D"),
        N,
        nelem_vert,
        zmax,
        param_set,
        model;
        hmax = FT(40),
        solver_type = ode_solver_type,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFLmax,
    )

    # --- Zero-out horizontal variations:
    vsp = vars_state(model, Prognostic(), FT)
    horizontally_average!(
        driver_config.grid,
        solver_config.Q,
        varsindex(vsp, :turbconv),
    )
    horizontally_average!(
        driver_config.grid,
        solver_config.Q,
        varsindex(vsp, :energy, :ρe),
    )
    vsa = vars_state(model, Auxiliary(), FT)
    horizontally_average!(
        driver_config.grid,
        solver_config.dg.state_auxiliary,
        varsindex(vsa, :turbconv),
    )
    # ---

    dgn_config = config_diagnostics(driver_config)

    # boyd vandeven filter
    cb_boyd = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            ("energy.ρe",),
            solver_config.dg.grid,
            BoydVandevenFilter(
                solver_config.dg.grid,
                1, #default=0
                4, #default=32
            ),
        )
        nothing
    end

    diag_arr = [single_stack_diagnostics(solver_config)]
    time_data = FT[0]

    # Define the number of outputs from `t0` to `timeend`
    n_outputs = 5
    # This equates to exports every ceil(Int, timeend/n_outputs) time-step:
    every_x_simulation_time = ceil(Int, timeend / n_outputs)

    cb_data_vs_time =
        GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
            diag_vs_z = single_stack_diagnostics(solver_config)

            nstep = getsteps(solver_config.solver)

            push!(diag_arr, diag_vs_z)
            push!(time_data, gettime(solver_config.solver))
            nothing
        end

    # Mass tendencies = 0 for Anelastic1D model,
    # so mass should be completely conserved:
    check_cons =
        (ClimateMachine.ConservationCheck("ρ", "3000steps", FT(0.00000001)),)

    cb_print_step = GenericCallbacks.EveryXSimulationSteps(100) do
        @show getsteps(solver_config.solver)
        nothing
    end

    if !isnothing(cl_args["cparam_file"])
        ClimateMachine.Settings.output_dir = cl_args["cparam_file"] * ".output"
    end
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        check_cons = check_cons,
        user_callbacks = (cb_boyd, cb_data_vs_time, cb_print_step),
        check_euclidean_distance = true,
    )

    diag_vs_z = single_stack_diagnostics(solver_config)
    push!(diag_arr, diag_vs_z)
    push!(time_data, gettime(solver_config.solver))

    return solver_config, diag_arr, time_data
end

# ArgParse in global scope to modify Clima Parameters
sbl_args = ArgParseSettings(autofix_names = true)
add_arg_group!(sbl_args, "StableBoundaryLayer")
@add_arg_table! sbl_args begin
    "--cparam-file"
    help = "specify CLIMAParameters file"
    arg_type = Union{String, Nothing}
    default = nothing

    "--surface-flux"
    help = "specify surface flux for energy and moisture"
    metavar = "prescribed|bulk|custom_sbl"
    arg_type = String
    default = "custom_sbl"
end

cl_args = ClimateMachine.init(parse_clargs = true, custom_clargs = sbl_args)
if !isnothing(cl_args["cparam_file"])
    filename = cl_args["cparam_file"]
    set_clima_parameters(filename)
end

solver_config, diag_arr, time_data = main(Float64, cl_args)

include(joinpath(@__DIR__, "report_mse_sbl_anelastic.jl"))

nothing


text:
using LinearAlgebra

function time_slice_spectrum()
    H_matrix = kinetic_matrix_imgtime * auxiliary_field_matrix() * kinetic_matrix_imgtime
end

function sub_partition_weight(L)
    matrix_product = Matrix{Complex{Float64}}(I, num_sites, num_sites)
    for i = 1 : L
        matrix_product *= time_slice_spectrum()
    end
    #return matrix_product
    #if flag_background_subtraction == 1
        #matrix_product *= exp(sum(β * (U/2) * background_dist.^2))
        #return matrix_product
    #end
    effective_spectrum, u = eigen(matrix_product) # spectrum element: e^(-βϵᵢ)
    u_inv = inv(u)
    U_array = complex(zeros(num_sites, num_sites, num_sites))
    for i = 1 : num_sites
        for j = 1 : num_sites
            for k = 1 : num_sites
                U_array[i, j, k] = u_inv[k, i] * u[j, k]
            end
        end
    end
    return effective_spectrum, U_array
end

function recursion_partition_function(effective_spectrum)
    Z = complex(zeros(num_bosons))
    zₖ = complex(zeros(num_bosons))
    for k = 1 : num_bosons
        zₖ[k] = sum(effective_spectrum .^ k)
    end
    for N = 1 : num_bosons
        for M = 1 : N
            if M == N
                Z[N] += 1 * zₖ[M]
            else
                Z[N] += Z[N-M] * zₖ[M]
            end
        end
        Z[N] = Z[N] / N
    end
    return Z
end

function recursion_OccNum(Z, effective_spectrum) # Calculating <nᵢ>
    n_i = complex(zeros(num_bosons, num_sites))
    n_i[1, :] = effective_spectrum
    for N = 2 : num_bosons
        n_i[N, :] = effective_spectrum .* (Z[N-1] .+ n_i[N-1, :])
    end
    return n_i
end

function recursion_double_OccNum_eq(Z, effective_spectrum, n_i) # Calculating <nᵢ^2>
    n²_i = complex(zeros(num_bosons, num_sites)) # <nᵢ^2> array
    n²_i[1, :] = effective_spectrum
    for N = 2 : num_bosons
        n²_i[N, :] = effective_spectrum .* (Z[N-1] .+ 2 * n_i[N-1, :] .+ n²_i[N-1, :])
    end
    return n²_i[num_bosons, :]
end

function recursion_double_OccNum_neq(Z, effective_spectrum, n_i) # Calculating <nᵢnⱼ>
    n²_i = complex(zeros(num_bosons, num_sites, num_sites)) # <nᵢnⱼ> array
    for i = 1 : num_sites
        for j = i + 1 : num_sites
            for N = 2 : num_bosons
                n²_i[N, i, j] = effective_spectrum[j] / 2 * (n_i[N-1, i] + n²_i[N-1, i, j]) +
                            + effective_spectrum[i] / 2 * (n_i[N-1, j] + n²_i[N-1, i, j])
            end
            n²_i[num_bosons, j, i] = n²_i[num_bosons, i, j]
        end
    end
    return n²_i[num_bosons, :, :]
end

#function measure_condensate_fraciton()
#end

function measure_TotEnergy(num_iterations)
    partition_func_array = Vector{Complex{Float64}}()
    Ek_array = Vector{Complex{Float64}}()
    Ep_array = Vector{Complex{Float64}}()
    E_tot_array = Vector{Complex{Float64}}()
    for N = 1 : num_iterations
        effective_spectrum, U_array = sub_partition_weight(L)
        Z = recursion_partition_function(effective_spectrum)
        n_i = recursion_OccNum(Z, effective_spectrum)
        n²_i_diag = recursion_double_OccNum_eq(Z, effective_spectrum, n_i)
        n²_i_offdiag = recursion_double_OccNum_neq(Z, effective_spectrum, n_i)
        E_onebody, E_twobody = 0, 0
        # one-body energy calculation
        for i = 1 : num_sites
            for j = 1 : num_sites
                if abs(i - j) == 1 || abs(i - j) == num_sites - 1
                    D_ij = sum(U_array[i, j, :] .* n_i[num_bosons, :])
                    E_onebody += -t * D_ij
                end
            end
        end
        # two-body energy calculation
        for i = 1 : num_sites
            U_matrix = complex(zeros(num_sites, num_sites))
            D_iiii = 0
            for p = 1 :num_sites
                for r = 1 : num_sites
                    if p != r
                        U_matrix[p, r] = U_array[i, i, p] * U_array[i, i, r]
                        D_iiii += U_matrix[p, r] * n_i[num_bosons, p]
                    end
                end
                D_iiii += U_array[i, i, p]^2 * n²_i_diag[p]
            end
            D_iiii += sum( 2 * U_matrix .* n²_i_offdiag)
            E_twobody += U/2 * D_iiii
        end
        push!(partition_func_array, Z[num_bosons])
        push!(Ek_array, E_onebody)
        push!(Ep_array, E_twobody)
        push!(E_tot_array, E_onebody + E_twobody)
    end
    Z_mean = mean(partition_func_array)
    Ek_mean = mean(Ek_array) / Z_mean
    Ek_error = std(Ek_array) / (sqrt(num_iterations) * real(Z_mean))
    Ep_mean = mean(Ep_array) / Z_mean
    Ep_error = std(Ep_array) / (sqrt(num_iterations) * real(Z_mean))
    E_tot_mean = mean(E_tot_array) / Z_mean
    E_tot_error_bar = std(E_tot_array) / (sqrt(num_iterations) * real(Z_mean))
    return Ek_mean, Ek_error, Ep_mean, Ep_error, E_tot_mean, E_tot_error_bar
end

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("Unit tests confirm correct algorithm implementation") do

    context("pending facts are pending") do
        @pending foo(12) --> "foo"
    end
end


text:
# Writer the byte b in w.
function write(w::WritableFile, b::UInt8)
    write(w, Ref(b))
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: pricing
using AWS.Compat
using AWS.UUIDs

"""
    describe_services()
    describe_services(params::Dict{String,<:Any})

Returns the metadata for one service or a list of the metadata for all services. Use this
without a service code to get the service codes for all services. Use it with a service
code, such as AmazonEC2, to get information specific to that service, such as the attribute
names available for that service. For example, some of the attribute names available for
EC2 are volumeType, maxIopsVolume, operation, locationType, and instanceCapacity10xlarge.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"FormatVersion"`: The format version that you want the response to be in. Valid values
  are: aws_v1
- `"MaxResults"`: The maximum number of results that you want returned in the response.
- `"NextToken"`: The pagination token that indicates the next set of results that you want
  to retrieve.
- `"ServiceCode"`: The code for the service whose information you want to retrieve, such as
  AmazonEC2. You can use the ServiceCode to filter the results in a GetProducts call. To
  retrieve a list of all services, leave this blank.
"""
describe_services(; aws_config::AbstractAWSConfig=global_aws_config()) = pricing("DescribeServices"; aws_config=aws_config)
describe_services(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = pricing("DescribeServices", params; aws_config=aws_config)

"""
    get_attribute_values(attribute_name, service_code)
    get_attribute_values(attribute_name, service_code, params::Dict{String,<:Any})

Returns a list of attribute values. Attibutes are similar to the details in a Price List
API offer file. For a list of available attributes, see Offer File Definitions in the AWS
Billing and Cost Management User Guide.

# Arguments
- `attribute_name`: The name of the attribute that you want to retrieve the values for,
  such as volumeType.
- `service_code`: The service code for the service whose attributes you want to retrieve.
  For example, if you want the retrieve an EC2 attribute, use AmazonEC2.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return in response.
- `"NextToken"`: The pagination token that indicates the next set of results that you want
  to retrieve.
"""
get_attribute_values(AttributeName, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = pricing("GetAttributeValues", Dict{String, Any}("AttributeName"=>AttributeName, "ServiceCode"=>ServiceCode); aws_config=aws_config)
get_attribute_values(AttributeName, ServiceCode, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = pricing("GetAttributeValues", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AttributeName"=>AttributeName, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    get_products()
    get_products(params::Dict{String,<:Any})

Returns a list of all products that match the filter criteria.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: The list of filters that limit the returned products. only products that
  match all filters are returned.
- `"FormatVersion"`: The format version that you want the response to be in. Valid values
  are: aws_v1
- `"MaxResults"`: The maximum number of results to return in the response.
- `"NextToken"`: The pagination token that indicates the next set of results that you want
  to retrieve.
- `"ServiceCode"`: The code for the service whose products you want to retrieve.
"""
get_products(; aws_config::AbstractAWSConfig=global_aws_config()) = pricing("GetProducts"; aws_config=aws_config)
get_products(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = pricing("GetProducts", params; aws_config=aws_config)


text:
for ff in names(@__MODULE__; all=true)
    fname = string(ff)
    if occursin("SetLinearSolver", fname) &&
        !occursin("#", fname) && # filter out compiler generated names
        !occursin("Dls", fname) && !occursin("Spils", fname) # filter out old names
        @eval $ff(mem, LS::SUNLinearSolver, A::Ptr, args...) = $ff(mem, LS, convert(SUNMatrix, A), args...)
    end
end

# Julia wrapper for header: sunlinsol_lapackband.h
# Automatically generated using Clang.jl

function SUNLinSol_LapackBand(y::N_Vector, A::SUNMatrix)
    ccall((:SUNLinSol_LapackBand, libsundials_sunlinsollapackband), SUNLinearSolver, (N_Vector, SUNMatrix), y, A)
end

function SUNLinSol_LapackBand(y, A)
    __y = convert(NVector, y)
    SUNLinSol_LapackBand(convert(N_Vector, __y), A)
end

function SUNLapackBand(y::N_Vector, A::SUNMatrix)
    ccall((:SUNLapackBand, libsundials_sunlinsollapackband), SUNLinearSolver, (N_Vector, SUNMatrix), y, A)
end

function SUNLapackBand(y, A)
    __y = convert(NVector, y)
    SUNLapackBand(convert(N_Vector, __y), A)
end

function SUNLinSolGetType_LapackBand(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_LapackBand, libsundials_sunlinsollapackband), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolGetID_LapackBand(S::SUNLinearSolver)
    ccall((:SUNLinSolGetID_LapackBand, libsundials_sunlinsollapackband), SUNLinearSolver_ID, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_LapackBand(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_LapackBand, libsundials_sunlinsollapackband), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetup_LapackBand(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_LapackBand, libsundials_sunlinsollapackband), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function SUNLinSolSolve_LapackBand(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_LapackBand, libsundials_sunlinsollapackband), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_LapackBand(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    SUNLinSolSolve_LapackBand(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolLastFlag_LapackBand(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_LapackBand, libsundials_sunlinsollapackband), sunindextype, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_LapackBand(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_LapackBand, libsundials_sunlinsollapackband), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_LapackBand(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_LapackBand, libsundials_sunlinsollapackband), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: sunlinsol_lapackdense.h
# Automatically generated using Clang.jl

function SUNLinSol_LapackDense(y::N_Vector, A::SUNMatrix)
    ccall((:SUNLinSol_LapackDense, libsundials_sunlinsollapackdense), SUNLinearSolver, (N_Vector, SUNMatrix), y, A)
end

function SUNLinSol_LapackDense(y, A)
    __y = convert(NVector, y)
    SUNLinSol_LapackDense(convert(N_Vector, __y), A)
end

function SUNLapackDense(y::N_Vector, A::SUNMatrix)
    ccall((:SUNLapackDense, libsundials_sunlinsollapackdense), SUNLinearSolver, (N_Vector, SUNMatrix), y, A)
end

function SUNLapackDense(y, A)
    __y = convert(NVector, y)
    SUNLapackDense(convert(N_Vector, __y), A)
end

function SUNLinSolGetType_LapackDense(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_LapackDense, libsundials_sunlinsollapackdense), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolGetID_LapackDense(S::SUNLinearSolver)
    ccall((:SUNLinSolGetID_LapackDense, libsundials_sunlinsollapackdense), SUNLinearSolver_ID, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_LapackDense(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_LapackDense, libsundials_sunlinsollapackdense), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetup_LapackDense(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_LapackDense, libsundials_sunlinsollapackdense), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function SUNLinSolSolve_LapackDense(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_LapackDense, libsundials_sunlinsollapackdense), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_LapackDense(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    SUNLinSolSolve_LapackDense(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolLastFlag_LapackDense(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_LapackDense, libsundials_sunlinsollapackdense), sunindextype, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_LapackDense(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_LapackDense, libsundials_sunlinsollapackdense), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_LapackDense(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_LapackDense, libsundials_sunlinsollapackdense), Cint, (SUNLinearSolver,), S)
end


text:
import ThreeJS

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
    d = Input(1.0)
    vbox(
        "ThreeJS Example",
        vskip(2em),
        hbox(
            vbox(
                "Sides",
                 #hbox("width",slider(1.0:5.0) >>> w),
                 #hbox("height",slider(1.0:5.0) >>> h),
                hbox("depth",slider(1.0:5.0) >>> d),
            ),
            # hskip(2em),
            # vbox(
            #     "Rotations",
            #     hbox("x", slider(0.0:5.0:360.0) >>> rx),
            #     hbox("y", slider(0.0:5.0:360.0) >>> ry),
            #     hbox("z", slider(0.0:5.0:360.0) >>> rz),
            # )
        ),
        vskip(2em),
        lift(d) do d
        ThreeJS.outerdiv() <<
            (ThreeJS.initscene() <<
                [
                 ThreeJS.mesh(0.0, 0.0, 0.0) <<
                    [
                     ThreeJS.box(1.,1., d), ThreeJS.material(Dict(:kind=>"lambert",:color=>"red"))
                    ],
                 ThreeJS.pointlight(10.0, 10.0, 10.0),
                 ThreeJS.camera(0.0, 0.0, 10.0)
                ]
            )
        end
    ) |> pad(2em)
end


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

"""
    BlockPartition(sides; neighbors=false)
    BlockPartition(side₁, side₂, ...; neighbors=false)

A method for partitioning spatial objects into blocks of given `sides`.
Optionally, compute the `neighbors` of a block as the metadata.
"""
struct BlockPartition{Dim,T} <: PartitionMethod
  sides::SVector{Dim,T}
  neighbors::Bool
end

BlockPartition(sides) = BlockPartition(sides, false)
BlockPartition(sides::NTuple; neighbors=false) =
  BlockPartition(SVector(sides), neighbors)
BlockPartition(sides::Vararg{<:Number}; neighbors=false) =
  BlockPartition(SVector(sides), neighbors)

function partition(object, method::BlockPartition)
  Dim = embeddim(object)

  psides = method.sides
  bbox = boundingbox(object)

  @assert all(psides .≤ sides(bbox)) "invalid block sides"

  # bounding box properties
  lo, up = coordinates.(extrema(bbox))
  ce = coordinates(center(bbox))

  # find number of blocks to left and right
  nleft  = @. ceil(Int, (ce - lo) / psides)
  nright = @. ceil(Int, (up - ce) / psides)

  start   = @. ce - nleft * psides
  nblocks = @. nleft + nright

  subsets   = [Int[] for i in 1:prod(nblocks)]

  # Cartesian to linear indices
  linear = LinearIndices(Dims(nblocks))

  for j in 1:nelements(object)
    coords = coordinates(centroid(object, j))

    # find block coordinates
    c = @. floor(Int, (coords - start) / psides) + 1
    c = @. clamp(c, 1, nblocks)
    bcoords = CartesianIndex(Tuple(c))

    # block index
    i = linear[bcoords]

    append!(subsets[i], j)
  end

  # intitialize metadata
  metadata = Dict()
  neighbors = [Int[] for i in 1:prod(nblocks)]

  # neighboring blocks metadata
  if method.neighbors == true
    bstart  = CartesianIndex(ntuple(i -> 1, Dim))
    boffset = CartesianIndex(ntuple(i -> 1, Dim))
    bfinish = CartesianIndex(Dims(nblocks))
    for (i, bcoords) in enumerate(bstart:bfinish)
      for b in (bcoords - boffset):(bcoords + boffset)
        if all(Tuple(bstart) .≤ Tuple(b) .≤ Tuple(bfinish)) && b ≠ bcoords
          push!(neighbors[i], linear[b])
        end
      end
    end

    metadata[:neighbors] = neighbors
  end

  # filter out empty blocks
  empty = isempty.(subsets)
  subsets = subsets[.!empty]
  neighbors = neighbors[.!empty]
  for i in findall(empty)
    for n in neighbors
      setdiff!(n, i)
    end
  end

  Partition(object, subsets, metadata)
end