text:
using ACEbonds, ACE, ACEbase, Test, StaticArrays, LinearAlgebra, JuLIP
using ACEbase.Testing
using ACEbonds: CylindricalCutoff
using ACE: discrete_jacobi, Rn1pBasis, scal1pbasis, Scal1pBasis,
           evaluate, evaluate_d, Trig1pBasis, @λ,
           Product1pBasis, Categorical1pBasis, SimpleSparseBasis,
           SymmetricBasis, Invariant, PIBasis

##

r0cut = 3.2
rcut = 4.0
zcut = 5.0

cutoff = CylindricalCutoff(r0cut, rcut, zcut)

maxdeg = 8
maxL = 5

Jz = discrete_jacobi(maxdeg; pcut = 2, pin = 2, xcut = zcut, xin = -zcut)
Jr = discrete_jacobi(maxdeg; pcut = 2, pin = 2, xcut = rcut, xin = -rcut)
Jr0 = discrete_jacobi(maxdeg; pcut = 2, pin = 2, xcut = r0cut, xin = -r0cut)

Cbe = Categorical1pBasis([:bond, :env], :be, :be, "Cbe")
Zm = Scal1pBasis(:z, nothing, :m, Jz, "Zm")
Rn = Scal1pBasis(:r, nothing, :n, Jr, "Rn")
El = Trig1pBasis(maxL; varsym = :θ, lsym = :l, label = "El")

Pk = Scal1pBasis(:rij, nothing, :k, Jr0, "Jk")

##
B1p = Product1pBasis((Cbe, Pk, Rn, El, Zm))
Bsel = ACEbonds.SparseBondBasis(; maxorder = 3,
         default_maxdeg = maxdeg,
         weight = Dict{Symbol, Float64}(:m => 1.0, :n => 1.0, :k => 1.0, :l => 1.0),
                                 )
# ACE.init1pspec!(B1p, Bsel)
# length(B1p)

basis = SymmetricBasis(ACE.Invariant(), B1p, ACE.NoSym(), Bsel; isreal=true)
@show length(basis)

model = ACE.LinearACEModel(basis)
θ = ACE.params(model)
θ = randn(length(θ)) ./ (1:length(θ)).^2
ACE.set_params!(model, θ)

##

using ACEbonds: ACEBondPotentialBasis, ACEBondPotential

zSi = AtomicNumber(:Si)
_bases = Dict((zSi, zSi) => basis)
_models = Dict((zSi, zSi) => model)
inds = Dict((zSi, zSi) => 1:length(basis))
pot = ACEBondPotential(_models, cutoff)
potbasis = ACEbonds.basis(pot)
# potbasis = ACEBondPotentialBasis(models, inds, cutoff)

at = rattle!(set_pbc!(bulk(:Si, cubic=true) * (2, 2, 1), false), 0.2)

##

@info("Testing energy pot vs energy basis")
println_slim(@test energy(pot, at) ≈ dot(energy(potbasis, at), θ))

@info("Testing forces pot vs basis")
dB = forces(potbasis, at)
println_slim(@test( sum( θ[n] * dB[n, :] for n = 1:length(θ) ) ≈ forces(pot, at) ))

@info("Finite-difference forces test")
println_slim(@test JuLIP.Testing.fdtest(pot, at))

##

@info("Finite-difference virial test")
at = set_pbc!(bulk(:Si, cubic=true) * (2, 1, 1), true)
JuLIP.set_variablecell!(at, true)
println_slim(@test JuLIP.Testing.fdtest(pot, at))

at = rattle!(set_pbc!(bulk(:Si, cubic=true) * (2, 1, 1), true), 0.2)
JuLIP.set_variablecell!(at, true)
println_slim(@test JuLIP.Testing.fdtest(pot, at))

@info("Virial pt vs basis")
at = rattle!(set_pbc!(bulk(:Si, cubic=true) * (2, 2, 2), true), 0.2)
vB = virial(potbasis, at)
println_slim(@test( sum( θ[n] * vB[n] for n = 1:length(potbasis) ) ≈ virial(pot, at) ))


text:
function write_output(	file_path::String,
					particle_type::String,
					R::Array{Float64, 2},
					Lx::Float64,
					Ly::Float64,
					Lz::Float64,
					phi::Float64,
					X::Array{Float64, 1},
					Y::Array{Float64, 1},
					Z::Array{Float64, 1},
					Q0::Array{Float64, 1},
					Q1::Array{Float64, 1},
					Q2::Array{Float64, 1},
					Q3::Array{Float64, 1},
					t_exec::Float64)

	file_stream::IOStream = open(file_path, "w")

	@printf(file_stream, "%s", "<output_generation>\n")

	write_key(file_stream, "particle_type", particle_type)
	write_key(file_stream, "R", R[:])
	write_key(file_stream, "domain_size_x", Lx)
	write_key(file_stream, "domain_size_y", Ly)
	write_key(file_stream, "domain_size_z", Lz)
	write_key(file_stream, "phi", phi)
	write_key(file_stream, "X", X)
	write_key(file_stream, "Y", Y)
	write_key(file_stream, "Z", Z)
	write_key(file_stream, "Q0", Q0)
	write_key(file_stream, "Q1", Q1)
	write_key(file_stream, "Q2", Q2)
	write_key(file_stream, "Q3", Q3)
	write_key(file_stream, "execution_time", t_exec)

	@printf(file_stream, "%s", "</output_generation>")

	close(file_stream)

	nothing
end

ext: jl

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
"""
    struct Model{F,argnames,defaultnames,missings,Targs,Tdefaults}
        name::Symbol
        f::F
        args::NamedTuple{argnames,Targs}
        defaults::NamedTuple{defaultnames,Tdefaults}
    end

A `Model` struct with model evaluation function of type `F`, arguments of names `argnames`
types `Targs`, default arguments of names `defaultnames` with types `Tdefaults`, and missing
arguments `missings`.

Here `argnames`, `defaultargnames`, and `missings` are tuples of symbols, e.g. `(:a, :b)`.

An argument with a type of `Missing` will be in `missings` by default. However, in
non-traditional use-cases `missings` can be defined differently. All variables in `missings`
are treated as random variables rather than observations.

The default arguments are used internally when constructing instances of the same model with
different arguments.

# Examples

```julia
julia> Model(f, (x = 1.0, y = 2.0))
Model{typeof(f),(:x, :y),(),(),Tuple{Float64,Float64},Tuple{}}(f, (x = 1.0, y = 2.0), NamedTuple())

julia> Model(f, (x = 1.0, y = 2.0), (x = 42,))
Model{typeof(f),(:x, :y),(:x,),(),Tuple{Float64,Float64},Tuple{Int64}}(f, (x = 1.0, y = 2.0), (x = 42,))

julia> Model{(:y,)}(f, (x = 1.0, y = 2.0), (x = 42,)) # with special definition of missings
Model{typeof(f),(:x, :y),(:x,),(:y,),Tuple{Float64,Float64},Tuple{Int64}}(f, (x = 1.0, y = 2.0), (x = 42,))
```
"""
struct Model{F,argnames,defaultnames,missings,Targs,Tdefaults} <:
       AbstractProbabilisticProgram
    name::Symbol
    f::F
    args::NamedTuple{argnames,Targs}
    defaults::NamedTuple{defaultnames,Tdefaults}

    @doc """
        Model{missings}(name::Symbol, f, args::NamedTuple, defaults::NamedTuple)

    Create a model of name `name` with evaluation function `f` and missing arguments
    overwritten by `missings`.
    """
    function Model{missings}(
        name::Symbol,
        f::F,
        args::NamedTuple{argnames,Targs},
        defaults::NamedTuple{defaultnames,Tdefaults},
    ) where {missings,F,argnames,Targs,defaultnames,Tdefaults}
        return new{F,argnames,defaultnames,missings,Targs,Tdefaults}(
            name, f, args, defaults
        )
    end
end

"""
    Model(name::Symbol, f, args::NamedTuple[, defaults::NamedTuple = ()])

Create a model of name `name` with evaluation function `f` and missing arguments deduced
from `args`.

Default arguments `defaults` are used internally when constructing instances of the same
model with different arguments.
"""
@generated function Model(
    name::Symbol, f::F, args::NamedTuple{argnames,Targs}, defaults::NamedTuple=NamedTuple()
) where {F,argnames,Targs}
    missings = Tuple(name for (name, typ) in zip(argnames, Targs.types) if typ <: Missing)
    return :(Model{$missings}(name, f, args, defaults))
end

"""
    (model::Model)([rng, varinfo, sampler, context])

Sample from the `model` using the `sampler` with random number generator `rng` and the
`context`, and store the sample and log joint probability in `varinfo`.

The method resets the log joint probability of `varinfo` and increases the evaluation
number of `sampler`.
"""
function (model::Model)(
    rng::Random.AbstractRNG,
    varinfo::AbstractVarInfo=VarInfo(),
    sampler::AbstractSampler=SampleFromPrior(),
    context::AbstractContext=DefaultContext(),
)
    return model(varinfo, SamplingContext(rng, sampler, context))
end

(model::Model)(context::AbstractContext) = model(VarInfo(), context)
function (model::Model)(varinfo::AbstractVarInfo, context::AbstractContext)
    if Threads.nthreads() == 1
        return evaluate_threadunsafe(model, varinfo, context)
    else
        return evaluate_threadsafe(model, varinfo, context)
    end
end

function (model::Model)(args...)
    return model(Random.GLOBAL_RNG, args...)
end

# without VarInfo
function (model::Model)(rng::Random.AbstractRNG, sampler::AbstractSampler, args...)
    return model(rng, VarInfo(), sampler, args...)
end

# without VarInfo and without AbstractSampler
function (model::Model)(rng::Random.AbstractRNG, context::AbstractContext)
    return model(rng, VarInfo(), SampleFromPrior(), context)
end

"""
    evaluate_threadunsafe(model, varinfo, context)

Evaluate the `model` without wrapping `varinfo` inside a `ThreadSafeVarInfo`.

If the `model` makes use of Julia's multithreading this will lead to undefined behaviour.
This method is not exposed and supposed to be used only internally in DynamicPPL.

See also: [`evaluate_threadsafe`](@ref)
"""
function evaluate_threadunsafe(model, varinfo, context)
    resetlogp!(varinfo)
    return _evaluate(model, varinfo, context)
end

"""
    evaluate_threadsafe(model, varinfo, context)

Evaluate the `model` with `varinfo` wrapped inside a `ThreadSafeVarInfo`.

With the wrapper, Julia's multithreading can be used for observe statements in the `model`
but parallel sampling will lead to undefined behaviour.
This method is not exposed and supposed to be used only internally in DynamicPPL.

See also: [`evaluate_threadunsafe`](@ref)
"""
function evaluate_threadsafe(model, varinfo, context)
    resetlogp!(varinfo)
    wrapper = ThreadSafeVarInfo(varinfo)
    result = _evaluate(model, wrapper, context)
    setlogp!(varinfo, getlogp(wrapper))
    return result
end

"""
    _evaluate(model::Model, varinfo, context)

Evaluate the `model` with the arguments matching the given `context` and `varinfo` object.
"""
@generated function _evaluate(
    model::Model{_F,argnames}, varinfo, context
) where {_F,argnames}
    unwrap_args = [:($matchingvalue(context, varinfo, model.args.$var)) for var in argnames]
    return :(model.f(model, varinfo, context, $(unwrap_args...)))
end

"""
    getargnames(model::Model)

Get a tuple of the argument names of the `model`.
"""
getargnames(model::Model{_F,argnames}) where {argnames,_F} = argnames

"""
    getmissings(model::Model)

Get a tuple of the names of the missing arguments of the `model`.
"""
getmissings(model::Model{_F,_a,_d,missings}) where {missings,_F,_a,_d} = missings

"""
    nameof(model::Model)

Get the name of the `model` as `Symbol`.
"""
Base.nameof(model::Model) = model.name

"""
    logjoint(model::Model, varinfo::AbstractVarInfo)

Return the log joint probability of variables `varinfo` for the probabilistic `model`.

See [`logjoint`](@ref) and [`loglikelihood`](@ref).
"""
function logjoint(model::Model, varinfo::AbstractVarInfo)
    model(varinfo, DefaultContext())
    return getlogp(varinfo)
end

"""
    logprior(model::Model, varinfo::AbstractVarInfo)

Return the log prior probability of variables `varinfo` for the probabilistic `model`.

See also [`logjoint`](@ref) and [`loglikelihood`](@ref).
"""
function logprior(model::Model, varinfo::AbstractVarInfo)
    model(varinfo, PriorContext())
    return getlogp(varinfo)
end

"""
    loglikelihood(model::Model, varinfo::AbstractVarInfo)

Return the log likelihood of variables `varinfo` for the probabilistic `model`.

See also [`logjoint`](@ref) and [`logprior`](@ref).
"""
function Distributions.loglikelihood(model::Model, varinfo::AbstractVarInfo)
    model(varinfo, LikelihoodContext())
    return getlogp(varinfo)
end

"""
    generated_quantities(model::Model, chain::AbstractChains)

Execute `model` for each of the samples in `chain` and return an array of the values
returned by the `model` for each sample.

# Examples
## General
Often you might have additional quantities computed inside the model that you want to
inspect, e.g.
```julia
@model function demo(x)
    # sample and observe
    θ ~ Prior()
    x ~ Likelihood()
    return interesting_quantity(θ, x)
end
m = demo(data)
chain = sample(m, alg, n)
# To inspect the `interesting_quantity(θ, x)` where `θ` is replaced by samples
# from the posterior/`chain`:
generated_quantities(m, chain) # <= results in a `Vector` of returned values
                               #    from `interesting_quantity(θ, x)`
```
## Concrete (and simple)
```julia
julia> using DynamicPPL, Turing

julia> @model function demo(xs)
           s ~ InverseGamma(2, 3)
           m_shifted ~ Normal(10, √s)
           m = m_shifted - 10

           for i in eachindex(xs)
               xs[i] ~ Normal(m, √s)
           end

           return (m, )
       end
demo (generic function with 1 method)

julia> model = demo(randn(10));

julia> chain = sample(model, MH(), 10);

julia> generated_quantities(model, chain)
10×1 Array{Tuple{Float64},2}:
 (2.1964758025119338,)
 (2.1964758025119338,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.043088571494005024,)
 (-0.16489786710222099,)
 (-0.16489786710222099,)
```
"""
function generated_quantities(model::Model, chain::AbstractChains)
    varinfo = VarInfo(model)
    iters = Iterators.product(1:size(chain, 1), 1:size(chain, 3))
    return map(iters) do (sample_idx, chain_idx)
        setval_and_resample!(varinfo, chain, sample_idx, chain_idx)
        model(varinfo)
    end
end


text:
using ChainComplexes
using Test

@testset "ChainComplexes.jl" begin
    true
end

max_stars_repo_name: osumray/ChainComplexes

max_issues_repo_name: osumray/ChainComplexes

max_forks_repo_name: osumray/ChainComplexes


text:
using BinDeps

@BinDeps.setup

version = "0.9.15"
url = "https://github.com/LMDB/lmdb/archive/LMDB_$(version).tar.gz"

liblmdb = library_dependency("liblmdb")

provides(Sources, URI(url), liblmdb, unpacked_dir="lmdb-LMDB_$version")

lmdbsrcdir = joinpath(BinDeps.srcdir(liblmdb),"lmdb-LMDB_$version","libraries","liblmdb")
lmdblibfile = joinpath(BinDeps.libdir(liblmdb),liblmdb.name*"."*BinDeps.shlib_ext)

provides(BuildProcess,
	(@build_steps begin
		GetSources(liblmdb)
		CreateDirectory(BinDeps.libdir(liblmdb))
		FileRule(lmdblibfile, @build_steps begin
			MakeTargets(lmdbsrcdir, ["liblmdb.so"])
			`mv $lmdbsrcdir/liblmdb.so $lmdblibfile`
		end)
	end), liblmdb)

@BinDeps.install Dict( :liblmdb => :liblmdb)


text:
## Residual Constrained Alternating Minimization
# Function Definitions

function AltMinBackslash(b,r,K)
# Performs RCAM on a single matricised monochromatic slice
# Calls: minFrob

    # Generate L0
    m,n = size(b)
    L = randn(m,r)

    local R

    # Begin Alterations
    for k = 1:K

        # Solve for R
        R = minFrob(L,b)

        # Solve for L
        L = minFrob(R,b')

    end

    return L, R
end
export AltMinBackslash

function minFrob(L,b)
# Solves for each row of input
# Calls: minFrobRsub

    m = size(b,2)
    r = size(L,2)

    # Initialize R
    R = zeros(m,r)
    if ~isreal(b)
        R = complex(R)
    end

    # Solve for each row of R
    for k=1:m

        R[k,:] = minFrobRsub(L,b,k)

    end

    R = conj(R)

    return R
end

function minFrobRsub(L,b,k)
# Solves a single row

    # Solve for the kth row of R
    bsub = b[:,k]
    ind = find(x -> x != 0, bsub)
    bsub = bsub[ind]
    L = L[ind,:]

    v = L\bsub

    return v

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
using Xpress

@testset "Thermal loads" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/thermal_load.json")

    @test round(results["ExistingBoiler"]["year_one_boiler_fuel_consumption_mmbtu"], digits=0) ≈ 2905

    data = JSON.parsefile("./scenarios/thermal_load.json")
    data["DomesticHotWaterLoad"]["fuel_loads_mmbtu_per_hour"] = repeat([0.5], 8760)
    data["SpaceHeatingLoad"]["fuel_loads_mmbtu_per_hour"] = repeat([0.5], 8760)
    s = Scenario(data)
    inputs = REoptInputs(s)
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, inputs)

    @test round(results["ExistingBoiler"]["year_one_boiler_fuel_consumption_mmbtu"], digits=0) ≈ 8760
end

@testset "CHP Sizing" begin
    # Sizing CHP with non-constant efficiency, no cost curve, no unavailability_periods
    data_sizing = JSON.parsefile("./scenarios/chp_sizing.json")
    s = Scenario(data_sizing)
    inputs = REoptInputs(s)
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "MIPRELSTOP" => 0.01, "OUTPUTLOG" => 0))
    results = run_reopt(m, inputs)

    @test round(results["CHP"]["size_kw"], digits=0) ≈ 468.7 atol=1.0
    @test round(results["Financial"]["lcc"], digits=0) ≈ 1.3476e7 atol=1.0e7
end

@testset "CHP Cost Curve and Min Allowable Size" begin
    # Fixed size CHP with cost curve, no unavailability_periods
    data_cost_curve = JSON.parsefile("./scenarios/chp_sizing.json")
    data_cost_curve["CHP"] = Dict()
    data_cost_curve["CHP"]["prime_mover"] = "recip_engine"
    data_cost_curve["CHP"]["size_class"] = 2
    data_cost_curve["CHP"]["fuel_cost_per_mmbtu"] = 8.0
    data_cost_curve["CHP"]["min_kw"] = 0
    data_cost_curve["CHP"]["min_allowable_kw"] = 555.5
    data_cost_curve["CHP"]["max_kw"] = 1000
    data_cost_curve["CHP"]["installed_cost_per_kw"] = 1800.0
    data_cost_curve["CHP"]["installed_cost_per_kw"] = [2300.0, 1800.0, 1500.0]
    data_cost_curve["CHP"]["tech_sizes_for_cost_curve"] = [100.0, 300.0, 1140.0]

    data_cost_curve["CHP"]["federal_itc_pct"] = 0.1
    data_cost_curve["CHP"]["macrs_option_years"] = 0
    data_cost_curve["CHP"]["macrs_bonus_pct"] = 0.0
    data_cost_curve["CHP"]["macrs_itc_reduction"] = 0.0

    expected_x = data_cost_curve["CHP"]["min_allowable_kw"]
    cap_cost_y = data_cost_curve["CHP"]["installed_cost_per_kw"]
    cap_cost_x = data_cost_curve["CHP"]["tech_sizes_for_cost_curve"]
    slope = (cap_cost_x[3] * cap_cost_y[3] - cap_cost_x[2] * cap_cost_y[2]) / (cap_cost_x[3] - cap_cost_x[2])
    init_capex_chp_expected = cap_cost_x[2] * cap_cost_y[2] + (expected_x - cap_cost_x[2]) * slope
    lifecycle_capex_chp_expected = init_capex_chp_expected -
        REopt.npv(data_cost_curve["Financial"]["offtaker_discount_pct"],
        [0, init_capex_chp_expected * data_cost_curve["CHP"]["federal_itc_pct"]])

    #PV
    data_cost_curve["PV"]["min_kw"] = 1500
    data_cost_curve["PV"]["max_kw"] = 1500
    data_cost_curve["PV"]["installed_cost_per_kw"] = 1600
    data_cost_curve["PV"]["federal_itc_pct"] = 0.26
    data_cost_curve["PV"]["macrs_option_years"] = 0
    data_cost_curve["PV"]["macrs_bonus_pct"] = 0.0
    data_cost_curve["PV"]["macrs_itc_reduction"] = 0.0

    init_capex_pv_expected = data_cost_curve["PV"]["max_kw"] * data_cost_curve["PV"]["installed_cost_per_kw"]
    lifecycle_capex_pv_expected = init_capex_pv_expected -
        REopt.npv(data_cost_curve["Financial"]["offtaker_discount_pct"],
        [0, init_capex_pv_expected * data_cost_curve["PV"]["federal_itc_pct"]])

    s = Scenario(data_cost_curve)
    inputs = REoptInputs(s)
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "MIPRELSTOP" => 0.01, "OUTPUTLOG" => 0))
    results = run_reopt(m, inputs)

    init_capex_total_expected = init_capex_chp_expected + init_capex_pv_expected
    lifecycle_capex_total_expected = lifecycle_capex_chp_expected + lifecycle_capex_pv_expected

    init_capex_total = results["Financial"]["initial_capital_costs"]
    lifecycle_capex_total = results["Financial"]["initial_capital_costs_after_incentives"]

    # Check initial CapEx (pre-incentive/tax) and life cycle CapEx (post-incentive/tax) cost with expect
    @test init_capex_total_expected ≈ init_capex_total atol=0.0001*init_capex_total_expected
    @test lifecycle_capex_total_expected ≈ lifecycle_capex_total atol=0.0001*lifecycle_capex_total_expected

    # Test CHP.min_allowable_kw - the size would otherwise be ~100 kW less by setting min_allowable_kw to zero
    @test results["CHP"]["size_kw"] ≈ data_cost_curve["CHP"]["min_allowable_kw"] atol=0.1
end

@testset "CHP Unavailability and Outage" begin
    """
    Validation to ensure that:
        1) CHP meets load during outage without exporting
        2) CHP never exports if chp.can_wholesale and chp.can_net_meter inputs are False (default)
        3) CHP does not "curtail", i.e. send power to a load bank when chp.can_curtail is False (default)
        4) CHP min_turn_down_pct is ignored during an outage
        5) **Not until cooling is added:** Cooling load gets zeroed out during the outage period
        6) Unavailability intervals that intersect with grid-outages get ignored
        7) Unavailability intervals that do not intersect with grid-outages result in no CHP production
    """
    # Sizing CHP with non-constant efficiency, no cost curve, no unavailability_periods
    data = JSON.parsefile("./scenarios/chp_unavailability_outage.json")

    # Add unavailability periods that 1) intersect (ignored) and 2) don't intersect with outage period
    data["CHP"]["unavailability_periods"] = [Dict([("month", 1), ("start_week_of_month", 2),
            ("start_day_of_week", 1), ("start_hour", 1), ("duration_hours", 8)]),
            Dict([("month", 1), ("start_week_of_month", 2),
            ("start_day_of_week", 3), ("start_hour", 9), ("duration_hours", 8)])]

    # Manually doing the math from the unavailability defined above
    unavail_1_start = 24 + 1
    unavail_1_end = unavail_1_start + 8 - 1
    unavail_2_start = 24*3 + 9
    unavail_2_end = unavail_2_start + 8 - 1

    # Specify the CHP.min_turn_down_pct which is NOT used during an outage
    data["CHP"]["min_turn_down_pct"] = 0.5
    # Specify outage period; outage timesteps are 1-indexed
    outage_start = unavail_1_start
    data["ElectricUtility"]["outage_start_time_step"] = outage_start
    outage_end = unavail_1_end
    data["ElectricUtility"]["outage_end_time_step"] = outage_end
    data["ElectricLoad"]["critical_load_pct"] = 0.25

    s = Scenario(data)
    inputs = REoptInputs(s)
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "MIPRELSTOP" => 0.01, "OUTPUTLOG" => 0))
    results = run_reopt(m, inputs)

    tot_elec_load = results["ElectricLoad"]["load_series_kw"]
    chp_total_elec_prod = results["CHP"]["year_one_electric_production_series_kw"]
    chp_to_load = results["CHP"]["year_one_to_load_series_kw"]
    chp_export = results["CHP"]["year_one_to_grid_series_kw"]
    #cooling_elec_load = results["LoadProfileChillerThermal"]["year_one_chiller_electric_load_kw"]

    # The values compared to the expected values
    #@test sum([(chp_to_load[i] - tot_elec_load[i]) for i in outage_start:outage_end])) == 0.0
    critical_load = tot_elec_load[outage_start:outage_end] * data["ElectricLoad"]["critical_load_pct"]
    @test sum(chp_to_load[outage_start:outage_end]) ≈ sum(critical_load) atol=0.1
    @test sum(chp_export) == 0.0
    @test sum(chp_total_elec_prod) ≈ sum(chp_to_load) atol=1.0e-5*sum(chp_total_elec_prod)
    #@test sum(cooling_elec_load[outage_start:outage_end]) == 0.0
    @test sum(chp_total_elec_prod[unavail_2_start:unavail_2_end]) == 0.0
end

@testset "CHP Supplementary firing and standby" begin
    """
    Test to ensure that supplementary firing and standby charges work as intended.  The thermal and
    electrical loads are constant, and the CHP system size is fixed; the supplementary firing has a
    similar cost to the boiler and is purcahsed and used when the boiler efficiency is set to a lower
    value than that of the supplementary firing. The test also ensures that demand charges are
    correctly calculated when CHP is and is not allowed to reduce demand charges.
    """
    data = JSON.parsefile("./scenarios/chp_supplementary_firing.json")
    data["CHP"]["supplementary_firing_capital_cost_per_kw"] = 10000
    data["ElectricLoad"]["loads_kw"] = repeat([800.0], 8760)
    data["DomesticHotWaterLoad"]["fuel_loads_mmbtu_per_hour"] = repeat([6.0], 8760)
    data["SpaceHeatingLoad"]["fuel_loads_mmbtu_per_hour"] = repeat([6.0], 8760)
    #part 1: supplementary firing not used when less efficient than the boiler and expensive
    m1 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    s = Scenario(data)
    inputs = REoptInputs(s)
    results = run_reopt(m1, inputs)
    @test results["CHP"]["size_kw"] == 800
    @test results["CHP"]["size_supplemental_firing_kw"] == 0
    @test results["CHP"]["year_one_electric_energy_produced_kwh"] ≈ 800*8760 rtol=1e-5
    @test results["CHP"]["year_one_thermal_energy_produced_mmbtu"] ≈ 800*(0.4418/0.3573)*8760/293.07107 rtol=1e-5
    @test results["ElectricTariff"]["lifecycle_demand_cost"] == 0

    #part 2: supplementary firing used when more efficient than the boiler and low-cost; demand charges not reduced by CHP
    data["CHP"]["supplementary_firing_capital_cost_per_kw"] = 10
    data["CHP"]["reduces_demand_charges"] = false
    data["ExistingBoiler"]["efficiency"] = 0.85
    m2 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    s = Scenario(data)
    inputs = REoptInputs(s)
    results = run_reopt(m2, inputs)
    @test results["CHP"]["size_supplemental_firing_kw"] ≈ 278.73 atol=0.1
    @test results["CHP"]["year_one_thermal_energy_produced_mmbtu"] ≈ 138624 rtol=1e-5
    @test results["ElectricTariff"]["lifecycle_demand_cost"] ≈ 5212.7 rtol=1e-5
end

#=
add a time-of-export rate that is greater than retail rate for the month of January,
check to make sure that PV does NOT export unless the site load is met first for the month of January.
=#
@testset "Do not allow_simultaneous_export_import" begin
    model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    data = JSON.parsefile("./scenarios/monthly_rate.json")

    # create wholesale_rate with compensation in January > retail rate
    jan_rate = data["ElectricTariff"]["monthly_energy_rates"][1]
    data["ElectricTariff"]["wholesale_rate"] =
        append!(repeat([jan_rate + 0.1], 31 * 24), repeat([0.0], 8760 - 31*24))
    data["ElectricTariff"]["monthly_demand_rates"] = repeat([0], 12)
    data["ElectricUtility"] = Dict("allow_simultaneous_export_import" => false)

    s = Scenario(data)
    inputs = REoptInputs(s)
    results = run_reopt(model, inputs)

    @test all(x == 0.0 for (i,x) in enumerate(results["ElectricUtility"]["year_one_to_load_series_kw"][1:744])
              if results["PV"]["year_one_to_grid_series_kw"][i] > 0)
end

@testset "Solar and Storage w/BAU" begin
    m1 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    m2 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt([m1,m2], "./scenarios/pv_storage.json")

    @test results["PV"]["size_kw"] ≈ 216.6667 atol=0.01
    @test results["PV"]["lcoe_per_kwh"] ≈ 0.0483 atol = 0.001
    @test results["Financial"]["lcc"] ≈ 1.240037e7 rtol=1e-5
    @test results["Financial"]["lcc_bau"] ≈ 12766397 rtol=1e-5
    @test results["Storage"]["size_kw"] ≈ 55.9 atol=0.1
    @test results["Storage"]["size_kwh"] ≈ 78.9 atol=0.1
    proforma_npv = REopt.npv(results["Financial"]["offtaker_annual_free_cashflows"] -
        results["Financial"]["offtaker_annual_free_cashflows_bau"], 0.081)
    @test results["Financial"]["npv"] ≈ proforma_npv rtol=0.0001
end

@testset "Outage with Generator, outate simulator, BAU critical load outputs" begin
    m1 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    m2 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    p = REoptInputs("./scenarios/generator.json")
    results = run_reopt([m1,m2], p)
    @test results["Generator"]["size_kw"] ≈ 8.13 atol=0.01
    @test (sum(results["Generator"]["year_one_to_load_series_kw"][i] for i in 1:9) +
           sum(results["Generator"]["year_one_to_load_series_kw"][i] for i in 13:8760)) == 0
    @test results["ElectricLoad"]["bau_critical_load_met"] == false
    @test results["ElectricLoad"]["bau_critical_load_met_time_steps"] == 0

    simresults = simulate_outages(results, p)
    @test simresults["resilience_hours_max"] == 11
end

@testset "Minimize Unserved Load" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/outage.json")

    @test results["Outages"]["expected_outage_cost"] ≈ 0
    @test sum(results["Outages"]["unserved_load_per_outage"]) ≈ 0
    @test value(m[:binMGTechUsed]["Generator"]) == 1
    @test value(m[:binMGTechUsed]["PV"]) == 0
    @test value(m[:binMGStorageUsed]) == 1
    @test results["Financial"]["lcc"] ≈ 7.3879557e7 atol=5e4

    #=
    Scenario with $0/kWh value_of_lost_load_per_kwh, 12x169 hour outages, 1kW load/hour, and min_resil_timesteps = 168
    - should meet 168 kWh in each outage such that the total unserved load is 12 kWh
    =#
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/nogridcost_minresilhours.json")
    @test sum(results["Outages"]["unserved_load_per_outage"]) ≈ 12

    # testing dvUnserved load, which would output 100 kWh for this scenario before output fix
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/nogridcost_multiscenario.json")
    @test sum(results["Outages"]["unserved_load_per_outage"]) ≈ 60

end

@testset "Multiple Sites" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    ps = [
        REoptInputs("./scenarios/pv_storage.json"),
        REoptInputs("./scenarios/monthly_rate.json"),
    ];
    results = run_reopt(m, ps)
    @test results[3]["Financial"]["lcc"] + results[10]["Financial"]["lcc"] ≈ 1.240037e7 + 437169.0 rtol=1e-5
end

@testset "MPC" begin
    model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    r = run_mpc(model, "./scenarios/mpc.json")
    @test maximum(r["ElectricUtility"]["to_load_series_kw"][1:15]) <= 98.0
    @test maximum(r["ElectricUtility"]["to_load_series_kw"][16:24]) <= 97.0
    @test sum(r["PV"]["to_grid_series_kw"]) ≈ 0
end

@testset "Complex Incentives" begin
    """
    This test was compared against the API test:
        reo.tests.test_reopt_url.EntryResourceTest.test_complex_incentives
    when using the hardcoded levelization_factor in this package's REoptInputs function.
    The two LCC's matched within 0.00005%. (The Julia pkg LCC is  1.0971991e7)
    """
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/incentives.json")
    @test results["Financial"]["lcc"] ≈ 1.094596365e7 atol=5e4
end

@testset verbose = true "Rate Structures" begin

    @testset "Tiered Energy" begin
        m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(m, "./scenarios/tiered_rate.json")
        @test results["ElectricTariff"]["year_one_energy_cost"] ≈ 2342.88
    end

    @testset "Lookback Demand Charges" begin
        m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(m, "./scenarios/lookback_rate.json")
        @test results["ElectricTariff"]["year_one_demand_cost"] ≈ 721.99
    end

    @testset "Blended tariff" begin
        model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(model, "./scenarios/no_techs.json")
        @test results["ElectricTariff"]["year_one_energy_cost"] ≈ 1000.0
        @test results["ElectricTariff"]["year_one_demand_cost"] ≈ 136.99
    end

    @testset "Coincident Peak Charges" begin
        model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(model, "./scenarios/coincident_peak.json")
        @test results["ElectricTariff"]["year_one_coincident_peak_cost"] ≈ 15.0
        @test results["ElectricTariff"]["lifecycle_coincident_peak_cost"] ≈ 15.0 * 12.94887 atol=0.1
    end

    @testset "URDB sell rate" begin
        #= The URDB contains at least one "Customer generation" tariff that only has a "sell" key in the energyratestructure (the tariff tested here)
        =#
        model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        p = REoptInputs("./scenarios/URDB_customer_generation.json")
        results = run_reopt(model, p)
        @test results["PV"]["size_kw"] ≈ p.max_sizes["PV"]
    end

    # # tiered monthly demand rate  TODO: expected results?
    # m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    # data = JSON.parsefile("./scenarios/tiered_rate.json")
    # data["ElectricTariff"]["urdb_label"] = "59bc22705457a3372642da67"
    # s = Scenario(data)
    # inputs = REoptInputs(s)
    # results = run_reopt(m, inputs)

    # TODO test for tiered TOU demand rates
end

@testset "Wind" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/wind.json")
    @test results["Wind"]["size_kw"] ≈ 3752 atol=0.1
    @test results["Financial"]["lcc"] ≈ 8.591017e6 rtol=1e-5
    #=
    0.5% higher LCC in this package as compared to API ? 8,591,017 vs 8,551,172
    - both have zero curtailment
    - same energy to grid: 5,839,317 vs 5,839,322
    - same energy to load: 4,160,683 vs 4,160,677
    - same city: Boulder
    - same total wind prod factor

    REopt.jl has:
    - bigger turbine: 3752 vs 3735
    - net_capital_costs_plus_om: 8,576,590 vs. 8,537,480

    TODO: will these discrepancies be addressed once NMIL binaries are added?
    =#
end

@testset "Multiple PVs" begin
    m1 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    m2 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt([m1,m2], "./scenarios/multiple_pvs.json")

    ground_pv = results["PV"][findfirst(pv -> pv["name"] == "ground", results["PV"])]
    roof_west = results["PV"][findfirst(pv -> pv["name"] == "roof_west", results["PV"])]
    roof_east = results["PV"][findfirst(pv -> pv["name"] == "roof_east", results["PV"])]

    @test ground_pv["size_kw"] ≈ 15 atol=0.1
    @test roof_west["size_kw"] ≈ 7 atol=0.1
    @test roof_east["size_kw"] ≈ 4 atol=0.1
    @test ground_pv["lifecycle_om_cost_bau"] ≈ 782.0 atol=0.1
    @test roof_west["lifecycle_om_cost_bau"] ≈ 782.0 atol=0.1
    @test ground_pv["average_annual_energy_produced_kwh_bau"] ≈ 8844.19 atol=0.1
    @test roof_west["average_annual_energy_produced_kwh_bau"] ≈ 7440.1 atol=0.1
    @test ground_pv["average_annual_energy_produced_kwh"] ≈ 26533.54 atol=0.1
    @test roof_west["average_annual_energy_produced_kwh"] ≈ 10416.52 atol=0.1
    @test roof_east["average_annual_energy_produced_kwh"] ≈ 6482.37 atol=0.1
end

## equivalent REopt API Post for test 2:
#   NOTE have to hack in API levelization_factor to get LCC within 5e-5 (Mosel tol)
# {"Scenario": {
#     "Site": {
#         "longitude": -118.1164613,
#         "latitude": 34.5794343,
#         "roof_squarefeet": 5000.0,
#         "land_acres": 1.0,
#     "PV": {
#         "macrs_bonus_pct": 0.4,
#         "installed_cost_per_kw": 2000.0,
#         "tilt": 34.579,
#         "degradation_pct": 0.005,
#         "macrs_option_years": 5,
#         "federal_itc_pct": 0.3,
#         "module_type": 0,
#         "array_type": 1,
#         "om_cost_per_kw": 16.0,
#         "macrs_itc_reduction": 0.5,
#         "azimuth": 180.0,
#         "federal_rebate_per_kw": 350.0,
#         "dc_ac_ratio": 1.1
#     },
#     "LoadProfile": {
#         "doe_reference_name": "RetailStore",
#         "annual_kwh": 10000000.0,
#         "city": "LosAngeles"
#     },
#     "Storage": {
#         "total_rebate_per_kw": 100.0,
#         "macrs_option_years": 5,
#         "can_grid_charge": true,
#         "macrs_bonus_pct": 0.4,
#         "macrs_itc_reduction": 0.5,
#         "total_itc_pct": 0,
#         "installed_cost_per_kw": 1000.0,
#         "installed_cost_per_kwh": 500.0,
#         "replace_cost_per_kw": 460.0,
#         "replace_cost_per_kwh": 230.0
#     },
#     "ElectricTariff": {
#         "urdb_label": "5ed6c1a15457a3367add15ae"
#     },
#     "Financial": {
#         "escalation_pct": 0.026,
#         "offtaker_discount_pct": 0.081,
#         "owner_discount_pct": 0.081,
#         "analysis_years": 20,
#         "offtaker_tax_pct": 0.4,
#         "owner_tax_pct": 0.4,
#         "om_cost_escalation_pct": 0.025
#     }
# }}}


text:
using Publish
using DataAugmentation
using Pkg.Artifacts

Publish.Themes.default() = artifact"flux-theme"

p = Publish.Project(DataAugmentation)

rm("dev", recursive = true, force = true)
rm(p.env["version"], recursive = true, force = true)

deploy(DataAugmentation; root = "/DataAugmentation.jl", force = true, label = "dev")


text:
#---------------------------------------------------------------------#
#This routine advances the solution in time using Split-Explicit using
#a simple General Order RK method.
#Note that for M=3 in 1st stage, M=2 in 2nd, and M=1 in 3rd we need Qi in
#RHS function instead of Qk to get 3rd order convergence.
#
#As is, we get 1st order convergence but this does not work in NUMA
#
#Written by F.X. Giraldo / P.R. Mugg on 9/21/19
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#---------------------------------------------------------------------#
function split_explicit_WS_rk3_original(q0,c,M,Δt,time)

    #1st RK Stage
    Qi = q0
    t = time
    (rhs,f𝑖,g𝑖) = rhs_function(Qi,c,t)
    #M=3
    Δτ = Δt/M
    Qk = q0
    for k = 1:M/3
        t += Δτ*(k-1)
        #(rhs,f𝑘,g𝑘) = rhs_function(Qi,c,t)    #Qk->Qi
        (rhs,f𝑘,g𝑘) = rhs_function(Qk,c,t)
        Qk += Δτ*(f𝑘 + g𝑖)
    end
    Qi = Qk

    #2nd RK Stage
    #M=2
    t = time + Δt/3
    (rhs,f𝑖,g𝑖) = rhs_function(Qi,c,t)
    Δτ = Δt/M
    Qk = q0
    for k = 1:M/2
        t += Δτ*(k-1)
        #(rhs,f𝑘,g𝑘) = rhs_function(Qi,c,t)    #Qk->Qi
        (rhs,f𝑘,g𝑘) = rhs_function(Qk,c,t)
        Qk += Δτ*(f𝑘 + g𝑖)
    end
    Qi = Qk

    #3rd RK Stage
    #M=1; #Emil
    t = time + Δt/2
    (rhs,f𝑖,g𝑖) = rhs_function(Qi,c,t)
    Δτ = Δt/M
    Qk = q0
    for k = 1:M
        t += Δτ*(k-1)
        #(rhs,f𝑘,g𝑘) = rhs_function(Qi,c,t)    #Qk->Qi
        (rhs,f𝑘,g𝑘) = rhs_function(Qk,c,t)
        Qk += Δτ*(f𝑘 + g𝑖)
    end

    return qp = Qk
end #end function


text:
@testset "implementation.jl" begin
    # https://en.wikipedia.org/wiki/Cronbach's_alpha#Formula_and_calculation
    C₁ = [10  6  6  6;
           6 11  6  6;
           6  6 12  6;
           6  6  6 13]
    C₂ = [10  4  5  7;
           4 11  6  8;
           5  6 12  9;
           7  8  9 13]

    @test cronbach(C₁) ≈ 0.8135 rtol=8e-5
    @test cronbach(C₂) ≈ 0.8387 rtol=2e-5

    # http://www.abepro.org.br/biblioteca/enegep2010_TN_STO_131_840_16412.pdf
    A₃ = [3 2 3 1 1 3 4 4 3 4 4 4 3 4 4 4 3 1;
          2 4 3 4 5 5 5 5 1 3 4 3 4 4 4 3 3 3;
          4 2 3 1 1 3 4 4 4 4 4 4 3 4 4 4 3 1;
          1 1 3 3 1 1 3 1 1 1 1 3 3 3 1 1 3 1;
          2 4 4 4 3 4 4 4 5 4 3 2 4 4 4 3 2 2;
          1 3 3 4 4 4 3 1 1 1 4 4 3 2 4 3 1 3;
          3 4 3 1 3 4 3 3 3 4 3 4 3 4 4 3 1 3]
    C₃ = cov(A₃)

    @test cronbach(C₃) ≈ 0.84 rtol=5e-3
end


text:
for N in [Float64, Float32, Rational{Int}]

    c = zeros(N, 2)
    E1 = Matrix(Diagonal(N[-1, 0.5]))
    E2 = N[1 1; 0.5 0.3]
    E = [E1, E2]
    F2 = hcat(N[-0.5, 1])
    F = [F2]
    G = Matrix(N(0.3)I, 2, 2)
    p = PolynomialZonotope(c, E, F, G)

    @test dim(p) == 2
    @test order(p) == 7//2
    @test polynomial_order(p) == 2

    # type-specific concrete methods
    scale(N(0.5), p)
    linear_map(N[1.0 2.0; 2.0 5.0], p)
    z = Zonotope(N[1., 2.], Matrix(N(1)I, 2, 2))
    minkowski_sum(p, z)
    minkowski_sum(z, p)
end


text:
@testset "build fν max" begin
    @test all(LadderDGA.build_fνmax_fast([4,3,2,1,2,3,4],2) .== [11, 19])
    @test all(LadderDGA.build_fνmax_fast([3,2,1,2,3],3) .== [1, 5, 11])
    @test all(LadderDGA.build_fνmax_fast([3,2,1,2,3],2) .== [5, 11])
    @test all(LadderDGA.build_fνmax_fast([3,2,1,2,3],1) .== [11])
    @test all(LadderDGA.build_fνmax_fast([3,2,1,1,2,3],3) .== [2, 6, 12])
    @test all(LadderDGA.build_fνmax_fast([3,2,1,1,2,3],2) .== [6, 12])
    @test all(LadderDGA.build_fνmax_fast([3,2,1,1,2,3],1) .== [12])
    a = [3 3 3 3 3; 3 2 2 2 3; 3 2 1 2 3; 3 2 2 2 3; 3 3 3 3 3]
    test = [1, 1+8*2, 1+8*2+16*3]
    @test all(LadderDGA.build_fνmax_fast(a, 3) .== test)
    @test all(LadderDGA.build_fνmax_fast(a, 2) .== test[2:3])
    @test all(LadderDGA.build_fνmax_fast(a, 1) .== test[3])
end

@testset "get_sum_helper" begin
#    LadderDGA.default_fit_range()
    LadderDGA.get_sum_helper(1:4, sP_1, :f) === typeof(Naive())
    LadderDGA.get_sum_helper(1:12, sP_2, :f) === typeof(Richardson(1:3,0:2))
end

@testset "sum_freq" begin
    test = [1/n^k for k in 2:2:8, n in 1:10]
    @test all(LadderDGA.sum_freq(test, [2], Naive(), Float64(π)) .≈ sum(test,dims=[2]) ./ π)
end

@testset "extend_γ" begin
    #data for beta=12,u=1 => h = 2π/β=0.5235987755982988
    test_γ = [875.0225752325874 - 0.047426272065926985im, -1725.9367975032499 + 0.13273183259918772im, 1106.7342694260153 - 0.09340186405920287im, -226.68570367819405 + 0.014352094227780764im, 0.9896265166414565 + 0.00043227905695820295im, 0.9978757876736383 - 0.0006585999770301414im, 0.9986930184922074 - 0.0010106467457844685im, 0.9986261046489859 - 0.0012958932989215043im, 0.9983993208183312 - 0.0016338629583344207im, 0.9980949142745263 - 0.002078738753328137im, 0.9976988866855971 - 0.002685458038367005im, 0.9972007583275463 - 0.0035383199841694316im, 0.9965368770069696 - 0.004768152090402092im, 0.995645502135919 - 0.006591805844899053im, 0.9943660393968459 - 0.00937175752297385im, 0.9923873308819153 - 0.013731243746690323im, 0.9889066965298301 - 0.020694180905500778im, 0.9816721358991439 - 0.03167341201329606im, 0.9640108858763057 - 0.04655766186896375im, 0.9244117285078746 - 0.04832854114431357im, 0.9244117318979092 + 0.04832854220003501im, 0.9640108843948348 + 0.046557660890996555im, 0.9816721358830809 + 0.03167341149714696im, 0.98890669634837 + 0.020694180217158756im, 0.992387331274425 + 0.013731243737292193im, 0.9943660396961701 + 0.009371757480616041im, 0.9956455024259752 + 0.006591805870092188im, 0.9965368769515343 + 0.004768151457081569im, 0.9972007584368632 + 0.0035383197349404417im, 0.9976988864192321 + 0.0026854567958323953im, 0.998094914275765 + 0.002078738237549894im, 0.9983993208192198 + 0.0016338624421033839im, 0.9986261046497896 + 0.0012958927826615702im, 0.9986930184916067 + 0.0010106462296497864im, 0.9978757876733757 + 0.0006585994612941748im, 0.989626516641224 - 0.00043227956742801im, -226.68570367820325 - 0.014351947062101895im, 1106.734269426076 + 0.09340114633426014im, -1725.9367975033356 - 0.13273071298071334im, 875.022575232618 + 0.04742570463505424im]
    test_γ_2 = copy(test_γ)
    test_γ_3 = copy(test_γ)
    ref = ones(size(test_γ)) .+ 0im

    @test all(LadderDGA.find_usable_γ(test_γ; threshold=50, prct_red=0.05) .== (7,34))
    LadderDGA.extend_γ!(test_γ, 0.5235987755982988)
    LadderDGA.extend_γ!(test_γ_3, ref)
    @test all(imag.(test_γ[[1 2 3 4 36 37 38 39]]) .< 0.01)
    @test all(abs.(real.(test_γ[[1 2 3 4 36 37 38 39]]) .- 1.00) .< 0.01)
    @test all(test_γ[7:34] .== test_γ_2[7:34])
    @test all(abs.(test_γ_3[[1 2 3 4 36 37 38 39]]) .== 1)
    @test all(test_γ_3[7:34] .== test_γ_2[7:34])
end


text:
module Day24
__revise_mode__ = :eval

using Advent: save_input

const OFFSETS_1 = Dict(
    'e' => CartesianIndex(+1, -1, +0),
    'w' => CartesianIndex(-1, +1, +0),
)

const OFFSETS_2 = Dict(
    "se" => CartesianIndex(+0, -1, +1),
    "sw" => CartesianIndex(-1, +0, +1),
    "nw" => CartesianIndex(+0, +1, -1),
    "ne" => CartesianIndex(+1, +0, -1),
)

const OFFSETS = union(values(OFFSETS_1), values(OFFSETS_2))
const OFFSETS_self = union(OFFSETS, [CartesianIndex(+0, +0, +0)])

function getcoord(path)
    i = 1
    coord = CartesianIndex(0, 0, 0)
    @inbounds while i <= length(path)
        if path[i] == 'e' || path[i] == 'w'
            coord += OFFSETS_1[path[i]]
            i += 1
        else
            coord += OFFSETS_2[path[i:i + 1]]
            i += 2
        end
    end
    coord
end

function parse_input(filename)
    black = Set{CartesianIndex{3}}()
    for line in eachline(filename)
        c = getcoord(line)
        c in black ? pop!(black, c) : push!(black, c)
    end
    black
end

function solve_a(filename)
    filename |> parse_input |> length
end

countoccupied(S, I) = @inbounds count(O -> I + O in S, OFFSETS)

function simulate(S)
    next = deepcopy(S)
    low, high = extrema(S)
    for I in S, O in OFFSETS_self
        n = I + O
        c = countoccupied(S, n)
        if n in next
            (c == 0 || c > 2) && pop!(next, n)
        else
            c == 2 && push!(next, n)
        end
    end
    next
end

function solve_b(filename)
    black = parse_input(filename)
    for _ in 1:100
        black = simulate(black)
    end
    length(black)
end

solve(filename = save_input(2020, 24)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
function read_gms_file(filename::AbstractString)

    println("Reading $(filename) ...")

    filepath = joinpath(Pkg.dir("toJuMP"),".gms","")
    filepath = string(filepath,filename,".gms")

    if isfile(filepath)
        f = open(filepath, "r")
    elseif isfile(filename)
        f = open(filename, "r")
    else
        error("No gms file detected.")
    end

    #=
        ----------------------- Line Header Space ------------------
            Example :: 4stufen.gms
        Unkown ::
        $offlisting
        $if set nostart $goto modeldef
        $label modeldef
        $if not set MINLP $set MINLP MINLP  Model Type

        *                 -> comment
        Null              -> emptyline
        Variables         -> All Variable List, separated by ","
        Positive Variable -> Non-negative variables list, separated by "," (assert)
        Binary Variable   -> Binary variables list, separated by "," (assert)
        Equations         -> Equations symbols, separated by ","
        EquSymbols..      -> Equation content
        VarSymbols.lo     -> Lower Bound
        VarSymbols.l      -> Upper Bound ??
        Model m / all /   -> ??
        m.limrow=0
        m.
        ------------------------------------------------------------
    =#

    # Null Trackers
    skip = 0
    comment = 0

    # Main data structure initialization
    gms = oneProblem()

    # -------------------------- Start Reading -------------------------- #
    for l in eachline(f)
        sl = lstrip(l, ' ')
        if isempty(sl) || length(sl) == 1  # Skip empty line
            skip += 1
        elseif sl[1] == "*"
            comment += 1
        else
            slt = split(sl, r" |,|\.")
            slt = [slt[i] for i in 1:length(slt) if length(slt[i]) > 0]
            if slt[1] in GMS_BLOCK_HEADER
                read_block(f, gms, l)           # Reading a general block
            elseif slt[1] in gms.rows #ismatch(r"e\d", sl[1])
                read_equation(f, gms, l)        # Reading an equation block
            elseif slt[1] in gms.cols #ismatch(r"\w\.\w", sl[1])
                read_bounds(f, gms, l)          # Reading a bound block
            elseif slt[1] == "\$"
                read_command(f, gms, l)         # Reading a logical command line
            elseif slt[1] == "Model"
                read_model(f, gms, sl)           # Reading a model command
            elseif slt[1] == "Solve"
                read_solve(f, gms, sl)           # Reading a solve command line
            else
                # println("[Unkown]DUMPING >> ", l)
            end
        end

    end
    close(f)
    return gms
end

function read_block(file::IOStream, gms::oneProblem, lInit::AbstractString; kwargs...)

    sl = split(get_one_line(file, one_line=lInit), r" |,")
    sl = [sl[i] for i in 1:length(sl) if !isempty(sl[i])]  # Eliminate empty entries

    blockIdentifier = sl[1]  # This should always be true

    for i in 1:length(sl)
        if blockIdentifier == "Variables" && (i > 1)
            push!(gms.cols, sl[i])
        elseif blockIdentifier == "Positive" && (i > 2)
            gms.colsType[sl[i]] = "Positive"
        elseif blockIdentifier == "Negative" && (i > 2)
            gms.colsType[sl[i]] = "Negative"
        elseif blockIdentifier == "Binary" && (i > 2)
            gms.colsType[sl[i]] = "Binary"
        elseif blockIdentifier == "Semicont"
            error("Currently don't support semi continous variables parsing.")
        elseif blockIdentifier == "Integer" && (i > 2)
            gms.colsType[sl[i]] = "Integer"
        elseif blockIdentifier == "Equations" && (i > 1)
            push!(gms.rows, sl[i])
        end
    end

    return
end

function read_equation(file::IOStream, gms::oneProblem, lInit::AbstractString; kwargs...)

    #=
        Example Equation
        #Single Line
        e1..    objvar - x145 - x146 - x147 - x148 - x149 - x150 =E= 3271.22725820856;
        # Multi Line
        e38..    109.0495*x1 + 100.462*x2 + 115.2937*x3 + 2.860271*x4 + 15.1404*x5
               + x44 - x556 =E= 101.1304;
    =#

    lhs = ""
    rhs = ""
    sense = ""
	eN = string(strip(split(lInit, r" |,")[1], ['.', '\r', '\n']))
    eS = split(eN, '.')
    if eN in gms.rows
        one_l = get_one_line(file, one_line=lInit)
        one_l = replace(one_l, "$(eN)..", "")
        if contains(one_l, "=E=")
            gms.rowsSense[eN] = "E"
        elseif contains(one_l, "=L=")
            gms.rowsSense[eN] = "L"
        elseif contains(one_l, "=G=")
            gms.rowsSense[eN] = "G"
		elseif contains(one_l, "=N=")
			gms.rowsSense[eN] = "N"
        else
            error("NO sense detected in equation. $(one_l)")
        end
        sl = split(one_l, "=$(gms.rowsSense[eN])=", keep=false)
        gms.rowsLHS[eN] = replace(sl[1], " ", "")
        gms.rowsRHS[eN] = Float64(parse(sl[2]))
    elseif length(eS) == 2 && eS[2] == "m"
        @assert eS[1] in gms.rows
        one_l = get_one_line(file, one_line=lInit)
        mVal = parse(split(one_l, "=", keep=false)[2])
        gms.m[eS[1]] = mVal
    else
        error("Non-indexed constraint detected. \n$(lInit)")
    end

    return
end

function read_bounds(file::IOStream, gms::oneProblem, lInit::AbstractString; kwargs...)

    #=
        .lo lower bound
        .l level or primal value
        .up upper bound
        .m marginal or dual value
    =#
	@assert rstrip(strip(lInit, ['\n','\r']), ' ')[end] == ';' # Assumption :: alwasy one line
    all_segs = split(lInit, ';', keep=false)
    all_segs = [i for i in all_segs if !isempty(strip(i,' '))]
    for lseg in all_segs
        lseg = strip(lseg, [';', '\n', ' ', '\r'])
        sl = split(lseg, "=")
        sl = [sl[i] for i in 1:length(sl) if !isempty(sl[i])]
        # There shouldn't be any space to begin with
        sl[1] = strip(sl[1], [';', '\n', ' ', '\r'])
		sl[2] = strip(sl[2], [';', '\n', ' ', '\r'])
        @assert ismatch(r"\w.\w", sl[1])
        boundVal = Float64(parse(sl[end]))
        var = split(sl[1],".")[1]
        attr = split(sl[1],".")[2]

        if attr == "lo"
            gms.lb[var] = boundVal
        elseif attr == "l"
            gms.l[var] = boundVal
        elseif attr == "fx"
            gms.fx[var] = boundVal
        elseif attr == "up"
            gms.ub[var] = boundVal
        elseif attr == "m"
            gms.m[var] = boundVal
        elseif attr == "scale"
            gms.scale[var] = boundVal
        elseif attr == "prior"
            gms.prior[var] = boundVal
        elseif attr == "stage"
            gms.prior[var] = boundVal
        else
            error("Unknown variable attribute. Consult gams hand-book for definition and report this issue.")
        end
    end

    return
end

function read_model(file::IOStream, gms::oneProblem, lInit::AbstractString; kwargs...)

    # Assume this will always be one line
	@assert strip(lInit, ['\n','\r'])[end] == ';'
    sl = split(strip(lInit,[';','\n','\r']), r" |,")
    @assert sl[1] == "Model"  # Already did outside
    # Not sure if the printlnrmation stored here is useful or not.
    for i in 1:length(sl)
        if sl[i] == "Model"
            gms.modelSymbol = sl[i+1]
        end
    end

    return
end

function read_solve(file::IOStream, gms::oneProblem, lInit::AbstractString; kwargs...)

	@assert strip(lInit, ['\n','\r'])[end] == ';'

    sl = split(strip(lInit, [';','\n','\r']), r" |,")
    @assert sl[1] == "Solve"

    for i in 1:length(sl)
        if sl[i] in PROBTYPE
            gms.modelType = sl[i]
        elseif sl[i] == "maximizing"
            gms.objSense = sl[i]
        elseif sl[i] == "minimizing"
            gms.objSense = sl[i]
        elseif sl[i] in gms.cols
            gms.objective = sl[i]
        end
    end
end

function gms2jump(gmsPath::AbstractString;
                    mode::AbstractString="index",
                    ending="m=m",
                    quadNL::Bool=false,
                    loopifpossible::Bool=true,
                    outdir::AbstractString="")

    probName = replace(splitdir(gmsPath)[end],".gms", "")
    gms = read_gms_file(gmsPath)
    write_julia_script(probName, gms, mode, loopifpossible=loopifpossible, quadNL=quadNL, outdir=outdir)

    return
end


text:
function div_thread(N :: Int64, nt :: Int64)

	npt = zeros(Int64, nt, 2)

	n = div(N, nt)

	for i in 1:nt
		npt[i, 1] = (i-1)*n + 1
		npt[i, 2] = npt[i, 1]-1+n
	end

	npt[nt, 2] += mod(N, nt)

	return npt
end


text:
#!/usr/bin/env julia --project
using CLIMA
CLIMA.init()
using CLIMA.Atmos
using CLIMA.ConfigTypes
using CLIMA.Diagnostics
using CLIMA.GenericCallbacks
using CLIMA.ODESolvers
using CLIMA.ColumnwiseLUSolver: ManyColumnLU
using CLIMA.Mesh.Filters
using CLIMA.Mesh.Grids
using CLIMA.MoistThermodynamics: air_temperature, internal_energy, air_pressure
using CLIMA.VariableTemplates

using Distributions: Uniform
using LinearAlgebra
using StaticArrays
using Random: rand
using Test

using CLIMAParameters
using CLIMAParameters.Planet: R_d, day, grav, cp_d, cv_d, planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

struct HeldSuarezDataConfig{FT}
    T_ref::FT
end

function init_heldsuarez!(bl, state, aux, coords, t)
    FT = eltype(state)

    # Set initial state to reference state with random perturbation
    rnd = FT(1.0 + rand(Uniform(-1e-3, 1e-3)))
    state.ρ = aux.ref_state.ρ
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.ρe = rnd * aux.ref_state.ρe

    nothing
end

function config_heldsuarez(FT, poly_order, resolution)
    # Set up a reference state for linearization of equations
    T_sfc::FT = 290 # surface temperature of reference state (K)
    ΔT::FT = 60     # temperature drop between surface and top of atmosphere (K)
    H_t::FT = 8e3   # height sclae over which temperature drops (m)
    temp_profile_ref = DecayingTemperatureProfile(T_sfc, ΔT, H_t)
    ref_state = HydrostaticState(temp_profile_ref, FT(0))

    # Set up a Rayleigh sponge to dampen flow at the top of the domain
    domain_height::FT = 30e3               # distance between surface and top of atmosphere (m)
    z_sponge::FT = 12e3                    # height at which sponge begins (m)
    α_relax::FT = 1 / 60 / 15              # sponge relaxation rate (1/s)
    exp_sponge = 2                         # sponge exponent for squared-sinusoid profile
    u_relax = SVector(FT(0), FT(0), FT(0)) # relaxation velocity (m/s)
    sponge = RayleighSponge{FT}(
        domain_height,
        z_sponge,
        α_relax,
        u_relax,
        exp_sponge,
    )

    # Set up the atmosphere model
    exp_name = "HeldSuarez"
    T_ref::FT = 255        # reference temperature for Held-Suarez forcing (K)
    τ_hyper::FT = 4 * 3600 # hyperdiffusion time scale in (s)
    c_smag::FT = 0.21      # Smagorinsky coefficient
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        param_set;
        ref_state = ref_state,
        turbulence = SmagorinskyLilly(c_smag),
        hyperdiffusion = StandardHyperDiffusion(τ_hyper),
        moisture = DryModel(),
        source = (Gravity(), Coriolis(), held_suarez_forcing!, sponge),
        init_state_conservative = init_heldsuarez!,
        data_config = HeldSuarezDataConfig(T_ref),
    )

    config = CLIMA.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_heldsuarez!;
        model = model,
    )

    return config
end

function held_suarez_forcing!(
    bl,
    source,
    state,
    diffusive,
    aux,
    t::Real,
    direction,
)
    FT = eltype(state)

    # Parameters
    T_ref = bl.data_config.T_ref

    # Extract the state
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe

    coord = aux.coord
    e_int = internal_energy(bl.moisture, bl.orientation, state, aux)
    T = air_temperature(bl.param_set, e_int)
    _R_d = FT(R_d(bl.param_set))
    _day = FT(day(bl.param_set))
    _grav = FT(grav(bl.param_set))
    _cp_d = FT(cp_d(bl.param_set))
    _cv_d = FT(cv_d(bl.param_set))

    # Held-Suarez parameters
    k_a = FT(1 / (40 * _day))
    k_f = FT(1 / _day)
    k_s = FT(1 / (4 * _day))
    ΔT_y = FT(60)
    Δθ_z = FT(10)
    T_equator = FT(315)
    T_min = FT(200)
    σ_b = FT(7 / 10)

    # Held-Suarez forcing
    φ = latitude(bl.orientation, aux)
    p = air_pressure(bl.param_set, T, ρ)

    #TODO: replace _p0 with dynamic surfce pressure in Δσ calculations to account
    #for topography, but leave unchanged for calculations of σ involved in T_equil
    _p0 = 1.01325e5
    σ = p / _p0
    exner_p = σ^(_R_d / _cp_d)
    Δσ = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil = (T_equator - ΔT_y * sin(φ)^2 - Δθ_z * log(σ) * cos(φ)^2) * exner_p
    T_equil = max(T_min, T_equil)
    k_T = k_a + (k_s - k_a) * height_factor * cos(φ)^4
    k_v = k_f * height_factor

    # Apply Held-Suarez forcing
    source.ρu -= k_v * projection_tangential(bl, aux, ρu)
    source.ρe -= k_T * ρ * _cv_d * (T - T_equil)
    return nothing
end

function config_diagnostics(FT, driver_config)
    interval = "100000steps" # chosen to allow a single diagnostics collection

    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(10), FT(10), FT(1000)) # in (deg, deg, m)
    interpol =
        CLIMA.InterpolationConfiguration(driver_config, boundaries, resolution)

    dgngrp = setup_dump_state_and_aux_diagnostics(
        interval,
        driver_config.name,
        interpol = interpol,
        project = true,
    )
    return CLIMA.DiagnosticsConfiguration([dgngrp])
end

function main()
    # Driver configuration parameters
    FT = Float32                             # floating type precision
    poly_order = 5                           # discontinuous Galerkin polynomial order
    n_horz = 5                               # horizontal element number
    n_vert = 5                               # vertical element number
    n_days = 120                             # experiment day number
    timestart = FT(0)                        # start time (s)
    timeend = FT(n_days * day(param_set))    # end time (s)

    # Set up driver configuration
    driver_config = config_heldsuarez(FT, poly_order, (n_horz, n_vert))

    # Set up experiment
    solver_config = CLIMA.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = 0.2,
        init_on_cpu = true,
        CFL_direction = HorizontalDirection(),
        diffdir = HorizontalDirection(),
    )

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    # Set up user-defined callbacks
    filterorder = 10
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            1:size(solver_config.Q, 2),
            solver_config.dg.grid,
            filter,
        )
        nothing
    end

    # Run the model
    result = CLIMA.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        check_euclidean_distance = true,
    )
end

main()


text:
"""
    ozofβ(β [, T = Float64])
return propagator exp(β σz)
"""
function ozofβ(β, T = Float64)
    σz = T[1 0; 0 -1]
    DTensor(exp(β*σz))
end

"""
    oxxofβ(β [, T = ComplexF64])
calculates the propagator for ∑_i σx_i σx_i+1, i.e.
    exp(β (∑_i σx_i σx_i+1 ))
according to Pirvu et al. 2010 _Matrix Product Operator Representations_
"""
function oxxofβ(β, T = Float64)
    id = T[1 0; 0 1]
    σx = T[0 1; 1 0]
    C0 = T[cosh(β) 0; 0 sinh(β)]
    C1 = sqrt(T(sinh(β)*cosh(β))) .* σx
    @tensor mpoxx[p1,v1,p2,v2] := C0[v1,v2] * id[p1,p2] + C1[v1,v2] * σx[p1,p2]
    return DTensor(mpoxx)
end

"""
    tfisingpropagator(β [, h = 1, T = ComplexF64])
return the MPO representing a second order trotter-suzuki expansion, i.e.
    U ≈ exp(β/2 h(σz⊗𝟙 + 𝟙⊗σz)) * exp(β σx ⊗ σx) * exp(β/2 h(σz⊗𝟙 + 𝟙⊗σz))
"""
function tfisingpropagator(β, h = 1, T = Float64)
    oz  = ozofβ(h*β/2,T)
    oxx = oxxofβ(β,T)
    @tensor op[1,2,3,4] := oz[1,-1] * oxx[-1,2,-3,4] * oz[-3,3]
    return op
end

"""
    tfisinghamiltonian([twobody = false, T = ComplexF64])
return the MPO representing the Hamiltonian
    H = -∑_i σx_i σx_i+1 - ∑_i σz_i
if `twobody=true`, return the two-body Hamiltonian
    H = -σx ⊗ σx - 1 ⊗ σz - σz ⊗ 1
"""
function tfisinghamiltonian(twobody = false, λ=1, T = Float64)
    id = T[1 0; 0 1]
    σx = T[0 1; 1 0]
    σz = T[1 0; 0 -1]

    h = DTensor(zeros(2,3,2,3))
    h[:,1,:,1] = h[:,3,:,3] =  id
    h[:,1,:,2] =  σx
    h[:,2,:,3] = -σx
    h[:,1,:,3] = -λ*σz
    twobody || return h

    rb = DTensor([0,0,1])
    lb = DTensor([1,0,0])
    @tensor op[1,2,3,4] := rb[-1] * h[1,-3,3,-1] * h[2,-2,4,-3] * lb[-2]
    return op
end

function mag(st)
    C, T = st.C, st.Ts[1]
    σx = DTensor([0 1; 1  0])
    @tensor begin
        env[1,2] := C[-1,-2] * T[-2,2,-3] * conj(C)[-3,-4] *
                    C[-4,-5] * conj(T)[-5,1,-6] * conj(C)[-6,-1]
        m = env[-1,-2] * σx[-1,-2]
        n = env[-1,-1]
    end
    return m/n
end

function energy(st)
    C, T = st.C, st.Ts[1]
    h = tfisinghamiltonian(true)
    @tensor begin
        env[1,2,3,4] := C[-1,-2] * T[-2,3,-3] * T[-3,4,-4] * conj(C)[-4,-5] *
                        C[-5,-6] * conj(T)[-6,2,-7] * conj(T)[-7,1,-8] *
                        conj(C)[-8,-1]
        ee = env[-1,-2,-3,-4] * h[-1,-2,-3,-4]
        n  = env[-1,-2,-1,-2]
    end
    return (ee/n)/2
end

function tfisingctm(β,χ, h = 1;kwargs...)
    iter = transconjctmiterable(tfisingpropagator(β,h),χ)
    return ctm_kernel(iter; kwargs...)
end


text:
if length(ARGS) != 0 && ARGS[1] == "local"
    include("local.jl"); exit()
end

names = [
    "test_operators_dense.jl",
    "test_superoperators.jl",
    "test_metrics.jl",
    "test_molecules.jl",
    "test_aggregate.jl",
    "test_evolution.jl",
    "test_schroedinger.jl",
    "test_liouville.jl",
    "test_interaction_picture.jl",
    "test_master.jl",
    "test_trace.jl",
    "test_initial_state.jl",
    "test_memory_kernel.jl",
    "test_master_ansatz.jl",
]

detected_tests =
    filter(name -> startswith(name, "test_") && endswith(name, ".jl"), readdir("."))

unused_tests = setdiff(detected_tests, names)
if length(unused_tests) != 0
    error("The following tests are not used:\n", join(unused_tests, "\n"))
end

unavailable_tests = setdiff(names, detected_tests)
if length(unavailable_tests) != 0
    error("The following tests could not be found:\n", join(unavailable_tests, "\n"))
end

if length(ARGS) != 0
    names = map(x -> "test_" * x * ".jl", ARGS)
end
# names = ["test_trace.jl"]
for name in names
    if startswith(name, "test_") && endswith(name, ".jl")
        include(name)
    end
end


text:
function compute_zenith_delay(obj::DiscreteTroposphericModel, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::FieldAbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{RealFieldElement}, (RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2)
end

function compute_zenith_delay(obj::DiscreteTroposphericModel, arg0::jdouble, arg1::Vector{jdouble}, arg2::AbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{jdouble}, (jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2)
end

function get_parameters(obj::MappingFunction)
    return jcall(obj, "getParameters", Vector{jdouble}, ())
end

function get_parameters(obj::MappingFunction, arg0::Field)
    return jcall(obj, "getParameters", Vector{RealFieldElement}, (Field,), arg0)
end

function get_parameters_drivers(obj::MappingFunction)
    return jcall(obj, "getParametersDrivers", List, ())
end

function mapping_factors(obj::MappingFunction, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{RealFieldElement}, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function mapping_factors(obj::MappingFunction, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{jdouble}, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function path_delay(obj::DiscreteTroposphericModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "pathDelay", RealFieldElement, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function path_delay(obj::DiscreteTroposphericModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "pathDelay", jdouble, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end


text:
module TestBasic
using Test
using RData

@testset "Basic/Utility RData functions" begin

# check for Float64 NA
@testset "Detect R floating-point NAs" begin
    @test !RData.isna_float64(reinterpret(UInt64, 1.0))
    @test !RData.isna(1.0)
    @test !RData.isna(NaN)
    @test !RData.isna(Inf)
    @test !RData.isna(-Inf)
    @test RData.isna_float64(RData.R_NA_FLOAT64)
    # check that alternative NA is also recognized (#10)
    @test RData.isna_float64(reinterpret(UInt64, RData.R_NA_FLOAT64 | ((Base.significand_mask(Float64) + 1) >> 1)))
end

@testset "Unescape R strings" begin
    @test RData.unescape_rstring("\\\"") == "\""
    @test RData.unescape_rstring("\\'") == "'"
end

@testset "Internals" begin
    @testset "sxtype()" begin
        v = RData.RVector{Int, RData.INTSXP}([1, 2], RData.emptyhash)
        @test RData.sxtype(typeof(v)) == RData.INTSXP
        @test RData.sxtype(v) == RData.INTSXP
        s = RData.RSymbol("abc")
        @test RData.sxtype(typeof(s)) == RData.SYMSXP
        @test RData.sxtype(s) == RData.SYMSXP
    end

    @testset "getdata()" begin
        v = RData.RVector{Int, RData.INTSXP}([1, 2, 3], RData.emptyhash)
        @test RData.getdata(v) == [1, 2, 3]
        @test RData.getdata(v) === v.data
        @test_throws UnsupportedROBJ RData.getdata(RData.RSymbol("abc"))
    end

    @testset "addattr()" begin
        v = RData.RVector{Int, RData.INTSXP}([1, 2, 3], RData.emptyhash)
        v2 = RData.addattr(v)
        @test v2 isa RData.RVector{Int, RData.INTSXP}
        @test v2.data === v.data
        @test v2.attr !== RData.emptyhash
    end

    @testset "getattr()" begin
        v = RData.addattr(RData.RVector{Int, RData.INTSXP}([1, 2, 3], RData.emptyhash))
        @test_throws KeyError RData.getattr(v, "names")
        @test RData.getattr(v, "names", ["a", "b", "c"]) == ["a", "b", "c"]
    end
end

end

end # TestBasic


text:
using Documenter
using DocumenterCitations

bib = CitationBibliography(joinpath(@__DIR__, "example.bib"))

makedocs(
    bib,
    sitename = "DocumenterCitations.jl",
    strict = true,
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Home"       => "index.md",
        "References" => "references.md"
    ]
)

deploydocs(
    repo = "github.com/ali-ramadhan/DocumenterCitations.jl.git",
    push_preview = true
)


text:
export vlplot_axes

function asvlplot_scale(scale)
  if scale in [:linear,:log,:sqrt,:time,:utc,:sequential]
    string(scale)
  else
    error("Unsupported scale type $scale")
  end
end

function vlplot_axes(data,args...;colors="reds",width=300,height=300,kwds...)
  df, axes = asplotable(data,args...;kwds...)
  if any(x -> x isa QualitativePlotAxis,axes)
    error("VegaLite interface does not yet support qualitative axes.")
  end

  vlplot_axes_(df,axes,names(df)[2:end]...;colors=colors,
               width=width,height=height)
end

function vlplot_axes_(df,axes,x;colors,width,height)
  df |>
    @vlplot(:line, width=width, height=height,
            x={field=x,typ="quantitative", bin={step=axes[1].step},
               scale={typ=asvlplot_scale(axes[1].scale)}},
            y={field=:value,typ="quantitative"})
end

function vlplot_axes_(df,axes,x,y;colors,width,height)
  df |>
  @vlplot(:rect, width=width, height=height,
          x={field=x,typ="quantitative", bin={step=axes[1].step},
             scale={typ=asvlplot_scale(axes[1].scale)}},
          y={field=y,typ="quantitative", bin={step=axes[2].step},
             scale={typ=asvlplot_scale(axes[2].scale)}},
          color={field=:value, aggregate="mean", typ="quantitative"},
          config={view={stroke="transparent"},
                  scale={bandPaddingInner=0, bandPaddingOuter=0},
                  range={heatmap={scheme=colors}}})
end

function vlplot_axes_(df,axes,x,y,z;colors,width,height)
  df |>
    @vlplot(:rect, width=width, height=height,
            x={field=x,typ="quantitative", bin={step=axes[1].step},
               scale={typ=asvlplot_scale(axes[1].scale)}},
            y={field=y,typ="quantitative", bin={step=axes[2].step},
               scale={typ=asvlplot_scale(axes[2].scale)}},
            color={field=:value, aggregate="mean", typ="quantitative"},
            column={field=z},
            config={view={stroke="transparent"},
                    scale={bandPaddingInner=0, bandPaddingOuter=0},
                    range={heatmap={scheme=colors}}})
end

function vlplot_axes_(df,axes,x,y,z,w;colors,width,height)
  df |>
    @vlplot(:rect, width=width, height=height,
            x={field=x,typ="quantitative", bin={step=axes[1].step},
               scale={typ=asvlplot_scale(axes[1].scale)}},
            y={field=y,typ="quantitative", bin={step=axes[2].step},
               scale={typ=asvlplot_scale(axes[2].scale)}},
            color={field=:value, aggregate="mean", typ="quantitative"},
            column={field=z},row={field=w},
            config={view={stroke="transparent"},
                    scale={bandPaddingInner=0, bandPaddingOuter=0},
                    range={heatmap={scheme=colors}}})
end

function vlplot_axes_(df,axes,args...;kwds...)
  error("Plotting data with $(length(axes)) dims along
        $(length(args)) axes is not supported.")
end


text:
using GeometricBase.Config
using GeometricBase.Utils
using GeometricIntegrators.Integrators
using GeometricIntegrators.Tableaus
using GeometricProblems.LotkaVolterra2d
using SimpleSolvers
using Test

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

const Δt = 0.01
const nt = 10
const q₀ = [1.0, 1.0]
const parameters = (a₁=1.0, a₂=1.0, b₁=-1.0, b₂=-2.0)

ode  = lotka_volterra_2d_ode(q₀; params=parameters)
iode = lotka_volterra_2d_iode(q₀; params=parameters)
lode = lotka_volterra_2d_lode(q₀; params=parameters)
int  = IntegratorFIRK(ode, TableauGauss(8), Δt)
sol  = integrate(ode, int, nt)
refx = sol.q[end]

@testset "$(rpad("Runge-Kutta integrators for implicit equations",80))" begin

    int = IntegratorFIRKimplicit(iode, TableauGauss(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-5

    int = IntegratorFIRKimplicit(iode, TableauGauss(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-7

    int = IntegratorFIRKimplicit(iode, TableauGauss(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-10

    int = IntegratorFIRKimplicit(iode, TableauGauss(4), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-14

    int = IntegratorSRKimplicit(iode, TableauGauss(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-5

    int = IntegratorSRKimplicit(iode, TableauGauss(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-7

    int = IntegratorSRKimplicit(iode, TableauGauss(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-10

    int = IntegratorSRKimplicit(iode, TableauGauss(4), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-14

    int = IntegratorPRKimplicit(iode, PartitionedTableauGauss(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    int = IntegratorPRKimplicit(iode, PartitionedTableauGauss(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-7

    int = IntegratorPRKimplicit(iode, PartitionedTableauGauss(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-12

    # int = IntegratorPRKimplicit(iode, PartitionedTableauGauss(4), Δt)
    # sol = integrate(iode, int, nt)
    # @test relative_maximum_error(sol.q, refx) < 2E-16

end

@testset "$(rpad("Special integrators",80))" begin

    flint = IntegratorFLRK(lode, TableauGauss(2), Δt)
    flsol = integrate(lode, flint, nt)
    @test relative_maximum_error(flsol.q, refx) < 4E-12

    flint = IntegratorFLRK(lode, TableauGauss(3), Δt)
    flsol = integrate(lode, flint, nt)
    @test relative_maximum_error(flsol.q, refx) < 4E-16

    flint = IntegratorFLRK(lode, TableauGauss(4), Δt)
    flsol = integrate(lode, flint, nt)
    @test relative_maximum_error(flsol.q, refx) < 8E-16

    pgint = IntegratorPGLRK(ode, CoefficientsPGLRK(2), Δt)
    pgsol = integrate(ode, pgint, nt)
    @test relative_maximum_error(pgsol.q, refx) < 4E-12

    pgint = IntegratorPGLRK(ode, CoefficientsPGLRK(3), Δt)
    pgsol = integrate(ode, pgint, nt)
    @test relative_maximum_error(pgsol.q, refx) < 4E-16

    pgint = IntegratorPGLRK(ode, CoefficientsPGLRK(4), Δt)
    pgsol = integrate(ode, pgint, nt)
    @test relative_maximum_error(pgsol.q, refx) < 4E-16

end


text:
using CSV
using DataFrames
using DiscreteChoiceModels

data = CSV.read("../data/biogeme_swissmetro.dat", DataFrame, delim='\t')
data = data[in.(data.PURPOSE, [Set([1, 3])]) .& (data.CHOICE .!= 0), :]

@assert nrow(data) == 6768

model = multinomial_logit(
    @utility(begin
        1 ~ :αtrain + :βtravel_time * TRAIN_TT / 100 + :βcost * (TRAIN_CO * (GA == 0)) / 100
        2 ~ :αswissmetro + :βtravel_time * SM_TT / 100 + :βcost * SM_CO * (GA == 0) / 100
        3 ~ :αcar + :βtravel_time * CAR_TT / 100 + :βcost * CAR_CO / 100

        :αswissmetro = 0f  # fix swissmetro ASC to zero
    end),
    data.CHOICE,
    data,
    availability=[
        1 => (data.TRAIN_AV .== 1) .& (data.SP .!= 0),
        2 => data.SM_AV .== 1,
        3 => (data.CAR_AV .== 1) .& (data.SP .!= 0),
    ]
)

println(DiscreteChoiceModels.summary(model))


text:
import Base: keys, getindex, +
export +, getindex

"""

Represents a sweep of one or more parameters of a `HarmonicEquation`.
During a sweep, the selected parameters vary linearly over some timespan and are constant elsewhere.

Sweeps of different variables can be combined using `+`.

# Fields
$(TYPEDFIELDS)

## Examples
```julia-repl
# create a sweep of parameter a from 0 to 1 over time 0 -> 100
julia> @variables a,b;
julia> sweep = ParameterSweep(a => [0., 1.], (0, 100));
julia> sweep[a](50)
0.5
julia> sweep[a](200)
1.0

# do the same, varying two parameters simultaneously
julia> sweep = ParameterSweep([a => [0.,1.], b => [0., 1.]], (0,100))
```
"""
struct ParameterSweep
    """Maps each swept parameter to a function."""
    functions::Dict{Num, Function}

    ParameterSweep(functions...) = new(Dict(functions...))
    ParameterSweep() = ParameterSweep([])

end

# overload so that ParameterSweep can be accessed like a Dict
keys(s::ParameterSweep) = keys(s.functions)
getindex(s::ParameterSweep, i) = getindex(s.functions, i)

# overload +
function +(s1::ParameterSweep, s2::ParameterSweep)
    common_params = intersect(keys(s1), keys(s2))
    !isempty(common_params) && error("cannot combine sweeps of the same parameter")
    return ParameterSweep(merge(s1.functions, s2.functions))

    # combine sweeps of the same parameter
    #interval_overlap(s1.timespan, s2.timespan) && error("cannot combine sweeps with overlapping timespans")
    #new_funcs = filter(x -> !in(x.first, common_params), all_params)
end

#interval_overlap((t1, t2), (t3, t4)) = t4 > t2 > t3 || t2 > t4 > t1


text:
import Statistics
import Random
import Zygote
import RestrictedBoltzmannMachines as RBMs

using Test: @test, @testset
using Random: bitrand, randn!
using Statistics: mean
using LinearAlgebra: logdet, diag, diagm, dot, isposdef
using LogExpFunctions: logsumexp
using QuadGK: quadgk
using RestrictedBoltzmannMachines: BinaryRBM, visible, hidden, weights
using RestrictedBoltzmannMachines: inputs_h_from_v, inputs_v_from_h, inputs_v_to_h, inputs_h_to_v
using RestrictedBoltzmannMachines: energy, interaction_energy, free_energy

@testset "batches, n=$n, m=$m, Bv=$Bv, Bh=$Bh" for n in (5, (5,2)), m in (2, (3,4)), Bv in ((), (3,2)), Bh in ((), (3,2))
    rbm = RBMs.BinaryRBM(randn(n...), randn(m...), randn(n..., m...))
    wmat = reshape(weights(rbm), length(visible(rbm)), length(hidden(rbm)))
    v = bitrand(n..., Bv...)
    h = bitrand(m..., Bh...)

    @test RBMs.batch_size(visible(rbm), v) == Bv
    @test RBMs.batch_size(hidden(rbm), h) == Bh

    @test size(inputs_h_from_v(rbm, v)) == (size(hidden(rbm))...,  Bv...)
    @test size(inputs_v_from_h(rbm, h)) == (size(visible(rbm))..., Bh...)
    @test inputs_h_from_v(rbm, v) == inputs_v_to_h(rbm, v)
    @test inputs_v_from_h(rbm, h) == inputs_h_to_v(rbm, h)

    if length(Bv) == length(Bh) == 0
        @test interaction_energy(rbm, v, h) isa Number
        @test interaction_energy(rbm, v, h) ≈ -vec(v)' * wmat * vec(h)
        @test energy(rbm, v, h) isa Number
    elseif length(Bv) == 0
        hmat = reshape(h, length(hidden(rbm)), :)
        E = -vec(v)' * wmat * hmat
        @test interaction_energy(rbm, v, h) isa AbstractArray
        @test size(interaction_energy(rbm, v, h)) == Bh
        @test interaction_energy(rbm, v, h) ≈ reshape(E, Bh)
        @test size(energy(rbm, v, h)) == Bh
    elseif length(Bh) == 0
        vmat = reshape(v, length(visible(rbm)), :)
        E = -vec(h)' * wmat' * vmat
        @test interaction_energy(rbm, v, h) isa AbstractArray
        @test size(interaction_energy(rbm, v, h)) == Bv
        @test interaction_energy(rbm, v, h) ≈ reshape(E, Bv)
        @test size(energy(rbm, v, h)) == Bv
    else
        vmat = reshape(v, length(visible(rbm)), :)
        hmat = reshape(h, length(hidden(rbm)), :)
        E = -dot.(eachcol(vmat), Ref(wmat), eachcol(hmat))
        @test interaction_energy(rbm, v, h) isa AbstractArray
        @test size(interaction_energy(rbm, v, h)) == Bv == Bh
        @test interaction_energy(rbm, v, h) ≈ reshape(E, Bv)
        @test size(energy(rbm, v, h)) == Bv == Bh
    end
    @inferred inputs_h_from_v(rbm, v)
    @inferred inputs_v_from_h(rbm, h)
    @inferred interaction_energy(rbm, v, h)
    @inferred energy(rbm, v, h)
    gs = Zygote.gradient(rbm) do rbm
        mean(energy(rbm, v, h))
    end
    ∂w = @inferred RBMs.∂interaction_energy(rbm, v, h)
    @test ∂w ≈ only(gs).w
end

@testset "singleton batch dims" begin # behavior useful for ConvolutionalRBMs.jl
    rbm = RBMs.BinaryRBM(randn(3), randn(2), randn(3,2))

    v = bitrand(3,1,2)
    h = bitrand(2,3,1)
    @test size(energy(rbm, v, h)) == @inferred(RBMs.batch_size(rbm, v, h)) == (3,2)
    @test energy(rbm, v, h) ≈ [energy(rbm, v[:,1,j], h[:,i,1]) for i = 1:3, j = 1:2]

    v = bitrand(3,1,2)
    h = bitrand(2,3)
    @test size(energy(rbm, v, h)) == @inferred(RBMs.batch_size(rbm, v, h)) == (3,2)
    @test energy(rbm, v, h) ≈ [energy(rbm, v[:,1,j], h[:,i]) for i = 1:3, j = 1:2]

    v = bitrand(3,1)
    h = bitrand(2,3,2)
    @test size(energy(rbm, v, h)) == @inferred(RBMs.batch_size(rbm, v, h)) == (3,2)
    @test energy(rbm, v, h) ≈ [energy(rbm, v[:,1], h[:,i,j]) for i = 1:3, j = 1:2]

    v = bitrand(3,1,2)
    h = bitrand(2,3,2)
    @test size(energy(rbm, v, h)) == @inferred(RBMs.batch_size(rbm, v, h)) == (3,2)
    @test energy(rbm, v, h) ≈ [energy(rbm, v[:,1,j], h[:,i,j]) for i = 1:3, j = 1:2]

    v = bitrand(3,3,1)
    h = bitrand(2,3,2)
    @test size(energy(rbm, v, h)) == @inferred(RBMs.batch_size(rbm, v, h)) == (3,2)
    @test energy(rbm, v, h) ≈ [energy(rbm, v[:,i,1], h[:,i,j]) for i = 1:3, j = 1:2]

    v = bitrand(3,2,2)
    h = bitrand(2,3,1)
    @test_throws Any RBMs.batch_size(rbm, v, h)
    @test_throws Any energy(rbm, v, h)

    v = bitrand(3,1,2)
    h = bitrand(2,3,3)
    @test_throws Any RBMs.batch_size(rbm, v, h)
    @test_throws Any energy(rbm, v, h)
end

@testset "sample_v_from_v and sample_h_from_h on binary RBM" begin
    rbm = RBMs.BinaryRBM(randn(3,2), randn(2,3), zeros(3,2,2,3))
    v = bitrand(size(visible(rbm))..., 10^6)
    v = RBMs.sample_v_from_v(rbm, v)
    @test RBMs.batchmean(visible(rbm), v) ≈ RBMs.mean_from_inputs(visible(rbm)) rtol=0.1

    h = bitrand(size(hidden(rbm))...,  10^6)
    h = RBMs.sample_h_from_h(rbm, h)
    @test RBMs.batchmean(hidden(rbm), h) ≈ RBMs.mean_from_inputs(hidden(rbm)) rtol=0.1

    randn!(weights(rbm))
    h = RBMs.sample_h_from_v(rbm, v)
    μ = RBMs.mean_from_inputs(hidden(rbm), inputs_h_from_v(rbm, v))
    @test RBMs.batchmean(hidden(rbm), h) ≈ RBMs.batchmean(hidden(rbm), μ) rtol=0.1

    randn!(weights(rbm))
    v = RBMs.sample_v_from_h(rbm, h)
    μ = RBMs.mean_from_inputs(visible(rbm), inputs_v_from_h(rbm, h))
    @test RBMs.batchmean(visible(rbm), v) ≈ RBMs.batchmean(visible(rbm), μ) rtol=0.1
end

@testset "rbm convenience constructors" begin
    rbm = RBMs.BinaryRBM(randn(5), randn(3), randn(5,3))
    @test visible(rbm) isa RBMs.Binary
    @test hidden(rbm) isa RBMs.Binary
    @test size(weights(rbm)) == (5,3)

    rbm = RBMs.HopfieldRBM(randn(5), randn(3), rand(3), randn(5,3))
    @test visible(rbm) isa RBMs.Spin
    @test hidden(rbm) isa RBMs.Gaussian
    @test size(RBMs.weights(rbm)) == (5,3)
    @test all(hidden(rbm).γ .> 0)

    rbm = RBMs.HopfieldRBM(randn(5), randn(5,3))
    @test RBMs.visible(rbm) isa RBMs.Spin
    @test iszero(RBMs.mean_from_inputs(hidden(rbm)))
    @test RBMs.var_from_inputs(hidden(rbm)) == ones(size(hidden(rbm)))
    @test size(RBMs.weights(rbm)) == (5,3)
end

@testset "Binary-Binary RBM" begin
    rbm = RBMs.BinaryRBM(randn(5,2), randn(4,3), randn(5,2,4,3))
    v = rand(Bool, size(visible(rbm))..., 7)
    h = rand(Bool, size(hidden(rbm))..., 7)
    @test size(@inferred interaction_energy(rbm, v, h)) == (7,)
    @test size(@inferred energy(rbm, v, h)) == (7,)

    Ew = -[sum(v[i,j,b] * weights(rbm)[i,j,μ,ν] * h[μ,ν,b] for i=1:5, j=1:2, μ=1:4, ν=1:3) for b=1:7]
    @test interaction_energy(rbm, v, h) ≈ Ew
    @test energy(rbm, v, h) ≈ energy(visible(rbm), v) + energy(hidden(rbm), h) + Ew

    @test (@inferred energy(rbm, v[:,:,1], h[:,:,1])) isa Real
    @test (@inferred energy(rbm, v[:,:,1], h)) isa AbstractVector{<:Real}
    @test (@inferred energy(rbm, v, h[:,:,1])) isa AbstractVector{<:Real}

    for b = 1:7
        @test energy(rbm, v[:,:,b], h[:,:,b]) ≈ energy(rbm, v, h)[b]
        @test energy(rbm, v[:,:,b], h) ≈ [energy(rbm, v[:,:,b], h[:,:,k]) for k=1:7]
        @test energy(rbm, v, h[:,:,b]) ≈ [energy(rbm, v[:,:,k], h[:,:,b]) for k=1:7]
    end

    @test size(@inferred RBMs.sample_h_from_v(rbm, v)) == size(h)
    @test size(@inferred RBMs.sample_v_from_h(rbm, h)) == size(v)
    @test size(@inferred RBMs.sample_h_from_v(rbm, v[:,:,1])) == size(hidden(rbm))
    @test size(@inferred RBMs.sample_v_from_h(rbm, h[:,:,1])) == size(visible(rbm))
    for k = 1:3
        @test size(@inferred RBMs.sample_v_from_v(rbm, v; steps=k)) == size(v)
        @test size(@inferred RBMs.sample_h_from_h(rbm, h; steps=k)) == size(h)
        @test size(@inferred RBMs.sample_v_from_v(rbm, v[:,:,1]; steps=k)) == size(visible(rbm))
        @test size(@inferred RBMs.sample_h_from_h(rbm, h[:,:,1]; steps=k)) == size(hidden(rbm))
    end

    @test size(@inferred RBMs.mean_h_from_v(rbm, v)) == size(h)
    @test size(@inferred RBMs.mean_v_from_h(rbm, h)) == size(v)
    @test size(@inferred RBMs.mean_h_from_v(rbm, v[:,:,1])) == size(hidden(rbm))
    @test size(@inferred RBMs.mean_v_from_h(rbm, h[:,:,1])) == size(visible(rbm))

    @test size(@inferred RBMs.mode_h_from_v(rbm, v)) == size(h)
    @test size(@inferred RBMs.mode_v_from_h(rbm, h)) == size(v)
    @test size(@inferred RBMs.mode_h_from_v(rbm, v[:,:,1])) == size(hidden(rbm))
    @test size(@inferred RBMs.mode_v_from_h(rbm, h[:,:,1])) == size(visible(rbm))

    @test size(@inferred RBMs.var_h_from_v(rbm, v)) == size(h)
    @test size(@inferred RBMs.var_v_from_h(rbm, h)) == size(v)
    @test size(@inferred RBMs.var_h_from_v(rbm, v[:,:,1])) == size(hidden(rbm))
    @test size(@inferred RBMs.var_v_from_h(rbm, h[:,:,1])) == size(visible(rbm))

    @test size(@inferred free_energy(rbm, v)) == (7,)
    @test size(@inferred RBMs.reconstruction_error(rbm, v)) == (7,)
    @test (@inferred free_energy(rbm, v[:,:,1])) isa Real
    @test (@inferred RBMs.reconstruction_error(rbm, v[:,:,1])) isa Real

    @inferred RBMs.mirror(rbm)
    @test RBMs.mirror(rbm).visible == hidden(rbm)
    @test RBMs.mirror(rbm).hidden == visible(rbm)
    @test energy(RBMs.mirror(rbm), h, v) ≈ energy(rbm, v, h)
end

@testset "binary free energy" begin
    rbm = RBMs.BinaryRBM(randn(1), randn(1), randn(1,1))
    for v in [[0], [1]]
        @test -free_energy(rbm, v) ≈ logsumexp(-energy(rbm, v, h) for h in [[0], [1]])
    end

    rbm = RBMs.BinaryRBM(randn(7), randn(2), randn(7,2))
    v = bitrand(7)
    hs = [[0,0], [0,1], [1,0], [1,1]]
    @test -free_energy(rbm, v) ≈ logsumexp(-energy(rbm, v, h) for h in hs)
end

@testset "Gaussian-Gaussian RBM, 1-dimension" begin
    rbm = RBMs.RBM(
        RBMs.Gaussian(randn(1), rand(1) .+ 1),
        RBMs.Gaussian(randn(1), rand(1) .+ 1),
        randn(1, 1) / 1e2
    )

    @assert isposdef([
        visible(rbm).γ -weights(rbm);
        -weights(rbm)'  hidden(rbm).γ
    ])

    Z, ϵ = quadgk(x -> exp(-only(free_energy(rbm, [x;;]))), -Inf, Inf)
    @test RBMs.log_partition(rbm) ≈ log(Z)
end

@testset "Gaussian-Gaussian RBM, multi-dimensional" begin
    n = (10, 3)
    m = (7, 2)
    rbm = RBMs.RBM(
        RBMs.Gaussian(randn(n...), rand(n...) .+ 0.5),
        RBMs.Gaussian(randn(m...), rand(m...) .+ 0.5),
        randn(n..., m...) / (10 * prod(n) * prod(m))
    )

    N = length(visible(rbm))
    M = length(hidden(rbm))

    θ = [
        vec(visible(rbm).θ);
        vec(hidden(rbm).θ)
    ]
    γv = vec(visible(rbm).γ)
    γh = vec(hidden(rbm).γ)
    w = reshape(weights(rbm), length(visible(rbm)), length(hidden(rbm)))
    A = [diagm(γv) -w;
         -w'  diagm(γh)]

    v = randn(n..., 1)
    h = randn(m..., 1)
    x = [reshape(v, N, 1); reshape(h, M, 1)]

    @test energy(rbm, v, h) ≈ x' * A * x / 2 - θ' * x
    @test RBMs.log_partition(rbm) ≈ (
        (N + M)/2 * log(2π) + θ' * inv(A) * θ / 2 - logdet(A)/2
    )
    @test RBMs.log_likelihood(rbm, v) ≈ -free_energy(rbm, v) .- RBMs.log_partition(rbm)

    Ev = sum(visible(rbm).γ .* v.^2 ./ 2 - visible(rbm).θ .* v)
    Fv = sum(-(hidden(rbm).θ .+ inputs_h_from_v(rbm, v)).^2 ./ 2hidden(rbm).γ)
    @test only(free_energy(rbm, v)) ≈ Ev + Fv - sum(log.(2π ./ hidden(rbm).γ)) / 2

    gs = Zygote.gradient(rbm) do rbm
        RBMs.log_partition(rbm)
    end
    ∂θv = vec(only(gs).visible.θ)
    ∂θh = vec(only(gs).hidden.θ)
    ∂γv = vec(only(gs).visible.γ)
    ∂γh = vec(only(gs).hidden.γ)
    ∂w = reshape(only(gs).w, length(visible(rbm)), length(hidden(rbm)))

    Σ = inv(A) # covariances
    μ = Σ * θ  # means
    C = Σ + μ * μ' # non-centered second moments

    @test [∂θv; ∂θh] ≈ μ
    @test -2∂γv ≈ diag(C[1:N, 1:N]) # <vi^2>
    @test -2∂γh ≈ diag(C[(N + 1):end, (N + 1):end]) # <hμ^2>
    @test ∂w ≈ C[1:N, (N + 1):end] # <vi*hμ>
end

@testset "zero hidden units" begin
    rbm = BinaryRBM(randn(5), randn(0), randn(5,0))
    v = bitrand(5)
    @test free_energy(rbm, v) ≈ energy(rbm.visible, v)
end


text:
"onesite tdvp"
@with_kw struct Tdvp <: Algorithm
    tol::Float64 = Defaults.tol
    tolgauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
end

"""
    function timestep(psi, operator, dt, alg,envs = environments(psi,operator))

time evolves psi by timestep dt using algorithm alg
"""
function timestep(state::InfiniteMPS, H::Hamiltonian, timestep::Number,alg::Tdvp,envs::Cache=environments(state,H))

    temp_ACs = similar(state.AC);
    temp_CRs = similar(state.CR);

    @sync for (loc,(ac,c)) in enumerate(zip(state.AC,state.CR))
        @Threads.spawn begin
            ($temp_ACs[loc],convhist) = exponentiate(x->ac_prime(x,$loc,$state,$envs) ,-1im*$timestep,$ac,Lanczos(tol=alg.tol))
            convhist.converged==0 && @info "time evolving ac($loc) failed $(convhist.normres)"
        end

        @Threads.spawn begin
            ($temp_CRs[loc],convhist) = exponentiate(x->c_prime(x,$loc,$state,$envs) ,-1im*$timestep,$c,Lanczos(tol=alg.tol))
            convhist.converged==0 && @info "time evolving a($loc) failed $(convhist.normres)"
        end
    end

    for loc in 1:length(state)

        #find Al that best fits these new Acenter and centers
        QAc,_ = leftorth!(temp_ACs[loc],alg=TensorKit.QRpos())
        Qc,_ = leftorth!(temp_CRs[loc],alg=TensorKit.QRpos())
        @plansor temp_ACs[loc][-1 -2;-3] = QAc[-1 -2;1]*conj(Qc[-3;1])
    end

    nstate = InfiniteMPS(temp_ACs,state.CR[end]; tol = alg.tolgauge, maxiter = alg.maxiter)
    recalculate!(envs,nstate)
    nstate,envs
end

function timestep!(state::Union{FiniteMPS,MPSComoving}, H::Operator, timestep::Number,alg::Tdvp,envs=environments(state,H))
    #left to right
    for i in 1:(length(state)-1)
        (state.AC[i],convhist) = exponentiate(@closure(x->ac_prime(x,i,state,envs)),-1im*timestep/2,state.AC[i],Lanczos(tol=alg.tolgauge))
        (state.CR[i],convhist) = exponentiate(@closure(x->c_prime(x,i,state,envs)),1im*timestep/2,state.CR[i],Lanczos(tol=alg.tolgauge))

    end

    (state.AC[end],convhist) = exponentiate(@closure(x->ac_prime(x,length(state),state,envs)),-1im*timestep/2,state.AC[end],Lanczos(tol=alg.tolgauge))

    #right to left
    for i in length(state):-1:2
        (state.AC[i],convhist) = exponentiate(@closure(x->ac_prime(x,i,state,envs)),-1im*timestep/2,state.AC[i],Lanczos(tol=alg.tolgauge))
        (state.CR[i-1],convhist) = exponentiate(@closure(x->c_prime(x,i-1,state,envs)),1im*timestep/2,state.CR[i-1],Lanczos(tol=alg.tolgauge))
    end

    (state.AC[1],convhist) = exponentiate(@closure(x->ac_prime(x,1,state,envs)),-1im*timestep/2,state.AC[1],Lanczos(tol=alg.tolgauge))

    return state,envs
end

"twosite tdvp (works for finite mps's)"
@with_kw struct Tdvp2 <: Algorithm
    tol::Float64 = Defaults.tol
    tolgauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    trscheme = truncerr(1e-3)
end

#twosite tdvp for finite mps
function timestep!(state::Union{FiniteMPS,MPSComoving}, H::Operator, timestep::Number,alg::Tdvp2,envs=environments(state,H);rightorthed=false)
    #left to right
    for i in 1:(length(state)-1)
        ac2 = _transpose_front(state.AC[i])*_transpose_tail(state.AR[i+1])

        (nac2,convhist) = exponentiate(@closure(x->ac2_prime(x,i,state,envs)),-1im*timestep/2,ac2,Lanczos())
        (nal,nc,nar) = tsvd(nac2,trunc=alg.trscheme, alg=TensorKit.SVD())

        state.AC[i] = (nal,complex(nc))
        state.AC[i+1] = (complex(nc),_transpose_front(nar))

        if(i!=(length(state)-1))
            (state.AC[i+1],convhist) = exponentiate(@closure(x->ac_prime(x,i+1,state,envs)),1im*timestep/2,state.AC[i+1],Lanczos())
        end

    end

    #right to left
    for i in length(state):-1:2
        ac2 = _transpose_front(state.AL[i-1])*_transpose_tail(state.AC[i])

        (nac2,convhist) = exponentiate(@closure(x->ac2_prime(x,i-1,state,envs)),-1im*timestep/2,ac2,Lanczos())
        (nal,nc,nar) = tsvd(nac2,trunc=alg.trscheme,alg=TensorKit.SVD())

        state.AC[i-1] = (nal,complex(nc))
        state.AC[i] = (complex(nc),_transpose_front(nar));

        if(i!=2)
            (state.AC[i-1],convhist) = exponentiate(@closure(x->ac_prime(x,i-1,state,envs)),1im*timestep/2,state.AC[i-1],Lanczos())
        end
    end

    return state,envs
end

#copying version
timestep(state::Union{FiniteMPS,MPSComoving},H,timestep,alg::Union{Tdvp,Tdvp2},envs=environments(state,H)) = timestep!(copy(state),H,timestep,alg,envs)


text:
module Box

using Crayons

export BLACK_FG,
RED_FG,
GREEN_FG,
YELLOW_FG,
BLUE_FG,
MAGENTA_FG,
CYAN_FG,
LIGHT_GRAY_FG,
DEFAULT_FG,
DARK_GRAY_FG,
LIGHT_RED_FG,
LIGHT_GREEN_FG,
LIGHT_YELLOW_FG,
LIGHT_BLUE_FG,
LIGHT_MAGENTA_FG,
LIGHT_CYAN_FG,
WHITE_FG,
BLACK_BG,
RED_BG,
GREEN_BG,
YELLOW_BG,
BLUE_BG,
MAGENTA_BG,
CYAN_BG,
LIGHT_GRAY_BG,
DEFAULT_BG,
DARK_GRAY_BG,
LIGHT_RED_BG,
LIGHT_GREEN_BG,
LIGHT_YELLOW_BG,
LIGHT_BLUE_BG,
LIGHT_MAGENTA_BG,
LIGHT_CYAN_BG,
WHITE_BG,
BOLD,
FAINT,
ITALICS,
UNDERLINE,
BLINK,
NEGATIVE,
CONCEAL,
STRIKETHROUGH

const BLACK_FG         = Crayon(foreground = :black         )
const RED_FG           = Crayon(foreground = :red           )
const GREEN_FG         = Crayon(foreground = :green         )
const YELLOW_FG        = Crayon(foreground = :yellow        )
const BLUE_FG          = Crayon(foreground = :blue          )
const MAGENTA_FG       = Crayon(foreground = :magenta       )
const CYAN_FG          = Crayon(foreground = :cyan          )
const LIGHT_GRAY_FG    = Crayon(foreground = :light_gray    )
const DEFAULT_FG       = Crayon(foreground = :default       )
const DARK_GRAY_FG     = Crayon(foreground = :dark_gray     )
const LIGHT_RED_FG     = Crayon(foreground = :light_red     )
const LIGHT_GREEN_FG   = Crayon(foreground = :light_green   )
const LIGHT_YELLOW_FG  = Crayon(foreground = :light_yellow  )
const LIGHT_BLUE_FG    = Crayon(foreground = :light_blue    )
const LIGHT_MAGENTA_FG = Crayon(foreground = :light_magenta )
const LIGHT_CYAN_FG    = Crayon(foreground = :light_cyan    )
const WHITE_FG         = Crayon(foreground = :white         )

const BLACK_BG         = Crayon(background = :black         )
const RED_BG           = Crayon(background = :red           )
const GREEN_BG         = Crayon(background = :green         )
const YELLOW_BG        = Crayon(background = :yellow        )
const BLUE_BG          = Crayon(background = :blue          )
const MAGENTA_BG       = Crayon(background = :magenta       )
const CYAN_BG          = Crayon(background = :cyan          )
const LIGHT_GRAY_BG    = Crayon(background = :light_gray    )
const DEFAULT_BG       = Crayon(background = :default       )
const DARK_GRAY_BG     = Crayon(background = :dark_gray     )
const LIGHT_RED_BG     = Crayon(background = :light_red     )
const LIGHT_GREEN_BG   = Crayon(background = :light_green   )
const LIGHT_YELLOW_BG  = Crayon(background = :light_yellow  )
const LIGHT_BLUE_BG    = Crayon(background = :light_blue    )
const LIGHT_MAGENTA_BG = Crayon(background = :light_magenta )
const LIGHT_CYAN_BG    = Crayon(background = :light_cyan    )
const WHITE_BG         = Crayon(background = :white         )

const BOLD             = Crayon(bold          = true)
const FAINT            = Crayon(faint         = true)
const ITALICS          = Crayon(italics       = true)
const UNDERLINE        = Crayon(underline     = true)
const BLINK            = Crayon(blink         = true)
const NEGATIVE         = Crayon(negative      = true)
const CONCEAL          = Crayon(conceal       = true)
const STRIKETHROUGH    = Crayon(strikethrough = true)

end