text:
#### Turbulence model kernels

"""
    thermo_variables(ts::ThermodynamicState)

A NamedTuple of thermodynamic variables,
computed from the given thermodynamic state.
"""
function thermo_variables(ts::ThermodynamicState)
    return (
        θ_dry = dry_pottemp(ts),
        θ_liq = liquid_ice_pottemp(ts),
        q_tot = total_specific_humidity(ts),
        T = air_temperature(ts),
        R_m = gas_constant_air(ts),
        q_vap = vapor_specific_humidity(ts),
        q_liq = liquid_specific_humidity(ts),
        q_ice = ice_specific_humidity(ts),
    )
end

"""
    compute_buoyancy_gradients(
        m::AtmosModel{FT},
        args,
        ts_gm,
        ts_en
    ) where {FT}
Returns the environmental buoyancy gradient following Tan et al. (JAMES, 2018)
and the effective environmental static stability following
Lopez-Gomez et al. (JAMES, 2020), given:
 - `m`, an `AtmosModel`
 - `args`, top-level arguments
 - `ts_gm`, grid-mean thermodynamic state
 - `ts_en`, environment thermodynamic state
"""
function compute_buoyancy_gradients(
    m::AtmosModel{FT},
    args,
    ts_gm,
    ts_en,
) where {FT}
    @unpack state, aux, diffusive = args
    # Alias convention:
    gm = state
    en_dif = diffusive.turbconv.environment
    N_up = n_updrafts(turbconv_model(m))
    param_set = parameter_set(m)

    _grav::FT = grav(param_set)
    _R_d::FT = R_d(param_set)
    _R_v::FT = R_v(param_set)
    ε_v::FT = 1 / molmass_ratio(param_set)
    p = air_pressure(ts_gm)

    q_tot_en = total_specific_humidity(ts_en)
    θ_liq_en = liquid_ice_pottemp(ts_en)
    lv = latent_heat_vapor(ts_en)
    T = air_temperature(ts_en)
    Π = exner(ts_en)
    q_liq = liquid_specific_humidity(ts_en)
    _cp_m = cp_m(ts_en)
    θ_virt = virtual_pottemp(ts_en)

    (ts_dry, ts_cloudy, cloud_frac) =
        compute_subdomain_statistics(m, args, ts_gm, ts_en)
    cloudy = thermo_variables(ts_cloudy)
    dry = thermo_variables(ts_dry)

    prefactor = _grav * (_R_d * gm.ρ / p * Π)

    ∂b∂θl_dry = prefactor * (1 + (ε_v - 1) * dry.q_tot)
    ∂b∂qt_dry = prefactor * dry.θ_liq * (ε_v - 1)

    if cloud_frac > FT(0)
        num =
            prefactor *
            (1 + ε_v * (1 + lv / _R_v / cloudy.T) * cloudy.q_vap - cloudy.q_tot)
        den = 1 + lv * lv / _cp_m / _R_v / cloudy.T / cloudy.T * cloudy.q_vap
        ∂b∂θl_cloudy = num / den
        ∂b∂qt_cloudy =
            (lv / _cp_m / cloudy.T * ∂b∂θl_cloudy - prefactor) * cloudy.θ_dry
    else
        ∂b∂θl_cloudy = FT(0)
        ∂b∂qt_cloudy = FT(0)
    end

    ∂b∂θl = (cloud_frac * ∂b∂θl_cloudy + (1 - cloud_frac) * ∂b∂θl_dry)
    ∂b∂qt = (cloud_frac * ∂b∂qt_cloudy + (1 - cloud_frac) * ∂b∂qt_dry)

    # Partial buoyancy gradients
    ∂b∂z_θl = en_dif.∇θ_liq[3] * ∂b∂θl
    ∂b∂z_qt = en_dif.∇q_tot[3] * ∂b∂qt
    ∂b∂z = ∂b∂z_θl + ∂b∂z_qt

    # Computation of buoyancy frequency based on θ_lv
    ∂θvl∂θ_liq = 1 + (ε_v - 1) * q_tot_en
    ∂θvl∂qt = (ε_v - 1) * θ_liq_en
    # apply chain-rule
    ∂θvl∂z = ∂θvl∂θ_liq * en_dif.∇θ_liq[3] + ∂θvl∂qt * en_dif.∇q_tot[3]

    ∂θv∂θvl = exp(lv * q_liq / _cp_m / T)
    λ_stb = cloud_frac

    Nˢ_eff =
        _grav / θ_virt *
        ((1 - λ_stb) * en_dif.∇θv[3] + λ_stb * ∂θvl∂z * ∂θv∂θvl)
    return ∂b∂z, Nˢ_eff
end;

"""
    ∇Richardson_number(
        ∂b∂z::FT,
        Shear²::FT,
        minval::FT,
        Ri_c::FT,
    ) where {FT}

Returns the gradient Richardson number, given:
 - `∂b∂z`, the vertical buoyancy gradient
 - `Shear²`, the squared vertical gradient of horizontal velocity
 - `maxval`, maximum value of the output, typically the critical Ri number
"""
function ∇Richardson_number(
    ∂b∂z::FT,
    Shear²::FT,
    minval::FT,
    Ri_c::FT,
) where {FT}
    return min(∂b∂z / max(Shear², minval), Ri_c)
end;

"""
    turbulent_Prandtl_number(
        Pr_n::FT,
        Grad_Ri::FT,
        ω_pr::FT
    ) where {FT}

Returns the turbulent Prandtl number, given:
 - `Pr_n`, the turbulent Prandtl number under neutral conditions
 - `Grad_Ri`, the gradient Richardson number
"""
function turbulent_Prandtl_number(Pr_n::FT, Grad_Ri::FT, ω_pr::FT) where {FT}
    if Grad_Ri > FT(0)
        factor =
            2 * Grad_Ri /
            (1 + ω_pr * Grad_Ri - sqrt((1 + ω_pr * Grad_Ri)^2 - 4 * Grad_Ri))
    else
        factor = FT(1)
    end
    return Pr_n * factor
end;


text:
########################################################################################################
# This file runs a number of model analyses from Budolfson et al. (2019) given user-specified settings.
########################################################################################################

# Activate the project for the paper and make sure all packages we need
# are installed.
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# Load required Julia packages.
using NLopt
using CSVFiles

# Load required files.
include("create_rice.jl")
include("create_fund.jl")
include("helper_functions.jl")

#------------------------------------------------------------------------------------------------------
# Model Parameters To Modify.
#------------------------------------------------------------------------------------------------------

# Pure rate of time preference.
ρ = 0.008

# Elasticity of marginal utility of consumption.
η = 1.5

# Should the RICE model be run with a social welfare function that uses Negish weights (true = remove Negishi weights).
# *Note: if using Negishi weights, the model uses the default RICE settings of η=1.5 and ρ=1.5%.
remove_negishi = true

#------------------------------------------------------------------------------------------------------
# Choices For Model Versions To Run.
#------------------------------------------------------------------------------------------------------

# Run an optiization with RICE using global carbon prices?
rice_cost_minimization = true

# Run an optization with RICE using regional carbon prices?
rice_utilitarian = true

# Run an optiization with RICE using global carbon prices?
fund_cost_minimization = false

# Run an optization with RICE using regional carbon prices?
fund_utilitarian = false

#------------------------------------------------------------------------------------------------------
# Optimization Settings to Modify.
#------------------------------------------------------------------------------------------------------

# Name of folder to save this set of model runs in (a folder will be created with this name).
results_folder = "My Results"

# Number of model periods to optimize over for RICE (after which model assumes full decarbonization).
# NOTE: FUND does not have a backstop price and optimizes from 2010-2200 by default, and then assumes a constant carbon tax.
n_opt_periods = 30

# Optimization algorithm (the type should be a Symbol, e.g. :LN_SBPLX). See options at http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms
optimization_algorithm = :LN_SBPLX

# Maximum time in seconds to run each model (NOTE: FUND takes much longer to optimize than RICE).
stop_time_rice = 500
stop_time_fund = 7500

# Relative tolerance criteria for convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
tolerance_rice = 1e-10
tolerance_fund = 1e-10

#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
# Run Everything and Save Key Results (*no need to modify code below this line).
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
println("Starting Analysis.")
println()

# Load and run an instance of RICE just to extract the backstop prices (needed for multiple analyses).
backstop_rice = create_rice(ρ, η, remove_negishi)
run(backstop_rice)
backstop_prices = backstop_rice[:emissions, :pbacktime] .* 1000
output_directory = joinpath(@__DIR__, "../", "results", results_folder, "rice_BAU")
mkpath(output_directory)
save(joinpath(output_directory, "Emissions.csv"), DataFrame(backstop_rice[:emissions, :EIND]))
save(joinpath(output_directory, "PerCapitaConsumption.csv"), DataFrame(backstop_rice[:neteconomy, :CPC]))
save(joinpath(output_directory, "Population.csv"), DataFrame(backstop_rice[:welfare, :pop]))
LandUse = load(joinpath(@__DIR__, "../", "data", "RICELandUse.csv")) |> DataFrame
save(joinpath(output_directory, "LandUse.csv"), LandUse)

#------------------------------------------------------------------------------------------------------
# Run RICE Cost-minimization Optimization.
#------------------------------------------------------------------------------------------------------
if rice_cost_minimization == true

    println("Starting RICE cost-minimization optimization...")

    # Optimize model.
    opt_output_rice_cost_minimization, opt_emissions_rice_cost_minimization, opt_mitigation_rice_cost_minimization, opt_tax_rice_cost_minimization, opt_model_rice_cost_minimization, convergence_rice_cost_minimization = optimize_rice(optimization_algorithm, n_opt_periods, stop_time_rice, tolerance_rice, backstop_prices, run_utilitarian=false, ρ=ρ, η=η, remove_negishi=remove_negishi)

    # Create folder to store some key results.
    output_directory = joinpath(@__DIR__, "../", "results", results_folder, "rice_costmin")
    mkpath(output_directory)

    # Save optimal CO₂ mitigation, carbon tax, per capita consumption, and global temperature anomaly.
    save(joinpath(output_directory, "Emissions.csv"), DataFrame(opt_emissions_rice_cost_minimization))
    save(joinpath(output_directory, "MitigationRate.csv"), DataFrame(opt_mitigation_rice_cost_minimization))
    save(joinpath(output_directory, "Carbon Tax.csv"), DataFrame(global_carbon_tax=opt_tax_rice_cost_minimization))
    save(joinpath(output_directory, "Temperature.csv"), DataFrame(temp=opt_model_rice_cost_minimization[:climatedynamics, :TATM]))
    save(joinpath(output_directory, "PerCapitaConsumption.csv"), DataFrame(opt_model_rice_cost_minimization[:neteconomy, :CPC]))
    save(joinpath(output_directory, "GDP.csv"), DataFrame(opt_model_rice_cost_minimization[:neteconomy, :YNET]))

    println("Optimization convergence result: ", convergence_rice_cost_minimization)
    println("RICE cost-minimization optimization complete.")
    println()
end

#------------------------------------------------------------------------------------------------------
# Run RICE Utilitarian Optimization.
#------------------------------------------------------------------------------------------------------
if rice_utilitarian == true

    println("Starting RICE utilitarian optimization...")

    # Optimize model.
    opt_output_rice_utilitarian, opt_emissions_rice_utilitarian, opt_mitigation_rice_utilitarian, opt_tax_rice_utilitarian, opt_model_rice_utilitarian, convergence_rice_utilitarian = optimize_rice(optimization_algorithm, n_opt_periods, stop_time_rice, tolerance_rice, backstop_prices, run_utilitarian=true, ρ=ρ, η=η, remove_negishi=remove_negishi)

    # Create folder to store some key results.
    output_directory = joinpath(@__DIR__, "../", "results", results_folder, "rice_utilitarian")
    mkpath(output_directory)

    # Save optimal CO₂ mitigation, carbon tax, per capita consumption, and global temperature anomaly.
    save(joinpath(output_directory, "MitigationRate.csv"), DataFrame(opt_mitigation_rice_utilitarian))
    save(joinpath(output_directory, "Emissions.csv"), DataFrame(opt_emissions_rice_utilitarian))
    save(joinpath(output_directory, "Carbon Tax.csv"), DataFrame(opt_tax_rice_utilitarian))
    save(joinpath(output_directory, "Temperature.csv"), DataFrame(temp=opt_model_rice_utilitarian[:climatedynamics, :TATM]))
    save(joinpath(output_directory, "GDP.csv"), DataFrame(opt_model_rice_utilitarian[:neteconomy, :YNET]))
    save(joinpath(output_directory, "PerCapitaConsumption.csv"), DataFrame(opt_model_rice_utilitarian[:neteconomy, :CPC]))

    println("Optimization convergence result: ", convergence_rice_utilitarian)
    println("RICE utilitarian optimization complete.")
    println()
end

#------------------------------------------------------------------------------------------------------
# Run FUND Cost-minimization Optimization.
#------------------------------------------------------------------------------------------------------
if fund_cost_minimization == true

    println("Starting FUND cost-minimization optimization...")

    # Optimize model.
    opt_output_fund_cost_minimization, opt_mitigation_fund_cost_minimization, opt_tax_fund_cost_minimization, opt_model_fund_cost_minimization, convergence_fund_cost_minimization = optimize_fund(optimization_algorithm, 20, stop_time_fund, tolerance_fund, run_utilitarian=false, ρ=ρ, η=η, welfare_year=2010, end_year=2300)

    # Create folder to store some key results.
    output_directory = joinpath(@__DIR__, "../", "results", results_folder, "fund_costmin")
    mkpath(output_directory)

    # Save optimal CO₂ mitigation, carbon tax, per capita consumption, and global temperature anomaly.
    save(joinpath(output_directory, "MitigationRate.csv"), DataFrame(opt_mitigation_fund_cost_minimization))
    # Just save one column of global tax array (each column is identical for all regions).
    save(joinpath(output_directory, "Carbon Tax.csv"), DataFrame(global_carbon_tax=opt_tax_fund_cost_minimization[:,1]))
    save(joinpath(output_directory, "Temperature.csv"), DataFrame(temperature=opt_model_fund_cost_minimization[:climatedynamics, :temp]))
    save(joinpath(output_directory, "PerCapitaConsumption.csv"), DataFrame(opt_model_fund_cost_minimization[:fund_welfare, :cpc]))

    println("Optimization convergence result: ", convergence_fund_cost_minimization)
    println("FUND cost_minimization optimization complete.")
    println()
end

#------------------------------------------------------------------------------------------------------
# Run FUND Utilitarian Optimization.
#------------------------------------------------------------------------------------------------------
if fund_utilitarian == true

    println("Starting FUND utilitarian optimization...")

    # Optimize model.
    opt_output_fund_utilitarian, opt_mitigation_fund_utilitarian, opt_tax_fund_utilitarian, opt_model_fund_utilitarian, convergence_fund_utilitarian = optimize_fund(optimization_algorithm, 20, stop_time_fund, tolerance_fund, run_utilitarian=true, ρ=ρ, η=η, welfare_year=2010, end_year=2300)

    # Create folder to store some key results.
    output_directory = joinpath(@__DIR__, "../", "results", results_folder, "fund_utilitarian")
    mkpath(output_directory)

    # Save optimal CO₂ mitigation, carbon tax, per capita consumption, and global temperature anomaly.
    save(joinpath(output_directory, "MitigationRate.csv"), DataFrame(opt_mitigation_fund_utilitarian))
    save(joinpath(output_directory, "Carbon Tax.csv"), DataFrame(opt_tax_fund_utilitarian))
    save(joinpath(output_directory, "Temperature.csv"), DataFrame(temperature=opt_model_fund_utilitarian[:climatedynamics, :temp]))
    save(joinpath(output_directory, "PerCapitaConsumption.csv"), DataFrame(opt_model_fund_utilitarian[:fund_welfare, :cpc]))

    println("Optimization convergence result: ", convergence_fund_utilitarian)
    println("FUND utilitiarian optimization complete.")
    println()
end

println("Analysis Complete.")

ext: jl

max_stars_repo_name: Environment-Research/Utilitarianism

max_issues_repo_name: Environment-Research/Utilitarianism

max_forks_repo_name: Environment-Research/Utilitarianism


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type TankL
	TankL()=begin
		PP=outers.PP
		NComp=outers.NComp
		[
			:(LI = Levelpercent_Initial),
			:(OutletLiquid.T = Temperature_Initial),
			:(OutletLiquid.z(1:NComp - 1) = Composition_Initial(1:NComp - 1)/sum(Composition_Initial)),
		],
		[
			"Initial level Percent","Initial Outlet Liquid Temperature","Initial Outlet Liquid Composition Normalized",
		],
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components",
				:Lower=>1
			)),
			acceleration(Dict{Symbol,Any}(
				:Brief=>"Gravity Acceleration",
				:Default=>9.81,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Low Flow",
				:Default=>1E-6,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"No Flow",
				:Default=>0,
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Constant for K factor pressure drop",
				:Unit=>"mol/(s*(Pa^0.5))",
				:Default=>1,
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"K factor for pressure drop",
				:Lower=>1E-8,
				:Default=>2
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Normal Flow",
				:Valid=>["on", "off"],
				:Default=>"on",
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Initial liquid height in Percent",
				:Default=>0.70
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Initial Liquid Temperature",
				:Default=>330
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Initial Composition",
				:Default=>0.10
			)),(NComp)),
			VesselVolume(Dict{Symbol,Any}(
				:Brief=>"Vessel Geometry",
				:Symbol=>" "
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed Stream",
				:PosX=>0.22,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.43,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat Duty",
				:PosX=>0.735,
				:PosY=>1,
				:Protected=>true,
				:Symbol=>"Q_{in}"
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the Vessel",
				:Protected=>true
			)),(NComp)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup in the Vessel",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume",
				:Protected=>true
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Static head at the bottom of the tank",
				:Protected=>true,
				:Symbol=>"P_{static}^{Liquid}"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature Indicator",
				:PosX=>0.525,
				:PosY=>0,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Pressure Indicator",
				:PosX=>0.368,
				:PosY=>0,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Level Indicator",
				:PosX=>1,
				:PosY=>0.6,
				:Protected=>true
			)),
			[
				:(diff(TotalHoldup)=Inlet.F*Inlet.z - OutletLiquid.F*OutletLiquid.z),
				:(diff(E) = Inlet.F*Inlet.h - OutletLiquid.F*OutletLiquid.h + InletQ),
				:(E = sum(TotalHoldup)*OutletLiquid.h),
				:(Pstatic = PP.LiquidDensity(OutletLiquid.T, Inlet.P, OutletLiquid.z) * Gconst * Geometry.Level),
				:(Inlet.P + Pstatic = OutletLiquid.P),
				:(vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(TotalHoldup = OutletLiquid.z*sum(TotalHoldup)),
				:(Geometry.Vfilled = sum(TotalHoldup) * vL),
				:(TI * "K" = OutletLiquid.T),
				:(PI * "atm" = OutletLiquid.P),
				:(LI*Geometry.Vtotal= Geometry.Vfilled),
			],
			[
				"Component Molar Balance","Energy Balance","Energy Holdup","Static Head","Mechanical Equilibrium","Liquid Volume","Molar Holdup","Liquid Level","Temperature indicator","Pressure indicator","Level indicator",
			],
			[:PP,:NComp,:Gconst,:low_flow,:zero_flow,:KfConst,:Kfactor,:NormalFlow,:Levelpercent_Initial,:Temperature_Initial,:Composition_Initial,],
			[:Geometry,:Inlet,:OutletLiquid,:InletQ,:TotalHoldup,:E,:vL,:Pstatic,:TI,:PI,:LI,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Gconst::acceleration
	low_flow::flow_mol
	zero_flow::flow_mol
	KfConst::DanaReal
	Kfactor::positive
	NormalFlow::DanaSwitcher
	Levelpercent_Initial::positive
	Temperature_Initial::temperature
	Composition_Initial::Array{fraction}
	Geometry::VesselVolume
	Inlet::stream
	OutletLiquid::liquid_stream
	InletQ::power
	TotalHoldup::Array{mol}
	E::energy
	vL::volume_mol
	Pstatic::pressure
	TI::control_signal
	PI::control_signal
	LI::control_signal
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export TankL
function set(in::TankL)
	Gconst = 9.81 * "m/(s^2)"
	 low_flow = 1E-6 * "kmol/h"
	 zero_flow = 0 * "kmol/h"
	 KfConst = 1*"mol/(s*(Pa^0.5))"

end
function initial(in::TankL)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function setEquationFlow(in::TankL)
	#
#switch NormalFlow
#
#case "on":
#	Inlet.F = Kfactor *sqrt(Pstatic)*KfConst;
#
#	when Inlet.F < low_flow switchto "off";
#
#case "off":
#	Inlet.F = zero_flow;
#
#	when Inlet.P > OutletLiquid.P switchto "on";
#
#end
#
addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
end
function atributes(in::TankL,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/TankL"
	fields[:Brief]="Model of a Tank."
	fields[:Info]="== ASSUMPTIONS ==
* liquid phase only;

== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;

== SPECIFY ==
* the Inlet stream;
* the OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the TankL model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the Tank temperature (OutletLiquid.T);
* Levelpercent_Initial : the Tank liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
"
	drive!(fields,_)
	return fields
end
TankL(_::Dict{Symbol,Any})=begin
	newModel=TankL()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module GaussQuadrature

# October 2013 by Bill McLean, School of Maths and Stats,
# The University of New South Wales.
#
# Based on earlier Fortran codes
#
# gaussq.f original version 20 Jan 1975 from Stanford
# gaussq.f modified 21 Dec by Eric Grosse
# gaussquad.f95 Nov 2005 by Bill Mclean
#
# This module provides functions to compute the abscissae x[j] and
# weights w[j] for the classical Gauss quadrature rules, including
# the Radau and Lobatto variants.  Thus, the sum
#
#           n
#          ---
#          \
#           |  w[j] f(x[j])
#          /
#          ---
#          j=1
#
# approximates
#
#          / hi
#          |
#          |    f(x) w(x) dx
#          |
#          / lo
#
# where the weight function w(x) and interval lo < x < hi are as shown
# in the table below.
#
# Name                      Interval     Weight Function
#
# Legendre                 -1 < x < 1          1
# Chebyshev (first kind)   -1 < x < 1     1 / sqrt(1-x^2)
# Chebyshev (second kind)  -1 < x < 1       sqrt(1-x^2)
# Jacobi                   -1 < x < 1   (1-x)^alpha (1+x)^beta
#
# For the Jacobi and Laguerre rules we require alpha > -1 and
# beta > -1, so that the weight function is integrable.
#
# Use the endpt argument to include one or both of the end points
# of the interval of integration as an abscissa in the quadrature
# rule, as follows.
#
# endpt = neither   Default      a < x[j] < b, j = 1:n.
# endpt = left      Left Radau   a = x[1] < x[j] < b, j = 2:n.
# endpt = right     Right Radau  a < x[j] < x[n] = b, j = 1:n-1.
# endpt = both      Lobatto      a = x[1] < x[j] < x[n] = b, j = 2:n-1.
#
# These labels make up an enumeration of type EndPt.
#
# The code uses the Golub and Welsch algorithm, in which the abscissae
# x[j] are the eigenvalues of a symmetric tridiagonal matrix whose
# entries depend on the coefficients in the 3-term recurrence relation
# for the othonormal polynomials generated by the weighted inner product.
#
# References:
#
#   1.  Golub, G. H., and Welsch, J. H., Calculation of Gaussian
#       quadrature rules, Mathematics of Computation 23 (April,
#       1969), pp. 221-230.
#   2.  Golub, G. H., Some modified matrix eigenvalue problems,
#       Siam Review 15 (april, 1973), pp. 318-334 (section 7).
#   3.  Stroud and Secrest, Gaussian Quadrature Formulas, Prentice-
#       Hall, Englewood Cliffs, N.J., 1966.

"""
Enumeration type used to specify which endpoints of the integration
interval should be included amongst the quadrature points: neither,
left, right or both.
"""
struct EndPt
    label :: Char
end

const neither = EndPt('N')
const left    = EndPt('L')
const right   = EndPt('R')
const both    = EndPt('B')

# Maximum number of QL iterations used by steig!.
# You might need to increase this.
maxiterations = Dict(Float32 => 30, Float64 => 30, BigFloat => 40)

"""
x, w = legendre(T, n, endpt=neither)
Returns points x and weights w for the n-point Gauss-Legendre rule
for the interval -1 < x < 1 with weight function w(x) = 1.
Use endpt=left, right, both for the left Radau, right Radau, Lobatto
rules.
"""
function legendre(::Type{T}, n::Integer, endpt::EndPt=neither) where {T<:AbstractFloat}
    a, b, muzero = legendre_coeff(T, n, endpt)
    return custom_gauss_rule(-one(T), one(T), a, b, muzero, endpt)
end

"""
x, w = legendre(n, endpt=neither)
Convenience function with type T = Float64.
"""
legendre(n, endpt=neither) = legendre(Float64, n, endpt)

function legendre_coeff(::Type{T}, n::Integer, endpt::EndPt) where {T<:AbstractFloat}
    muzero = convert(T, 2.0)
    a = zeros(T, n)
    b = zeros(T, n)
    for i = 1:n
        b[i] = i / sqrt(convert(T, 4*i^2-1))
    end
    return a, b, muzero
end

"""
x, w = chebyshev(T, n, kind=1, endpt=neither)
Returns points x and weights w for the n-point Gauss-Chebyshev rule
for the interval -1 < x < 1 with weight function
    w(x) = 1 / sqrt(1-x^2) if kind=1
    w(x) = sqrt(1-x^2)     if kind=2.
Use endpt=left, right, both for the left Radau, right Radau, Lobatto
rules.
"""
function chebyshev(::Type{T}, n::Integer, kind::Integer=1, endpt::EndPt=neither) where {T<:AbstractFloat}
    a, b, muzero = chebyshev_coeff(T, n, kind, endpt)
    return custom_gauss_rule(-one(T), one(T), a, b, muzero, endpt)
end

"""
x, w = chebyshev(n, kind=1, endpt=neither)
Convenience function with type T = Float64.
"""
chebyshev(n, kind=1, endpt=neither) = chebyshev(Float64, n, kind, endpt)

function chebyshev_coeff(::Type{T}, n::Integer, kind::Integer, endpt::EndPt) where {T<:AbstractFloat}
    muzero = convert(T, pi)
    half = convert(T, 0.5)
    a = zeros(T, n)
    b = fill(half, n)
    if kind == 1
        b[1] = sqrt(half)
    elseif kind == 2
        muzero /= 2
    else
        error("Unsupported value for kind")
    end
    return a, b, muzero
end

function custom_gauss_rule(lo::T, hi::T,
         a::Array{T,1}, b::Array{T,1}, muzero::T, endpt::EndPt,
         maxits::Integer=maxiterations[T]) where {T<:AbstractFloat}
    #
    # On entry:
    #
    # a, b hold the coefficients (as given, for instance, by
    # legendre_coeff) in the three-term recurrence relation
    # for the orthonormal polynomials p_0, p_1, p_2, ... , that is,
    #
    #    b[j] p (x) = (x-a[j]) p   (x) - b[j-1] p   (x).
    #          j                j-1              j-2
    #
    # muzero holds the zeroth moment of the weight function, that is
    #
    #              / hi
    #             |
    #    muzero = | w(x) dx.
    #             |
    #             / lo
    #
    # On return:
    #
    # x, w hold the points and weights.
    #
    n = length(a)
    @assert length(b) == n
    if endpt == left
        if n == 1
            a[1] = lo
        else
            a[n] = solve(n, lo, a, b) * b[n-1]^2 + lo
        end
    elseif endpt == right
        if n == 1
            a[1] = hi
        else
            a[n] = solve(n, hi, a, b) * b[n-1]^2 + hi
        end
    elseif endpt == both
        if n == 1
            error("Must have at least two points for both ends.")
        end
        g = solve(n, lo, a, b)
        t1 = ( hi - lo ) / ( g - solve(n, hi, a, b) )
        b[n-1] = sqrt(t1)
        a[n] = lo + g * t1
    end
    w = zero(a)
    steig!(a, b, w, maxits)
    for i = 1:n
        w[i] = muzero * w[i]^2
    end
    idx = sortperm(a)
    return a[idx], w[idx]
end

function solve(n::Integer, shift::T, a::Array{T,1}, b::Array{T,1}) where {T<:AbstractFloat}
    #
    # Perform elimination to find the nth component s = delta[n]
    # of the solution to the nxn linear system
    #
    #     ( J_n - shift I_n ) delta = e_n,
    #
    # where J_n is the symmetric tridiagonal matrix with diagonal
    # entries a[i] and off-diagonal entries b[i], and e_n is the nth
    # standard basis vector.
    #
    t = a[1] - shift
    for i = 2:n-1
        t = a[i] - shift - b[i-1]^2 / t
    end
    return one(t) / t
end

function steig!(d::Array{T,1}, e::Array{T,1}, z::Array{T,1}, maxits::Integer) where {T<:AbstractFloat}
    #
    # Finds the eigenvalues and first components of the normalised
    # eigenvectors of a symmetric tridiagonal matrix by the implicit
    # QL method.
    #
    # d[i]   On entry, holds the ith diagonal entry of the matrix.
    #        On exit, holds the ith eigenvalue.
    #
    # e[i]   On entry, holds the [i+1,i] entry of the matrix for
    #        i = 1, 2, ..., n-1.  (The value of e[n] is not used.)
    #        On exit, e is overwritten.
    #
    # z[i]   On exit, holds the first component of the ith normalised
    #        eigenvector associated with d[i].
    #
    # maxits The maximum number of QL iterations.
    #
    # Martin and Wilkinson, Numer. Math. 12: 377-383 (1968).
    # Dubrulle, Numer. Math. 15: 450 (1970).
    # Handbook for Automatic Computation, Vol ii, Linear Algebra,
    #        pp. 241-248, 1971.
    #
    # This is a modified version of the Eispack routine imtql2.
    #
    n = length(z)
    z[1] = 1
    z[2:n] = 0
    e[n] = 0

    if n == 1 # Nothing to do for a 1x1 matrix.
        return
    end
    for l = 1:n
        for j = 1:maxits
            # Look for small off-diagonal elements.
            m = n
            for i = l:n-1
                if abs(e[i]) <= eps(T) * ( abs(d[i]) + abs(d[i+1]) )
                    m = i
                    break
                end
            end
            p = d[l]
            if m == l
                continue
            end
            if j == maxits
                msg = @sprintf("No convergence after %d iterations", j)
                msg *= " (try increasing maxits)"
                error(msg)
            end
            # Form shift
            g = ( d[l+1] - p ) / ( 2 * e[l] )
            r = hypot(g, one(T))
            g = d[m] - p + e[l] / ( g + copysign(r, g) )
            s = one(T)
            c = one(T)
            p = zero(T)
            for i = m-1:-1:l
                f = s * e[i]
                b = c * e[i]
                if abs(f) <  abs(g)
                    s = f / g
                    r = hypot(s, one(T))
                    e[i+1] = g * r
                    c = one(T) / r
                    s *= c
                else
                    c = g / f
                    r = hypot(c, one(T))
                    e[i+1] = f * r
                    s = one(T) / r
                    c *= s
                end
                g = d[i+1] - p
                r = ( d[i] - g ) * s + 2 * c * b
                p = s * r
                d[i+1] = g + p
                g = c * r - b
                # Form first component of vector.
                f = z[i+1]
                z[i+1] = s * z[i] + c * f
                z[i]   = c * z[i] - s * f
            end # loop over i
            d[l] -= p
            e[l] = g
            e[m] = zero(T)
        end # loop over j
    end # loop over l
end

function orthonormal_poly(x::Array{T,1},
                         a::Array{T,1}, b::Array{T,1}, muzero::T) where {T<:AbstractFloat}
    # p[i,j] = value at x[i] of orthonormal polynomial of degree j-1.
    m = length(x)
    n = length(a)
    p = zeros(T, m, n+1)
    c = one(T) / sqrt(muzero)
    rb = one(T) / b[1]
    for i = 1:m
        p[i,1] = c
        p[i,2] = rb * ( x[i] - a[1] ) * c
    end
    for j = 2:n
       rb = one(T) / b[j]
       for i = 1:m
           p[i,j+1] = rb * ( (x[i]-a[j]) * p[i,j]
                                - b[j-1] * p[i,j-1] )
       end
    end
    return p
end

end


text:
import Base.timedwait
import Base.isapprox
import PrettyPrint

using Test
using Scruff
using Scruff.Utils
using Scruff.RTUtils
using Scruff.Models
using Scruff.SFuncs
using Scruff.Operators
import Scruff.Algorithms: LSFI, probability_bounds, prepare, refine, ve

using Scruff.MultiInterface

import Scruff.Operators: Support, bounded_probs, support_quality, support
import Scruff.SFuncs: Dist

struct MyFlip <: Dist{Bool} end

solve_count = 0

@impl begin
    struct MyFlipBoundedProbs end
    function bounded_probs(sf::MyFlip,
            range::__OptVec{Bool},
            parranges::NTuple{N,Vector})::Tuple{Vector{<:AbstractFloat}, Vector{<:AbstractFloat}} where {N}

        global solve_count += 1
        return ([0.1, 0.9], [0.1, 0.9])
    end
end

@impl begin
    struct MyFlipSupport end

    function support(sf::MyFlip,
                    parranges::NTuple{N,Vector},
                    size::Integer,
                    curr::Vector{Bool}) where {N}
        [false, true]
    end
end

@impl begin
    struct MyFlipSupportQuality end
    function support_quality(::MyFlip, parranges)
        :CompleteSupport
    end
end

@testset "lsfi" begin

    @testset "expander" begin

        @testset "Construction" begin
            x1 = Cat([1,2], [0.1, 0.9])()
            net1 = InstantNetwork(Variable[x1(:out)], VariableGraph())

            x2 = Cat([1,2], [0.2, 0.8])()
            net2 = InstantNetwork(Variable[x2(:out)], VariableGraph())

            f(b) = b ? net1 : net2
            b3v = Flip(0.5)()(:b3)
            e3v = Expander(f, Tuple{Bool}, Int)()(:out)
            net3 = InstantNetwork(Variable[e3v,b3v], VariableGraph(e3v=>[b3v]))

            run3 = Runtime(net3)
            ensure_expansion_state!(run3)

            @test isempty(get_state(run3, :subnets))
            @test isempty(get_state(run3, :subruntimes))
        end

        @testset "Expand" begin
            x1 = Cat([1,2], [0.1, 0.9])()
            net1 = InstantNetwork(Variable[x1(:out)], VariableGraph())

            x2 = Cat([1,2], [0.2, 0.8])()
            net2 = InstantNetwork(Variable[x2(:out)], VariableGraph())

            f(b) = b ? net1 : net2
            vb3 = Flip(0.5)()(:b3)
            ve3 = Expander(f, Tuple{Bool}, Int)()(:out)
            net3 = InstantNetwork(Variable[ve3,vb3], VariableGraph(ve3=>[vb3]))

            run3 = Runtime(net3)
            @testset "The first time" begin
                expand!(run3, ve3, true)
                @test expanded(run3, ve3, true)
                @test expansion(run3, ve3, true) == net1
                run1 = subruntime(run3, ve3, net1)
                @test run1.network == net1
                @test !expanded(run3, ve3, false)
            end
            @testset "The second time" begin
                expand!(run3, ve3, true)
                @test expanded(run3, ve3, true)
                @test !expanded(run3, ve3, false)
            end
        end

        @testset "Getting range" begin
            x1 = Cat([1,2], [0.1, 0.9])()
            net1 = InstantNetwork(Variable[x1(:out)], VariableGraph())

            x2 = Cat([1,2], [0.2, 0.8])()
            net2 = InstantNetwork(Variable[x2(:out)], VariableGraph())

            f₀(b) = (b) = b ? net1 : net2
            vb3 = Flip(0.5)()(:b3)
            ve3 = Expander(f₀ , Tuple{Bool}, Int)()(:out)
            net3 = InstantNetwork(Variable[ve3,vb3], VariableGraph(ve3=>[vb3]))

            @testset "With size requiring all expansions" begin
                run3 = Runtime(net3)
                ensure_all!(run3)
                instb3 = current_instance(run3, vb3)
                set_range!(run3, instb3, [false, true], 2)
                (rng, _) = expander_range(run3, ve3, 4, 2)
                @test expanded(run3, ve3, false)
                @test expanded(run3, ve3, true)
                @test rng == [1, 2]
            end

            @testset "With size not requiring all expansions" begin
                run3 = Runtime(net3)
                ensure_all!(run3)
                instb3 = current_instance(run3, vb3)
                set_range!(run3, instb3, [false, true], 2)
                (rng, _) = expander_range(run3, ve3, 1, 2)
                @test expanded(run3, ve3, false)
                @test !expanded(run3, ve3, true)
                @test rng == [1, 2]
            end

            @testset "With recursion" begin
                b1 = Cat([1,2], [0.1, 0.9])()
                net1 = InstantNetwork(Variable[b1(:out)], VariableGraph())

                b2 = Cat([1,2], [0.2, 0.8])()
                net2 = InstantNetwork(Variable[b2(:out)], VariableGraph())

                f₁(b) = b ? net1 : net2
                vb3 = Flip(0.5)()(:b3)
                ve3 = Expander(f₁, Tuple{Bool}, Int)()(:out)
                net3 = InstantNetwork(Variable[ve3,vb3], VariableGraph(ve3=>[vb3]))

                f₂(b) = b ? net1 : net2
                vb3 = Flip(0.7)()(:b3)
                ve3 = Expander(f₂, Tuple{Bool}, Int)()(:out)
                net3 = InstantNetwork(Variable[vb3,ve3], VariableGraph(ve3=>[vb3]))

                g(b) = net3
                vb4 = Flip(0.6)()(:b4)
                ve4 = Expander(g, Tuple{Bool}, Int)()(:out)
                net4 = InstantNetwork(Variable[vb4,ve4], VariableGraph(ve4=>[vb4]))
                run4 = Runtime(net4)

                ensure_all!(run4)
                expand!(run4, ve4, false)
                expand!(run4, ve4, true)

                run3 = subruntime(run4, ve4, net3)
                expand!(run3, ve3, false)
                expand!(run3, ve3, true)

                instb4 = current_instance(run4, vb4)
                set_range!(run4, instb4, [false, true], 2)

                (rng, _) = expander_range(run4, ve4, 4, 2)
                @test rng == [] # recursion depth too shallow

                (rng, _) = expander_range(run4, ve4, 4, 3)
                @test Set(rng) == Set([1, 2])
            end

        end

        @testset "Expanding range" begin
            vx1 = Normal(-0.1, 1.0)()(:out)
            net1 = InstantNetwork(Variable[vx1], VariableGraph())

            vx2 = Normal(0.4, 1.0)()(:x2)
            vy2 = LinearGaussian((0.7,), 0.0, 1.0)()(:out)
            net2 = InstantNetwork(Variable[vx2,vy2], VariableGraph(vy2=>[vx2]))

            g(b) = b ? net1 : net2
            vb3 = Flip(0.5)()(:b3)
            ve3 = Expander(g, Tuple{Bool}, Float64)()(:out)
            net3 = InstantNetwork(Variable[vb3,ve3], VariableGraph(ve3=>[vb3]))

            run3 = Runtime(net3)

            instb3 = instantiate!(run3, vb3, 0)
            set_range!(run3, instb3, [false, true], 2)
            (r1,_) = expander_range(run3, ve3, 5, 2)
            (r2,_) = expander_range(run3, ve3, 10, 2)

            @test length(r2) >= 10
            @test issorted(r2)
            @test length(Set(r2)) == length(r2)
            @test issubset(Set(r1), Set(r2))
        end

        @testset "Computing probability bounds" begin
            vx1 = Cat([1,2], [0.1, 0.9])()(:out)
            net1 = InstantNetwork(Variable[vx1], VariableGraph())

            vx2 = Cat([1,2], [0.2, 0.8])()(:out)
            net2 = InstantNetwork(Variable[vx2], VariableGraph())

            f(b) = b ? net1 : net2
            vb3 = Flip(0.5)()(:b3)
            ve3 = Expander(f, Tuple{Bool}, Int)()(:out)
            net3 = InstantNetwork(Variable[vb3,ve3], VariableGraph(ve3=>[vb3]))

            @testset "With a fully expanded model" begin
                run3 = Runtime(net3)
                ensure_all!(run3)
                instb3 = current_instance(run3, vb3)
                set_range!(run3, instb3, [false, true])
                ord = topsort(get_initial_graph(net3))
                set_ranges!(run3, Dict{Symbol, Score}(), 4, 2, ord)
                f₃(runtime, order, query_vars, depth) = ve(runtime, order, query_vars; depth = depth, bounds = true)
                (l, u) = expander_probs(run3, f₃, ve3, 2)

                @test l == [0.2, 0.8, 0.1, 0.9]
                @test u == [0.2, 0.8, 0.1, 0.9]
            end

            @testset "With a partially expanded model" begin
                run3 = Runtime(net3)
                ensure_all!(run3)
                instb3 = current_instance(run3, vb3)
                set_range!(run3, instb3, [false, true])
                ord = topsort(get_initial_graph(net3))
                set_ranges!(run3, Dict{Symbol, Score}(), 1, 2, ord)
                f₄(runtime, order, query_vars, depth) = ve(runtime, order, query_vars; depth = depth, bounds = true)
                (l, u) = expander_probs(run3, f₄, ve3, 2)

                @test l == [0.2, 0.8, 0, 0]
                @test u == [0.2, 0.8, 1, 1]
            end
        end

    end

    @testset "inside lsfi" begin

        @testset "Completeness" begin
            vx1 = Cat([1,2], [0.1, 0.9])()(:out)
            net1 = InstantNetwork(Variable[vx1], VariableGraph())

            vx2 = Cat([1,2], [0.2, 0.8])()(:out)
            net2 = InstantNetwork(Variable[vx2], VariableGraph())

            vx3 = Normal(-0.1, 1.0)()(:out)
            net3 = InstantNetwork(Variable[vx3], VariableGraph())

            vx4 = Normal(0.4, 1.0)()(:x4)
            vy4 = LinearGaussian((0.7,), 0.0, 1.0)()(:out)
            net4 = InstantNetwork(Variable[vx4,vy4], VariableGraph(vy4=>[vx4]))

            f5(b) = b ? net1 : net2
            vb5 = Flip(0.5)()(:b5)
            ve5 = Expander(f5, Tuple{Bool}, Int)()(:out)
            net5 = InstantNetwork(Variable[vb5,ve5], VariableGraph(ve5=>[vb5]))

            f6(b) = b ? net3 : net4
            vb6 = Flip(0.6)()(:b6)
            ve6 = Expander(f6, Tuple{Bool}, Float64)()(:out)
            net6 = InstantNetwork(Variable[vb6,ve6], VariableGraph(ve6=>[vb6]))

            run5 = Runtime(net5)
            ord5 = topsort(get_initial_graph(net5))
            run6 = Runtime(net6)
            ord6 = topsort(get_initial_graph(net6))
            ensure_all!(run5, 0)
            ensure_all!(run6, 0)
            expand!(run5, ve5, false)
            expand!(run5, ve5, true)
            expand!(run6, ve6, false)
            expand!(run6, ve6, true)

            @test !RTUtils._complete(run5)

            set_ranges!(run5, Dict{Symbol, Score}(), 1, 3, ord5) # range is too small to be complete
            @test !RTUtils._complete(run5)              # because subnetwork not expanded

            set_ranges!(run5, Dict{Symbol, Score}(), 4, 3, ord5) # range is large enough
            @test RTUtils._complete(run5)

            set_ranges!(run6, Dict{Symbol, Score}(), 10, 3, ord6)
            @test !RTUtils._complete(run6)

            set_ranges!(run6, Dict{Symbol, Score}(), 20, 3, ord6) # continuous range not complete
            @test !RTUtils._complete(run6)
        end

        @testset "Completeness with recursion" begin
            f(b) = net1
            vb1 = Flip(0.1)()(:b1)
            ve1 = Expander(f, Tuple{Bool}, Bool)()(:out)
            net1 = InstantNetwork(Variable[vb1,ve1], VariableGraph(ve1=>[vb1]))

            run1 = Runtime(net1)
            expand!(run1, ve1, false)
            expand!(run1, ve1, true)

            @test !RTUtils._complete(run1)
        end

        @testset "Run to completion" begin
            v1 = Cat([1,2], [0.1, 0.9])()(:out)
            net1 = InstantNetwork(Variable[v1], VariableGraph())

            x2 = Cat([1,2], [0.2, 0.8])()(:out)
            net2 = InstantNetwork(Variable[x2], VariableGraph())

            vb5 = Flip(0.4)()(:b5)
            f5(b) = b ? net2 : net1
            ve5 = Expander(f5, Tuple{Bool}, Int)()(:out)
            net5 = InstantNetwork(Variable[vb5,ve5], VariableGraph(ve5=>[vb5]))

            run5 = Runtime(net5)
            ord5 = topsort(get_initial_graph(net5))
            ensure_all!(run5, 0)
            alg = LSFI([ve5]; start_depth = 2)
            prepare(alg, run5)
            refine(alg, run5)
            inst = current_instance(run5, ve5)
            range = get_range(run5, inst)
            (ls,us) = probability_bounds(alg, run5, inst, range)
            pa = 0.4 * 0.2 + 0.6 * 0.1
            pb = 0.4 * 0.8 + 0.6 * 0.9
            @test isapprox(ls, [pa, pb])
            @test isapprox(us, [pa, pb])
        end

        @testset "With dynamic programming" begin
            expand_count = 0

            global function PrettyPrint.pp_impl(io, frame::StackTraces.StackFrame, indent::Int)
                line = frame.inlined ? "[inlined]" : "$(frame.line)"
                print(io, "$(frame.func) at $(frame.file):$(line)")
                indent
            end

            b1 = MyFlip()()(:out)
            net1 = InstantNetwork(Variable[b1], VariableGraph())

            function f1(x)
                if x expand_count += 1 end
                net1
            end

            vb2 = Flip(0.2)()(:b2)
            ve2 = Expander(f1, Tuple{Bool}, Bool)()(:out)
            net2 = InstantNetwork(Variable[vb2,ve2], VariableGraph(ve2=>[vb2]))

            vb3 = Flip(0.3)()(:b3)
            ve3 = Expander(x -> net2, Tuple{Bool}, Bool)()(:out)
            net3 = InstantNetwork(Variable[vb3,ve3], VariableGraph(ve3=>[vb3]))

            vb4 = Flip(0.4)()(:b4)
            ve4 = Expander(x -> net2, Tuple{Bool}, Bool)()(:out)
            net4 = InstantNetwork(Variable[vb4,ve4], VariableGraph(ve4=>[vb4]))

            vb5 = Flip(0.5)()(:b5)
            ve5 = Expander(b -> b ? net3 : net4, Tuple{Bool}, Bool)()(:out)
            net5 = InstantNetwork(Variable[vb5,ve5], VariableGraph(ve5=>[vb5]))

            run5 = Runtime(net5)
            ensure_all!(run5, 0)
            order = topsort(get_initial_graph(net5))
            alg = LSFI([ve5]; start_depth = 10)
            prepare(alg, run5)
            refine(alg, run5)
            inst = current_instance(run5, ve5)
            range = get_range(run5, inst)
            (lower, upper) = probability_bounds(alg, run5, inst, range)
            @test isapprox(lower, [0.1, 0.9], atol = 0.0000001)
            @test isapprox(upper, lower)
        end

        @testset "With recursion" begin
            bx1 = Flip(1.0)()(:out)
            netx1 = InstantNetwork(Variable[bx1], VariableGraph())

            vx2 = Flip(0.9)()(:b2)
            fx(b) = b ? netx2 : netx1
            vex2 = Expander(fx, Tuple{Bool}, Bool)()(:out)
            netx2 = InstantNetwork(Variable[vx2,vex2], VariableGraph(vex2=>[vx2]))

            runtime = Runtime(netx2)
            ensure_all!(runtime, 0)
            alg = LSFI([vex2]; start_depth = 11)
            prepare(alg, runtime)
            refine(alg, runtime)
            @test true # This test succeeds if the above call terminates
        end

    end

    @testset "Using the LazyInference interface" begin
        @testset "LSFI" begin
            m1 = Cat([1,2], [0.1, 0.9])()
            v1 = m1(:out)
            net1 = InstantNetwork(Variable[v1], VariableGraph())

            m2 = Cat([1,2], [0.2, 0.8])()
            v2 = m2(:out)
            net2 = InstantNetwork(Variable[v2], VariableGraph())

            f(b) = b ? net1 : net2
            m3 = Flip(0.9)()
            m4 = Expander(f, Tuple{Bool}, Int)()
            m5 = DiscreteCPT([1,2], Dict((1,) => [0.3, 0.7], (2,) => [0.4, 0.6]))()
            v3 = m3(:v3)
            v4 = m4(:v4)
            v5 = m5(:out)
            net3 = InstantNetwork(Variable[v3,v4,v5], VariableGraph(v4 => [v3], v5 => [v4]))

            alg = LSFI([get_node(net3, :out)]; increment = 3, max_iterations = 20, start_size=3)
                runtime = Runtime(net3)
            ensure_all!(runtime, 0)
            prepare(alg, runtime)
            refine(alg, runtime)
            @test has_instance(runtime, v5)
            i5 = current_instance(runtime, v5)
            state = alg.state
            @test state.next_size == 6
            @test state.next_iteration == 2
            @test state.next_depth == 2
            @test !state.is_complete
            # At first expansion, the bounds for i5 are (0,1)
            (ls1, us1) = probability_bounds(alg, runtime, i5, [1,2])
            @test isapprox(ls1[1], 0.0)
            @test isapprox(us1[1], 1.0)

            q2 = 0.9 * (0.1 * 0.3 + 0.9 * 0.4) + 0.1 * (0.2 * 0.3 + 0.8 * 0.4)
            refine(alg, runtime)
            @test has_instance(runtime, v5)
            i5 = current_instance(runtime, v5)
            state = alg.state
            @test state.next_size == 9
            @test state.next_iteration == 3
            @test state.next_depth == 3
            @test state.is_complete
            # At second expansion, lower bounds = upper bounds = q2
            (ls2, us2) = probability_bounds(alg, runtime, i5, [1,2])
            @test isapprox(ls2[1], q2)
            @test isapprox(us2[1], q2)
        end
    end

end


text:
# Youth_Suicide_Deaths_in_Washington_State_by_Gender_Age_0-17_Years__from_2008-2012.csv
# State of Washington — Youth Suicide Deaths in Washington State by Gender Age 0-17 Years, from 2008-2012
# https://catalog.data.gov/dataset/youth-suicide-deaths-in-washington-state-by-gender-age-0-17-years-from-2008-2012
#
# CLEAN
# 1. Too much unnecessary data
#
# TIDY
# 1. Values are in the columns
# 2. Multiple values in the same column (after stack)
#
# VISUALIZE
# 1. Bar chart - King county male youth suicides
# 2. Grouped Bar chart - King county youth suicides
# 3. Line chart - five year suicide trend by gender
# 4. Heatmap - county/year suicides

using DataFrames, CSV, Pipe, TabularDisplay

df = DataFrame(CSV.File("data/youth_suicide.csv"))

# Examine column names
names(df)
displaytable(stdout, names(df); index = true)
df

# Both Male/Female and Year are encoded in the columns
# Normally, time series are represented in rows as they're separate observations.
# Let's convert to long form first.
sdf = stack(df, 2:19; variable_name = :type_year)

# Now, spread it into two separate columns
sdf.type = [split(x, " ")[1] for x in sdf.type_year]
sdf.year = [split(x, " ")[2] for x in sdf.type_year]
deletecols!(sdf, :type_year)

# Let's check if the Total type really contains the total
# Pick a county.
sdf[sdf.County .== "Yakima", :]

# Remove the multi-year data since they're derivable
sdf = @pipe sdf |>
    filter(:type => !=("Total"), _) |>
    filter(:year => !=("(2008-2012)"), _)

# Check Yakima again
sdf[sdf.County .== "Yakima", :]

# Check unique values of year again
unique(sdf.year)

# Convert the year column to integer
sdf.year = [parse(Int, replace(x, r"[()]" => "")) for x in sdf.year]
sdf[sdf.County .== "Yakima", :]

# Examine
describe(sdf, :eltype, :min, :max, :nunique, :nmissing)
sdf

# Male and Female variables can be put back as columns now
df = unstack(sdf, :type, :value)
df[df.County .== "Yakima", :]

# Now, it's all tidy!

using Plots, StatsPlots

# Plot King county's yearly rate

filter(:County => ==("King"), df)

# line chart
@pipe df |>
    filter(:County => ==("King"), _) |>
    plot(_.year, _.Male,
        title = "King County Youth Suicides (Male)",
        legend = :none,
        size = (300, 300))

# bar chart
@pipe df |>
    filter(:County => ==("King"), _) |>
    bar(_.year, _.Male,
        title = "King County Youth Suicides (Male)",
        legend = :none,
        size = (300, 300))

# Prepare to plot both Male and Female together
df_stacked = stack(df, 3:4; variable_name = "gender")

filter(:County => ==("King"), df_stacked)

@pipe df_stacked |>
    filter(:County => ==("King"), _) |>
    groupedbar(    # StatsPlots.jl
        _.year,    # x-axis
        _.value;   # y-axis
        group = _.gender,
        bar_position = :stack,
        bar_width = 0.7,
        title = "King County Suicides",
        size = (300, 300),
        legend = :topleft)

# Let's try grouping function
df
@pipe df |>
    groupby(_, :County)

# 5 year totals
@pipe df |>
    groupby(_, :County) |>
    combine(_, :Female => sum, :Male => sum)

# 5 year totals with male/female combined
@pipe df |>
    groupby(_, :County) |>
    combine(_, :Female => sum, :Male => sum) |>
    select(_, :County, [:Female_sum, :Male_sum] => ByRow(+) => :Total)

# Just for fun, can we transpose the original data frame?
# @pipe df |>
#     groupby(_, :year) |>
#     combine(_, :Female => sum => :Female, :Male => sum => :Male) |>
#     stack(_, 2:3; variable_name = :gender) |>
#     unstack(_, :year, :value)

# Let's plot more charts

let
    data = @pipe df |>
        groupby(_, :year) |>
        combine(_, :Female => sum => :Female, :Male => sum => :Male) |>
        stack(_, 2:3; variable_name = :gender, value_name = :suicides)
    plot(data.year, data.suicides; groups = data.gender,
        title = "Youth Suicides Trend",
        legend = :topleft,
        linewidth = 3,
        size = (300, 300))
end

# Heatmap to show suicide stats by (year, county)
let
    data = @pipe df |>
        select(_, :County, :year, [:Female, :Male] => ByRow(+) => :total) |>
        unstack(_, :year, :total)
    values = Matrix(data[:, 2:end])
    StatsPlots.heatmap(2008:2012, data.County, values;
        title = "Youth Suicide Heatmap",
        xticks = :all,
        yticks = :all,
        size = (400, 600))
end


text:
include("../src/boardrobots.jl")
include("../src/solver.jl")
include("../src/readchallenges.jl")

function main()
    board = read_challenge((@__DIR__)*"/../instances/"*"challenge39.txt", autodetectformat=true)
    solutions = find_solution(board)
    print_solutions(board, solutions)
    return 0
end

@time main()

max_stars_repo_name: remi-garcia/LunarLockout

max_issues_repo_name: remi-garcia/LunarLockout

max_forks_repo_name: remi-garcia/LunarLockout


text:
using Plots

include(string(pwd(),"/src/simulation.jl"))
include(string(pwd(),"/gali_bayesian.jl"))
include(string(pwd(),"/src/likelihood.jl"))

dados,choques = simulate_dsge(GAMMA_0,GAMMA_1,PSI,PI,5000)

plot(dados[:,1], label = "Quartely CPI")
plot!(dados[:,2], label = "GDP")
plot!(dados[:,3], label = "i")

plot(dados[:,4], label = "v")
plot!(dados[:,3], label = "i")

param = [alfa,bet,epsilon,theta,sig,1,phi,phi_pi,phi_y,rho_v]

ss, filtered,loglike = log_like_dsge(param,dados[:,2])

plot(dados[1:100,1])
plot!(filtered[1:100,1])

param1 = [alfa*2,bet,epsilon,theta,sig,1,phi,phi_pi,phi_y,rho_v]

a2,ll = log_like_dsge(param1,dados[:,2])

bett = 0.70:0.005:1.15

vals_bet = collect(bett)

for i in 1:length(bett)
    param1 = [alfa,bett[i],epsilon,theta,sig,1,phi,phi_pi,phi_y,rho_v]
    ss,filtered,vals_bet[i] = log_like_dsge(param1,dados[:,2])
    println(i)
end

plot(bett,vals_bet)

alff = 0.10:0.01:0.9

vals_alfa = collect(alff)

for i in 1:length(alff)
    param1 = [alff[i],bet,epsilon,theta,sig,1,phi,phi_pi,phi_y,rho_v]
    ss,filtered,vals_alfa[i] = log_like_dsge(param1,dados[:,2])
end

plot(alff,vals)

bett = 0.70:0.005:1.15
alff = 0.10:0.01:0.9

vals = zeros(length(alff),length(bett))

for i in 1:length(alff), j in 1:length(bett)
    param1 = [alff[i],bett[j],epsilon,theta,sig,1,phi,phi_pi,phi_y,rho_v]
    ss,filtered,vals[i,j] = log_like_dsge(param1,dados[:,2])
end

plotly()

plot(alff,bett,vals, st=:surface, camera = (-45,45))
xaxis!("alfa")
yaxis!("beta")

using JLD2

@save "surface.jld2" alff bett vals

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
mutable struct Topography

    topoz_sT :: Array{Float64, 3}
    Nz_bot_sT :: Array{Int64, 3}

    mask_T      :: AbstractArray{Float64, 3}
    sfcmask_sT  :: AbstractArray{Float64, 3}
    deepmask_T  :: AbstractArray{Float64, 3}

    function Topography(
        Nz_bot        :: Union{AbstractArray{Int64,2}, Nothing},
        Nx            :: Int64,
        Ny            :: Int64,
        z_w           :: AbstractArray{Float64, 1};
        deep_depth    :: Float64 = 500.0
    )

        Nz = length(z_w) - 1

        sT_pts = Nx * Ny

        topoz_sT = zeros(Float64, 1, Nx, Ny)

        if Nz_bot == nothing
            Nz_bot_sT = zeros(Int64, 1, Nx, Ny)
            Nz_bot_sT .= length(z_w) - 1
        else
            if size(Nz_bot) != (Nx, Ny)
                throw(ErrorException(format("Size of `Nz_bot` variable is not the same as given (Nx, Ny) = ({:d}, {:d})", Nx, Ny)))
            end

            Nz_bot_sT = copy(reshape(Nz_bot, 1, Nx, Ny))

        end

        for j=1:Ny, i=1:Nx
            if Nz_bot_sT[1, i, j] > Nz
                Nz_bot_sT[1, i, j] = Nz
            end
            topoz_sT[1, i, j] = z_w[Nz_bot_sT[1, i, j]+1]
        end

        deep_Nz = lookupNz(- deep_depth, z_w)
        println(format("Water depth of {:f} m is considered as deep. The found deep_Nz = {:d}", deep_depth, deep_Nz))
        _lookupNz = (z) -> lookupNz(z, z_w)
        @. Nz_bot_sT = _lookupNz(topoz_sT)

        mask_T = ones(Float64, Nz, Nx, Ny)

        # turn the grids below topography to 0
        for j=1:Ny, i=1:Nx

            if Nz_bot_sT[1, i, j] < Nz
                mask_T[Nz_bot_sT[1, i, j]+1:end, i, j] .= 0.0
            end
        end

        sfcmask_sT = view(mask_T, 1:1, :, :)

        deepmask_T = copy(mask_T)
        # turn the entire column to 0 if it is too shallow
        for j=1:Ny, i=1:Nx
            if Nz_bot_sT[1, i, j] < deep_Nz
                deepmask_T[:, i, j] .= 0.0
            end
        end

        return new(
            topoz_sT,
            Nz_bot_sT,
            mask_T,
            sfcmask_sT,
            deepmask_T,
        )

    end

end

function lookupNz(z, z_w)

    max_Nz = length(z_w) - 1

    if z == 0.0
        return 0
    end

    for i=1:max_Nz-1
        if z_w[i] > z >= z_w[i+1]
            return i
        end
    end

    return max_Nz

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
T=PrepTest()

function KMeansModel(A::Array{Float64,2},NumQueries::Int,Clusters::Int)
    DQ=NewDataMatrix()
    DQ.A=A
    DQ.NumQueries=NumQueries
    (ma,na)=size(DQ.A)
    D=DQ.A[:,DQ.NumQueries+1:na]
    Q=DQ.A[:,1:NumQueries]
    KM=kmeans(D,Clusters)
    C=KM.centers
    (P,R)=qr(C)
    DS = P'*D #Using the P we find the projections of the documents
    (md,nd)=size(D)
    (mq,nq)=size(Q)
    Costheta_KM=zeros(nq,nd)
    #Set a tol level
    tol=linspace(0.01,0.99,20)
    lt=length(tol)
    KM_TDM_Result=TMResults()
    KM_TDM_Result=InitResults!(KM_TDM_Result,nq,lt)
    qindex=zeros(nq)
    for z=1:1:nq
        qindex[z]=z
        q=P'*Q[:,z]
        for i=1:nd
            Costheta_KM[z:z,i:i]=q'*DS[:,i]/(norm(q,2)*norm(DS[:,i],2))
        end
        for k=1:1:lt
            (KM_TDM_Result.Dr[z,k],
             KM_TDM_Result.Dt[z,k],
             KM_TDM_Result.Nr[z,k],
             KM_TDM_Result.Rec[z,k],
             KM_TDM_Result.Prec[z,k])=
             FindResults(k,z,Costheta_KM[z,:],tol[k])
        end
    end
    (KM_TDM_Result.RecFinal,KM_TDM_Result.PrecFinal)=average_RecPrec(KM_TDM_Result.Rec,KM_TDM_Result.Prec,qindex)
    return KM_TDM_Result.RecFinal,KM_TDM_Result.PrecFinal
end

function KMeansModel(QueryNum::Int,A::Array{Float64,2},NumQueries::Int,Clusters::Int)
    DQ=NewDataMatrix()
    DQ.A=A
    DQ.NumQueries=NumQueries
    (ma,na)=size(DQ.A)
    D=DQ.A[:,DQ.NumQueries+1:1063]
    Q=DQ.A[:,QueryNum]
    KM=kmeans(D,50)
    C=KM.centers
    (P,R)=qr(C)
    DS = P'*D #Using the P we find the projections of the documents
    (md,nd)=size(D)
    (mq,nq)=size(Q)
    Costheta_KM=zeros(nq,nd)
    #Set a tol level
    tol=linspace(0.01,0.99,20)
    lt=length(tol)
    KM_1Query_Result=TMResults()
    KM_1Query_Result=InitResults!(KM_TDM_Result,nq,lt)
    qindex=zeros(nq)
    for z=1:1:nq
        qindex[z]=z
        q=P'*Q[:,z]
        for i=1:nd
            Costheta_KM[z:z,i:i]=q'*DS[:,i]/(norm(q,2)*norm(DS[:,i],2))
        end
        for k=1:1:lt
            (KM_1Query_Result.Dr[z,k],
             KM_1Query_Result.Dt[z,k],
             KM_1Query_Result.Nr[z,k],
             KM_1Query_Result.Rec[z,k],
             KM_1Query_Result.Prec[z,k])=
             FindResults(k,z,Costheta_KM[z,:],tol[k])
        end
    end
    (KM_1Query_Result.RecFinal,KM_1Query_Result.PrecFinal)=average_RecPrec(KM_1Query_Result.Rec,KM_1Query_Result.Prec,qindex)
    return KM_1Query_Result.RecFinal,KM_1Query_Result.PrecFinal
end

function KMeansModel(Q_C::Corpus,D_C::Corpus,Clusters::Int)
    #Q_c and D_c are the query and Document corpuses.
    PreProcess!(Q_C)
    PreProcess!(D_C)
    #Compulsory update the lexicon, to identify the keywords
    update_lexicon!(Q_C)
    update_lexicon!(D_C)
    #Get the lexicon of the documents.
    DocLex=lexicon(D_C)
    #Cant use the above as query vectors since, the queries must also be wrt
    #Documents lexicon, and hence the dimensions will match
    D=full(tdm(D_C))
    (md,nd)=size(D)
    nq=length(Q_C.documents)
    Q=zeros(md,nq)
    for j=1:nq
        Q[:,j]=dtv(Q_C.documents[j],DocLex)'
    end
    QD=[Q D]
    r=KmeansModel(QD,nq,Clusters)
    return r
end


text:
################################################################################
#
# Document
# http://abelsiqueira.github.io/blog/test-driven-development-in-julia/
#
################################################################################
# The package directory needs to be "instantiate" in order to run from the
# package directory.
#
# Run test from command line
#   $ cd OilData/
#   $ julia --project=@. test/runtests.jl
#
# Run test from REPL
#   $ cd OilData/
#   $ julia --project=@.
#   julia> include("test/runtests.jl")
#
# The Pkg test will call "instantiate" and creates the "Manifest.toml" file.
#
# Run test from REPL using Pkg
#   $ cd OilData/
#   $ julia --project=@.
#   julia> import Pkg
#   julia> Pkg.test("OilData")
#
# Run test from Package mode
#   $ cd OilData/
#   $ julia --project=@.
#   julia> ]
#   (OilData) pkg> test OilData
#
################################################################################
# Add package to test/.
#   cd OilData/
#   julia --project=test
#   ]
#   add DataFrames
################################################################################
# Run from REPL, in development environment
#   $ cd OilDataDev/
#   $ julia --project=@.
#   julia> include("dev/OilData/test/runtests.jl")
#   julia> include("dev/OilData/test/test_parser.jl")
#   julia> include("dev/OilData/test/test_utility.jl")
# One can run in the development environment
#   julia> using Revise
# to speed up the development cycle.
################################################################################

# In a development parent project such as
#   OilDataDev/
# where OilData.jl is checkout in dev/
#   OilDataDev/dev/OilData/
#
#   cd OilDataDev/
#   julia --project=@.
#   using Revise
#   include("dev/OilData/test/runtests.jl")
# Test individual files
#   include("dev/OilData/test/test_parser.jl")
#   include("dev/OilData/test/test_utility.jl")

# Not working
# Is this the best way for development?
#   cd OilData/
#   julia --project=test
#

using Test
using Dates

# see https://discourse.julialang.org/t/writing-tests-in-vs-code-workflow-autocomplete-and-tooltips/57488
# see https://github.com/julia-vscode/julia-vscode/issues/800
if isdefined(@__MODULE__, :LanguageServer)
    # invoked by VS Code
    include("../src/OilData.jl")
    using .OilData
else
    # invoked during test
    using OilData
end

include("common.jl")

# This lets test_*.jl to detect if it is run
# as includes or standalone.
local_6cfddc95df03456583406d57fd5963ac = false

include("test_parser.jl")
include("test_utility.jl")


text:
"""
    Predicate{N}

The abstract supertype of predicates.
"""
abstract type Predicate{N} end

"""
    evaluate(p::Predicate, args...)

Evaluate a predicate for given arguments.

### Input

- `p`    -- predicate
- `args` -- argument list

### Output

A `Bool` corresponding to the evaluation of `p(args)`.

### Notes

We offer `p(args)` as a short-hand version of `evaluate(p, args)`.
"""
function evaluate end

function assert_same_length(::Predicate{Val{n}}, args...) where {n}
    @assert length(args) == n "incompatible number of arguments " *
        "($(length(args)) vs. $n)"
end


text:
"""
    reproject(source::GeoFormat, target::GeoFormat, dim::Dimension, val)

`reproject` uses ArchGDAL.reproject, but implemented for a reprojecting
a value array of values, a single dimension at a time.
"""
reproject(source, target, dim::Dimension, val) = val
function reproject(source::GeoFormat, target::GeoFormat, dim::XDim, val::Number)
    AG.reproject((val, zero(val)), source, target; order=:trad)[1]
end
function reproject(source::GeoFormat, target::GeoFormat, dim::YDim, val::Number)
    AG.reproject((zero(val), val), source, target; order=:trad)[2]
end
function reproject(source::GeoFormat, target::GeoFormat, ::XDim, vals::AbstractArray)
    rep = AG.reproject(map(v -> (v, zero(v)), vals), source, target; order=:trad)
    map(r -> r[1], rep)
end
function reproject(source::GeoFormat, target::GeoFormat, dim::YDim, vals::AbstractArray)
    rep = AG.reproject(map(v -> (zero(v), v), vals), source, target; order=:trad)
    map(r -> r[2], rep)
end
function reproject(source::GeoFormat, target::GeoFormat, ::XDim, vals::Tuple)
    reps = AG.reproject([(v, zero(v)) for v in vals], source, target; order=:trad)
    Tuple(r[1] for r in reps)
end
function reproject(source::GeoFormat, target::GeoFormat, dim::YDim, vals::Tuple)
    reps = AG.reproject([(zero(v), v) for v in vals], source, target; order=:trad)
    Tuple(r[2] for r in reps)
end


text:
using CryptoGroups
using Test

### Testing prime groups
using Primes

function testgroup(G)
    @test G*G^2 == G^3
    @test (G^2)^2 == G^4
    @test G^(order(G) + 1) == G
end

G = PrimeGroup(5,23,totient(23),4)
@show G^2
@show v = value(G)
@show security(G)
@show order(G)

@test PrimeGroup(binary(G^2),G)==G^2
@test typeof(G)(binary(G^2),G)==G^2

testgroup(G)

### Testing EllitpicGroup

G = CryptoGroups.Scep256k1Group()
@show G^2
@show value(G)
@show security(G)
@show order(G)

@test EllipticGroup(binary(G^2),G)==G^2
@test typeof(G)(binary(G^2),G)==G^2

testgroup(G)

###  RFC standart

CryptoGroups.FirstOakleyGroup()
CryptoGroups.SecondOakleyGroup()

G = CryptoGroups.MODP160Group()
testgroup(G)

G = CryptoGroups.MODP224Group()
testgroup(G)

G = CryptoGroups.MODP256Group()
#testgroup(G)

### Group generation algorithms

using Paillier
rng = Paillier.default_rng()

G = CryptoGroups.SophieGermainGroup(rng,2,100)
testgroup(G)

### Seems to be problem with large numbers
G = CryptoGroups.DSAStandartGroup(rng,10,10)
testgroup(G)


text:
using MathLink
using Test

import MathLink: WExpr, WSymbol

w = W"Factorial"(30)
@test_throws MathLink.MathLinkError weval(Int, w)
@test weval(BigInt, w) == factorial(big(30))

@test weval(W"Factorial"(20)) === factorial(20)

w = W"Sqrt"(2.0)
@test weval(w) == sqrt(2.0)

@test weval(W"Function"(W"x",W"Times"(W"x", 2))(100)) == 200

@test weval(W"Integrate"(W"Log"(W"x"), (W"x", 1, W"E"))) == 1


text:
module KLIEPInference

using Distributions, LinearAlgebra, Statistics, SparseArrays, Random
using StatsBase
using ProximalBase, CoordinateDescent

export
    IsingSampler,
    chain,
    kary_tree,
    removeEdges!,
    addEdges!,

    # KLIEP solvers
    KLIEPSolver,
    CD_KLIEP,
    KLIEP, KLIEP!,
    spKLIEP, spKLIEP!,
    Hinv_row,

    # inference
    BootstrapEstimates,
    boot_SparKLIE,
    boot_quantile,
    boot_quantile_studentized,

    # other
    Ψising,
    unpack,
    pack,
    KLIEP_Hessian,
    stderr_SparKLIE

include("utils.jl")
include("bootstrap.jl")
include("sampler.jl")
include("solver.jl")
include("debias.jl")
include("stderr.jl")

end # module


text:
################################################################################
#
#  Genus symbole
#
################################################################################

function _genus_symbol_kirschmer(L::HermLat, p; uniformizer = zero(order(p)))
  @assert order(p) == base_ring(base_ring(L))

  B, G, S = jordan_decomposition(L, p)
  R = base_ring(L)
  E = nf(R)
  K = base_field(E)
  if !isdyadic(p) || !isramified(R, p)
    sym = Tuple{Int, Int, Bool}[ (nrows(B[i]), S[i], islocal_norm(E, coeff(det(G[i]), 0), p)) for i in 1:length(B)]
  else
    P = prime_decomposition(R, p)[1][1]
    pi = E(K(Hecke.uniformizer(p)))
    sym = Tuple{Int, Int, Bool, Int, elem_type(K)}[]
    for i in 1:length(B)
      normal = (_get_norm_valuation_from_gram_matrix(G[i], P)::Int == S[i])::Bool
      GG = diagonal_matrix(dense_matrix_type(E)[pi^(max(0, S[i] - S[j])) * G[j] for j in 1:length(B)])::dense_matrix_type(E)
      v = _get_norm_valuation_from_gram_matrix(GG, P)::Int
      #_n = norm(lattice(hermitian_space(E, GG), identity_matrix(E, nrows(GG))))
      #vv = valuation(R * norm(lattice(hermitian_space(E, GG), identity_matrix(E, nrows(GG)))), P)::Int
      s = (nrows(B[i]), S[i], normal, v, coeff(det(diagonal_matrix([G[j] for j in 1:i])), 0))
      push!(sym, s)
    end
  end
  return sym
end

################################################################################
#
#  Local isometry
#
################################################################################

function _islocally_isometric_kirschmer(L1::HermLat, L2::HermLat, p)
  R = base_ring(L1)
  E = nf(R)
  S1 = _genus_symbol_kirschmer(L1, p)
  S2 = _genus_symbol_kirschmer(L2, p)
  if !isdyadic(p) || !isramified(R, p)
    return S1 == S2
  end

  t = length(S1)
  if t != length(S2)
    return false
  end
  for i in 1:t
    for k in 1:4
      if S1[i][k] != S2[i][k]
        return false
      end
    end
  end

  if !islocal_norm(E, S1[t][5]//S2[t][5], p)
    return false
  end

  for i in 1:(t-1)
    @assert valuation(S1[i][5], p) == valuation(S2[i][5], p)
    x = S1[i][5]//S2[i][5]
    n = 2 * normic_defect(E, x, p)
    if n < (S1[i][4] + S1[i + 1][4]) - 2 * S1[i][2]
      return false
    end
  end
  return true
end


text:
#TODO: nonlinear smatrix routines, possibly using Optim?
"""
    S, flux_sct, flux_tot, absorption = smatrix(sim, k; channels, is_linear=true, F=[1], disp_opt=true, fileName="")

scattering matrix `S` at frequences `k` on `channels`
"""
function smatrix(sim::Simulation, k; is_linear=true, disp_opt=true, file_name="", is_parallel=nprocs()>1, num_channel_blocks=1)

    if is_linear && !is_parallel
        S, scattered_flux, total_flux, absorption = smatrix_l(sim, k, 1:length(sim.sct.channels), disp_opt, file_name)
    elseif is_linear && is_parallel
        S, scattered_flux, total_flux, absorption = smatrix_lp(sim, k, disp_opt, file_name, num_channel_blocks)
    else
        throw("no nonlinear method has yet been defined.")
    end

    return S, scattered_flux, total_flux, absorption
end

"""
    smatrix_l(sim, k, channels, F, disp_opt, file_name)
"""
function smatrix_l(sim::Simulation, k, channels, disp_opt, file_name)

    build_dispersions!(sim)

    nc = length(channels)
    nk = length(k)
    M = length(sim.sct.channels)

    S = Array{ComplexF64}(undef, nk, nc, M)
    flux_sct = Array{Float64}(undef, nk, nc)
    flux_tot = Array{Float64}(undef, nk, nc)
    absorption = Array{Float64}(undef, nk, nc)
    a = zeros(ComplexF64,M)

    if disp_opt
        pg = Progress(length(k)*M*nc, PROGRESS_UPDATE_TIME::Float64, "smatrix ")
    end

    for i ∈ eachindex(k)

        ζ = factorize_scattering_operator(sim, k[i])

        for m in 1:nc
            a = 0a
            a[channels[m]] = 1.
            ψ_sct, ψ_tot, ψ_out, ζ = scattering(sim, k[i], a; H=ζ, is_linear=true)
            for m′ ∈ 1:M
                if isempty(sim.sys.waveguides)
                    c = analyze_field(sim, k[i], ψ_out, m′)
                else
                    if sim.sct.channels[m′].waveguide == sim.sct.channels[m].waveguide
                        c = analyze_field(sim, k[i], ψ_sct, m′)
                    else
                        c = analyze_field(sim, k[i], ψ_tot, m′)
                    end
                end
                S[i, m, m′] = c
                if disp_opt
                    next!(pg)
                end
            end
            flux_sct[i, m] = surface_flux(sim, ψ_sct)[1][1]
            flux_tot[i, m] = surface_flux(sim, ψ_tot)[1][1]
            absorption[i, m] = bulk_absorption(sim, k[i], ψ_tot)[1]
        end

        if !isempty(file_name)
            fid = open(file_name,"w")
            serialize(fid, (S, flux_sct, flux_tot, absorption, sim, i))
            close(fid)
        end
    end

    return S, flux_sct, flux_tot, absorption
end

################################################################################
### S-MATRIX LINEAR ROUTINES
################################################################################

# end of fuction smatrix_l

#
# ################################################################################
# ### S-MATRIX NONLINEAR ROUTINES
# ################################################################################
# """
# S =  smatrix_nl(input::InputStruct; N=10, N_type="D", isNonLinear=false, F=1.,
#     dispOpt = true, ψ_init = [], fileName = "")
#
#     N is the number of steps to go from D0 = 0 to given D0
# """
# function smatrix_nl(input1::InputStruct, k::Array{ComplexF64,1};
#     N::Int=1, N_type::String="D", isNonLinear::Bool=false, F::Array{Float64,1}=[1.],
#     dispOpt::Bool=true, ψ_init::Array{ComplexF64,1}=ComplexF64[],
#     fileName::String = "")::Array{ComplexF64,4}
#
#     if !isempty(input1.bc ∩ ["o", "open", "pml_in"])
#         input = deepcopy(input1)
#         for i in 1:4
#             if input.bc[i] in ["o", "open", "pml_in"]
#                 input.bc[i] = "pml_out"
#                 updateInputs!(input, :bc, input.bc)
#             end
#         end
#     else
#         input = input1
#     end
#
#     M = length(input.channels)
#     D₀ = input.D₀
#     A = input.a
#
#     ψ₋ = Array{ComplexF64}(prod(input.N_ext))
#     ψ  = Array{ComplexF64}(prod(input.N_ext))
#
#     # if isempty(ψ_init) && (N>1)
#     S = NaN*ones(ComplexF64,length(input.k),N,M,M)
#     # else
#         # S = NaN*ones(ComplexF64,length(input.k),1,2,2)
#     # end
#
#     if N > 1
#         D = linspace(0,D₀,N)
#         A_vec = linspace(0.001,1,N)
#     else
#         D = [input.D₀]
#         A_vec = [1]
#     end
#
#     for ii in 1:length(input.k)
#
#         k = input.k[ii]
#
#         if (ii/1 == round(ii/1)) & dispOpt
#             if typeof(k)<:Real
#                 printfmtln("Solving for frequency {1:d} of {2:d}, ω = {3:2.3f}.",ii,length(input.k),k)
#             else
#                 printfmtln("Solving for frequency {1:d} of {2:d}, ω = {3:2.3f}{4:+2.3f}i.",ii,length(input.k),real(k),imag(k))
#             end
#         end
#
#         if isempty(ψ_init)
#
#             for j in 1:N
#
#                 if N_type == "D"
#                     updateInputs!(input, :D₀, D[j])
#                 elseif N_type == "A"
#                     updateInputs!(input, :a, A_vec[j]*A)
#                 end
#
#                 if isNonLinear
#                     if isnan(ψ[1]) | j==1
#                         ψ, W = computePsi(input,k,isNonLinear = true, F = F)
#                     else
#                         ψ, W = computePsi(input,k,isNonLinear = true, F = F, ψ_init = ψ)
#                     end
#                 else
#                     ψ = 0.
#                 end
#
#                 # compute linearized S-matrix
#                 ζ = lufact(speye(1,1)) # initialize lufact variable (speeds up repeated inhomogeneous solves)
#                 for m1 in 1:M
#                     # set flux-normalized input amplitudes to unity
#                     at = zeros(ComplexF64,M)
#                     at[m1] = 1.
#                     updateInputs!(input,:a,at)
#                     # solve for scattered field
#                     ψ₋, dummy1, ζ, inputs_s = computePsi(input, k; isNonLinear=false, F=F./(1+abs2.(γ(k,input)*ψ)), A=ζ)
#                     # analyze into channels
#                     println("here 1")
#                     for m2 in m1:M
#                         dt=0
#                         if input.channelBoundaries[m2] in [1,2]
#                             t = inputs_s.x₂_ext
#                             u = inputs_s.x₁_inds
#                             dt = input.dx̄[2]
#                         else
#                             t = inputs_s.x₁_ext
#                             u = inputs_s.x₂_inds
#                             dt = input.dx̄[1]
#                         end
#                         println("here 2")
#                         ϕ = zeros(ComplexF64,length(t))
#                         for q in 1:length(t)
#                              ϕ[q] = conj(input.incidentWave(k, m2, input.∂R[input.channelBoundaries[m2]], t[ii], input.∂R, input.geometry, input.geoParams)[1])
#                         end
#                         println("here 3")
#                         P = 0
#                         if input.channelBoundaries[m2] in [1,3]
#                             println("here 4")
#                             println(size(ψ₋))
#                             println(size(t))
#                             println(size(u))
#                             P = reshape(ψ₋,:,length(t))[u[1],:]
#                             println("here 5")
#                         else
#                             println("here 6")
#                             P = reshape(ψ₋,:,length(t))[u[end],:]
#                         end
#                         println("here 7")
#                         println(size(P))
#                         println(size(S))
#                         println(j)
#                         println(ii)
#                         println(size(ϕ))
#                         S[ii,j,m1,m2] = sum(P.*ϕ)*dt
#                         S[ii,j,m2,m1] = S[ii,j,m1,m2]
#                     end
#                 end
#     println("here 8")
#                 updateInputs!(input,:D₀, D₀)
#                 updateInputs!(input,:a , A )
#
#                 if !isempty(fileName)
#                     fid = open(fileName,"w")
#                     serialize(fid,(input,D,S,ii,j))
#                     close(fid)
#                 end
#
#             end
#
#         else
#             if isNonLinear
#                 ψ = scattering(input,k,isNonLinear = true, F = F, ψ_init = ψ_init)
#             else
#                 ψ = 0.
#             end
#
#             input.a = [1.0,0.0]
#             ψ₊, W = scattering(input,k,isNonLinear = false, F = F./(1+abs2.(γ(input,k)*ψ)))
#
#             input.a = [0.0,1.0]
#             ψ₋, W = scattering(input,k,isNonLinear = false, F = F./(1+abs2.(γ(input,k)*ψ)), A = W)
#
#             S[1,1,ii,1] = ψ₊[x_inds[1]]*exp(+1im*input.dx*k)
#             S[2,1,ii,1] = ψ₊[x_inds[end]]
#             S[1,2,ii,1] = ψ₋[x_inds[1]]
#             S[2,2,ii,1] = ψ₋[x_inds[end]]*exp(-1im*input.dx*k)
#
#             input.D₀ = D₀
#             input.a = A
#
#             if !isempty(fileName)
#                 foo2(fid) = serialize(fid,(input,D,S,ii,1))
#                 open(foo2,fileName,"w")
#             end
#
#         end
#
#     end
#     println("here 9")
#     return S
# end # end of fuction smatrix_nl


text:
if length(ARGS)∉(7,8)
    println("usage: (assuming you are running from DFQExamples.jl directory)")
    println("  julia -p n_procs empirical_cdfs.jl ARGS[1] ... ARGS[7] [ARGS[8]]")
    println("  where:")
    println("    n_procs: the number of workers to run (and hence number of output files as each worker writes its own data)")
    println("             n_procs+1 cpus will be used in total")
    println("    ARGS[1]: path to include model definition")
    println("    ARGS[2]: path to include simulation definitions (initial_condition_generator, stopping time)")
    println("    ARGS[3]: n_sims, number of sims to do")
    println("    ARGS[4]: n_boot, number of bootstrap resamples to do")
    println("    ARGS[5]: seed for rng")
    println("    ARGS[6]: filepath to write CDF to")
    println("    ARGS[7]: filepath to write bootstrap sample of CDFs to")
    println("    ARGS[8]: (opional) true/false, to specify whether to record time or position of fluid at the stopping time of the sim, default false")
    throw(ArgumentError("Wrong number of arguments"))
end

using Distributed
const n_procs = nprocs()-1
@show ps = procs()

@show pwd()
@everywhere begin
    using Pkg; Pkg.activate(".")
    using DiscretisedFluidQueues
    using FluidFluidQueues
    using Random
    using StatsBase
    using JSON
    using TickTock
end

@everywhere include($(ARGS[1])) # define model
@everywhere include($(ARGS[2])) # define initial_condition_generator, stopping time, and points to evaluate the CDF,
const n_sims = parse(Int, ARGS[3]) # define number of sims to do
const n_boot = parse(Int, ARGS[4]) # number of bootstrap sample to do
const seed = parse(Int, ARGS[5]) # set rng seed
const rng = MersenneTwister(seed)
const record_time = length(ARGS)>7 ? parse(Bool, ARGS[8]) : false # specify whether to record time or position of sim

@everywhere begin
    function _simulate(model::FluidFluidQueue, tau, ic)
        rest, ~ = simulate(model, tau, ic)
        X = rest.X
        φ = Int.(rest.φ)
        t = rest.t
        return (X=X, φ=φ, t=t)
    end
    _simulate(model, tau, ic) = simulate(model, tau, ic)
end
@everywhere function empirical_cdf(rng, ic_generator, model, n_sims, tau, x::AbstractRange, record_time=false)
    println("Simulating...")
    cdf = zeros(Int, length(x)+1, typeof(model)===BoundedFluidQueue ? 2 : 4)
    f, h = first(x), step(x)
    n_laps = 0
    total_laps = n_sims÷100_000_000
    mean_lap_time = 0.0
    tick()
    t0 = peektimer()
    for n in 1:n_sims
        (n%500_000_000 == 0) && begin
            n_laps += 1
            t1 = peektimer()
            @info "msg @ time "*string(TickTock.now())
            println("  Simulating... progress on this worker: ", n/n_sims * 100, "%, elapsed: ", t1, " secs, last lap: ", t1-t0, " secs.")
            mean_lap_time = ((n_laps-1)*mean_lap_time + t1-t0)/(n_laps)
            println("  Simulating... predicted time remaining on this worker: ", (total_laps-n_laps)*mean_lap_time, " seconds.")
            println()
            t0 = t1
        end
        # simulate from model
        ic = ic_generator(rng,1,model)
        sim_tau = _simulate(model, tau, ic)
        Xτ, φτ, tτ = only(sim_tau.X), only(sim_tau.φ), only(sim_tau.t)
        quantity_of_interest = record_time ? tτ : Xτ
        # add to cdf
        idx = ceil(Int, (quantity_of_interest - f) / h)+1
        # idx = Xτ<x[idx] ? idx-1 : idx
        cdf[idx, φτ] += 1
    end
    @info "msg @ time "*string(TickTock.now())
    tock()
    return cumsum(cdf,dims=1)./n_sims
end

@everywhere begin
    @inline function batch_boot!(cdf,rng,pdf_ish,n)
        idxs = wsample(rng, 1:length(pdf_ish), pdf_ish[:], n)
        for i in idxs
            cdf[i] += 1
        end
        return nothing
    end

    function boot1(rng,pdf_ish,n)
        cdf_boot = zeros(Int, size(pdf_ish))
        # process in batches of 1_000_000 samples
        for _ in 1:(n÷1_000_000)
            batch_boot!(cdf_boot, rng, pdf_ish, 1_000_000)
        end
        batch_boot!(cdf_boot, rng, pdf_ish, n%1_000_000)
        return cumsum(cdf_boot, dims=1)./n
    end

    function boot(filepath, rng, cdf_paths, n_boot, ssize)
        jsons = [JSON.parsefile(pth) for pth in cdf_paths]
        cdf = zeros(length(jsons[1][1]),length(jsons[1]))
        for data in jsons
            cdf .+= hcat(data...)
        end
        cdf ./= length(jsons)
        println("Bootstrapping...")
        boot_samples = Array{Array{Array{Float64,1},1},1}(undef,n_boot)
        pdf_ish = [cdf[1,:]'; diff(cdf,dims=1)]
        n_laps = 0
        total_laps = n_boot÷50
        mean_lap_time = 0.0
        tick()
        t0 = peektimer()
        for n in 1:n_boot
            (n%5 == 0) && begin
                n_laps += 1
                t1 = peektimer()
                @info "msg @ time "*string(TickTock.now())
                println("  Bootstrapping... progress on this worker: ", n/n_boot * 100, "%, elapsed: ", t1, " secs, last lap: ", t1-t0, " secs.")
                mean_lap_time = ((n_laps-1)*mean_lap_time + t1-t0)/(n_laps)
                println("  Bootstrapping... predicted time remaining on this worker: ", (total_laps-n_laps)*mean_lap_time, " seconds.")
                println()
                t0 = t1
            end
            samp = boot1(rng, pdf_ish, ssize)
            boot_samples[n] = [samp[:,i] for i in 1:size(cdf,2)]
        end
        tock()
        @info "msg @ time "*string(TickTock.now())
        println("Writing bootstrap on this worker...")
        JSON.write(filepath,JSON.json(boot_samples))
        println("...write complete.")
        return nothing
    end
end

@everywhere function sim_and_write(write_path, rng, initial_condition_generator, model, n_sims, tau, x, record_time)
    ecdf = empirical_cdf(rng, initial_condition_generator, model, n_sims, tau, x, record_time)
    JSON.write(write_path, JSON.json(ecdf))
    return nothing
end

@time @sync begin
    for n in n_procs+1:-1:3
        local_rng = MersenneTwister(rand(rng, UInt32)) # make sure rngs on each thread are distinct
        local_write_file = "proc_"*string(ps[n])*"_"*ARGS[6]
        @spawnat ps[n] sim_and_write(local_write_file, rng, initial_condition_generator, model, n_sims÷n_procs, tau, x, record_time)
    end
    local_rng = MersenneTwister(rand(rng, UInt32)) # make sure rngs on each thread are distinct
    local_write_file = "proc_"*string(ps[2])*"_"*ARGS[6]
    @spawnat ps[2] sim_and_write(local_write_file, local_rng, initial_condition_generator, model, n_sims-(n_procs-1)*(n_sims÷n_procs), tau, x, record_time)
end

local_cdf_paths = "proc_".*string.(ps[2:end]).*"_".*ARGS[6]
println("Bootstrapping on workers...")
@time @sync begin
    for n in n_procs+1:-1:3
        local_rng = MersenneTwister(rand(rng, UInt32)) # make sure rngs on each thread are distinct
        local_write_file = "proc_"*string(ps[n])*"_"*ARGS[7]
        @spawnat ps[n] boot(local_write_file, local_rng, local_cdf_paths, n_boot÷n_procs, n_sims)
    end
    local_rng = MersenneTwister(rand(rng, UInt32)) # make sure rngs on each thread are distinct
    local_write_file = "proc_"*string(ps[2])*"_"*ARGS[7]
    @spawnat ps[2] boot(local_write_file, local_rng, local_cdf_paths, n_boot-(n_procs-1)*(n_boot÷n_procs), n_sims)
end
println("...workers complete.")


text:
#############################
# An abstract subdictionary
#############################

# We start with a generic description of subdicts and imlement as much
# functionality as possible at the abstract level of Subdictionary.
# This is followed by specialized concrete types of subdicts.

"""
A `Subdictionary` is an abstract subset of a dictionary. It is characterized by
the underlying larger dictionary, and a subcollection of its indices.
"""
abstract type Subdictionary{S,T} <: Dictionary{S,T}
end

# We assume that the underlying set is stored in a field called superdict
superdict(dict::Subdictionary) = dict.superdict

# We assume that the underlying indices are stored in a field called superindices
superindices(dict::Subdictionary) = dict.superindices
superindices(dict::Subdictionary, idx::Int) = dict.superindices[idx]
superindices(dict::Subdictionary, idx::DefaultNativeIndex) = superindices(dict, value(idx))

function similar(d::Subdictionary, ::Type{T}, n::Int) where {T}
    @assert n == length(d)
    similar_subdict(d, similar(superdict(d), T), superindices(d))
end

mapped_dict(s::Subdictionary, map) = similar_subdict(s, mapped_dict(superdict(s), map), superindices(s))

show(io::IO, mime::MIME"text/plain", d::Subdictionary) = composite_show(io, mime, d)
Display.object_parentheses(d::Subdictionary) = false
Display.stencil_parentheses(d::Subdictionary) = false
Display.displaystencil(d::Subdictionary) = [superdict(d), "[", superindices(d), "]"]

support(s::Subdictionary) = support(superdict(s))

support(s::Subdictionary, i) = support(superdict(s),superindices(s,i))

size(s::Subdictionary) = size(superindices(s))

##
# For various properties it is not a priori clear whether a subdict has them or not,
# even if the underlying set does. For each function foo, we invoke a new function
# subdict_foo, that has the underlying set and the range of indices as additional
# arguments. This function call can be intercepted and implemented for specific kinds
# of subdicts, if it is known that foo applies.

resize(s::Subdictionary, n) = subdict_resize(s, n, superdict(s), superindices(s))

hasinterpolationgrid(s::Subdictionary) = subdict_hasinterpolationgrid(s, superdict(s), superindices(s))
hasderivative(s::Subdictionary) = subdict_hasderivative(s, superdict(s), superindices(s))
hasantiderivative(s::Subdictionary) = subdict_hasantiderivative(s, superdict(s), superindices(s))
hastransform(s::Subdictionary) = subdict_hastransform(s, superdict(s), superindices(s))

derivative_dict(s::Subdictionary, order; options...) = subdict_derivative_dict(s, order, superdict(s), superindices(s); options...)
antiderivative_dict(s::Subdictionary, order; options...) = subdict_antiderivative_dict(s, order, superdict(s), superindices(s); options...)

interpolation_grid(s::Subdictionary) = subdict_interpolation_grid(s, superdict(s), superindices(s))

subdict_interpolation_grid(s::Subdictionary, gb::GridBasis, indices) =
    interpolation_grid(gb)[indices]

# By default, we have none of the properties
subdict_hasderivative(s::Subdictionary, superdict, superindices) = false
subdict_hasantiderivative(s::Subdictionary, superdict, superindices) = false
subdict_hastransform(s::Subdictionary, superdict, superindices) = false
subdict_hasinterpolationgrid(s::Subdictionary, superdict, superindices) = false

for op in (:isreal, :isbasis)
    @eval $op(dict::Subdictionary) = $op(superdict(dict))
end

for op in (:isreal, :isbiorthogonal, :isorthonormal, :isorthogonal)
    @eval $op(dict::Subdictionary, measure::Weight) = $op(superdict(dict), measure::Weight)
end

for op in (:moment, :norm)
    @eval $op(d::Subdictionary, i) = $op(superdict(d), superindices(d, i))
end

dict_in_support(d::Subdictionary, i, x) = in_support(superdict(d), superindices(d, i), x)

unsafe_eval_element(d::Subdictionary, i, x) = unsafe_eval_element(superdict(d), superindices(d, i), x)

unsafe_eval_element_derivative(d::Subdictionary, i, x, order) = unsafe_eval_element_derivative(superdict(d), superindices(d, i), x, order)

hasmeasure(dict::Subdictionary) = hasmeasure(superdict(dict))
measure(dict::Subdictionary) = measure(superdict(dict))

innerproduct1(d1::Subdictionary, i, d2, j, measure; options...) =
    innerproduct(superdict(d1), superindices(d1,i), d2, j, measure; options...)
innerproduct2(d1, i, d2::Subdictionary, j, measure; options...) =
    innerproduct(d1, i, superdict(d2), superindices(d2, j), measure; options...)

#####################
# Concrete subdicts
#####################

"""
A `DenseSubdict` is a large subset of a dictionary. Operators associated with the
subset are implemented in terms of corresponding operators on the underlying dictionary.
This often leads to an explicit extension to the full set, but it can take
advantage of possible fast implementations for the underlying dictionary. For
large subsets this is more efficient than iterating over the individual elements.
"""
struct DenseSubdict{SET,IDX,S,T} <: Subdictionary{S,T}
    superdict       ::  SET
    superindices    ::  IDX

    function DenseSubdict{SET,IDX,S,T}(dict::Dictionary{S,T}, idx::IDX) where {SET,IDX,S,T}
        checkbounds(dict, idx)
        new(dict, idx)
    end
end

DenseSubdict(dict::Dictionary{S,T}, superindices) where {S,T} =
    DenseSubdict{typeof(dict),typeof(superindices),S,T}(dict, superindices)

name(dict::DenseSubdict) = "Dense subdictionary"

similar_subdict(d::DenseSubdict, dict, superindices) = DenseSubdict(dict, superindices)

interpolation_grid(d::DenseSubdict) = interpolation_grid(superdict(d))

function extension(::Type{T}, src::DenseSubdict, dest::Dictionary; options...) where {T}
    @assert dest == superdict(src)
    IndexExtension{T}(src, dest, superindices(src); options...)
end

function restriction(::Type{T}, src::Dictionary, dest::DenseSubdict; options...) where {T}
    @assert src == superdict(dest)
    IndexRestriction{T}(src, dest, superindices(dest); options...)
end

# In general, the derivative set of a subdict can be the whole derivative set
# of the underlying set. We can not know generically whether the derivative set
# can be indexed as well. Same for antiderivative.
# Yet, we can generically define a differentiation operator by extending the subdict
# to the whole set and then invoking the differentiation operator of the latter,
# and we choose that to be the default.
subdict_hasderivative(s::DenseSubdict, superdict, superindices) = hasderivative(superdict)
subdict_hasantiderivative(s::DenseSubdict, superdict, superindices) = hasantiderivative(superdict)

function subdict_derivative_dict(Φ::DenseSubdict, order, superdict, superindices; options...)
    # We test explicitly for zero order derivative, in order to avoid changing
    # the dictionary in this case.
    if orderiszero(order)
        Φ
    else
        derivative_dict(superdict, order; options...)
    end
end
function subdict_antiderivative_dict(Φ::DenseSubdict, order, superdict, superindices; options...)
    if orderiszero(order)
        Φ
    else
        antiderivative_dict(superdict, order; options...)
    end
end

function differentiation(::Type{T}, src::DenseSubdict, dest::Dictionary, order; options...) where {T}
    @assert dest == derivative_dict(src, order)
    D = differentiation(T, superdict(src), dest, order; options...)
    E = extension(T, src, superdict(src); options...)
    D*E
end

function antidifferentiation(::Type{T}, src::DenseSubdict, dest::Dictionary, order; options...) where {T}
    @assert dest == antiderivative_dict(src, order)
    D = antidifferentiation(T, superdict(src), dest, order; options...)
    E = extension(T, src, superdict(src); options...)
    D*E
end

"""
A `SparseSubdict` is a subset of a dictionary with a small number of indices.
The difference with a regular function subset is that operators on a small set
are implemented by iterating explicitly over the indices, and not in terms of
an operator on the full underlying set.
"""
struct SparseSubdict{SET,IDX,S,T} <: Subdictionary{S,T}
    superdict       ::  SET
    superindices    ::  IDX

    function SparseSubdict{SET,IDX,S,T}(dict::Dictionary{S,T}, indices::IDX) where {SET,IDX,S,T}
        checkbounds(dict, indices)
        new(dict, indices)
    end
end

SparseSubdict(dict::Dictionary{S,T}, superindices) where {S,T} =
    SparseSubdict{typeof(dict),typeof(superindices),S,T}(dict, superindices)

name(dict::SparseSubdict) = "Sparse subdictionary"

similar_subdict(d::SparseSubdict, dict, superindices) = SparseSubdict(dict, superindices)

#########################################
# The default logic of creating subdicts
#########################################

# First some additions to the checkbounds ecosystem

checkbounds(::Type{Bool}, dict::Dictionary, indices::CartesianIndices) =
    checkbounds(Bool, dict, first(indices)) && checkbounds(Bool, dict, last(indices))

function checkbounds(::Type{Bool}, dict::Dictionary, indices::Array)
    result = true
    for idx in indices
        result &= checkbounds(Bool, dict, idx)
    end
    result
end

# We intercept getindex only if we are sure the indices are a set. We have no
# way of knowing what are all the user defined index types.
getindex(dict::Dictionary, idx::AbstractUnitRange) = sub(dict, idx)
getindex(dict::Dictionary, idx::CartesianIndices) = sub(dict, idx)
getindex(dict::Dictionary, idx::AbstractArray) = sub(dict, idx)
getindex(dict::Dictionary, idx::Colon) = sub(dict, idx)

# By default we generate a large subdict
sub(dict::Dictionary, superindices) = defaultsub(dict, superindices)
defaultsub(dict::Dictionary, superindices) = DenseSubdict(dict, superindices)

# For an array, which as far as we know has no additional structure, we create
# a small subdict by default. For large arrays, a large subdict may be more efficient.
sub(dict::Dictionary, superindices::Array) = SparseSubdict(dict, superindices)

sub(dict::Dictionary, ::Colon) = dict

# Avoid creating nested subdicts
sub(dict::Subdictionary, idx) = sub(superdict(dict), superindices(dict)[idx])

# A gram operator of a subdict orthonormal/orthogonal dictionary is also Identity/Diagonal
function gram1(T, s::Subdictionary, m; options...)
    if isorthonormal(s, m)
        return IdentityOperator{T}(s)
    elseif isorthogonal(s, m)
        return DiagonalOperator{T}(s, diag(gram(T, superdict(s), m; options...))[superindices(s)])
    end
    default_gram(T, s, m; options...)
end


text:
export SHalley_dir_AD

function SHalley_dir_AD(nlp            :: AbstractNLPModel,
                       nlp_at_x       :: NLPAtX,
                       x              :: AbstractVector,
                       Nwtdirection   :: Function)

    n = length(x)

    dₙ = Nwtdirection(nlp_at_x.Hx, nlp_at_x.gx, verbose = false)

    e₁ = zeros(n)
    e₁[1] = 1.0
    ∇³fdₙ = ∇f³xuv(nlp, x, dₙ, e₁)
    for i = 2:n
        eᵢ = zeros(n)
        eᵢ[i] = 1.0
        ∇³fdₙ = hcat(∇³fdₙ, ∇f³xuv(nlp, x, dₙ, eᵢ))
    end

    hess_plus_tenseur = (nlp_at_x.Hx + ∇³fdₙ)

    dV = Nwtdirection(hess_plus_tenseur, nlp_at_x.gx)

    grad_plus_tuv = -(nlp_at_x.gx + 0.5 .* ∇f³xuv(nlp, x, dₙ, dV))

    dSH = Nwtdirection(nlp_at_x.Hx, -grad_plus_tuv)

    return dSH
end