text:
#Flip x and y axis. Currently only works for XY plots
#Doesn't work for scatter, boxplot brittle
function flip!(ec::EChart; rotatedims::Bool = false)

    ec.xAxis, ec.yAxis = ec.yAxis, ec.xAxis

	#Rotate dimensions to keep same aspect ratio
	if rotatedims
		ec.ec_width, ec.ec_height = ec.ec_height, ec.ec_width
	end

	# If boxplot, invert outlier values
	# This could be more robust to not use numeric indices
	ec.ec_charttype == "box" &&  length(ec.series) > 1 ? ec.series[2].data = [[x[2], x[1]] for x in ec.series[2].data] : nothing

    return ec

end


text:
# Robust Optimization
# Version: 2.0
# Author: Edward J. Xu, edxu96@outlook.com
# Date: April 8th, 2019
# 0. Set working path, install packages and use module -----------------------------------------------------------------
push!(LOAD_PATH, "$(homedir())/Desktop/Production Planning, Robust Optimization, DTU02435")
cd("$(homedir())/Desktop/Production Planning, Robust Optimization, DTU02435")
using ExcelReaders
using RobustOptim
using LinearAlgebra
using JuMP
using GLPKMathProgInterface
# 1. Import Original Data ----------------------------------------------------------------------------------------------
include("data.jl")
# 2. Tranform to standard form -----------------------------------------------------------------------------------------
# Question 1: MILP with Box Uncertainty and Budget of Uncertainty
function sectionBudgetMILP(gammaCap)
    num_x = 40
    num_y = 4
    vec_min_y = hcat(zeros(4))
    vec_max_y = hcat(repeat([100], num_y))
    (vec_c, vec_f, vec_b, mat_a, mat_b, mat_a_bar, mat_a_hat, mat_b_bar, vec_b_bar) = paras_budgetMILP(num_x, num_y)
    # Define vec_gammaCap
    vec_gammaCap = zeros(4)
    for m = 1:4
        vec_gammaCap[m] = gammaCap * sum([mat_prodOnMachine[p, m]] for p = 1: 10)[1]
    end
    # Solve the model
    (obj, vec_y, vec_x, vec_result_z, mat_mu, vec_lambda) = RobustOptim.milpBoxBudget(
        num_x=num_x, num_y=num_y, vec_min_y=vec_min_y, vec_max_y=vec_max_y, vec_c=vec_c, vec_f=vec_f,
        vec_b=vec_b, mat_a=mat_a, mat_b=mat_b, mat_a_bar=mat_a_bar, mat_a_hat=mat_a_hat, mat_b_bar=mat_b_bar,
        vec_b_bar=vec_b_bar, vec_gammaCap=vec_gammaCap)
    # return ((obj, vec_y, vec_x, vec_result_z, mat_mu, vec_lambda)
end
# Question 2: Two-Stage Stochastic LP with Box Uncertainty
function sectionStochasticLP(penalty)
    num_x = 40
    num_y = 4
    vec_min_y = hcat(zeros(4))
    vec_max_y = hcat(repeat([100], num_y))
    (vec_c, vec_f, vec_b, mat_a, mat_b, mat_a_bar, mat_a_hat, mat_b_bar, vec_b_bar, vec_g, mat_d, mat_d_bar) =
        paras_stochasticLP(num_x, num_y, penalty)
    # Solve the model
    (obj, vec_y, vec_x, vec_gamma, vec_theta1, vec_theta2) =
        RobustOptim.milpAdjustBox(num_x=num_x, num_y=num_y, num_z=10, vec_min_y=vec_min_y, vec_max_y=vec_max_y,
        vec_c=vec_c, vec_f=vec_f, vec_g=vec_g, vec_b=vec_b, mat_a=mat_a, mat_b=mat_b, mat_d=mat_d,
        mat_a_bar=mat_a_bar, mat_a_hat=mat_a_hat, mat_b_bar=mat_b_bar, mat_d_bar=mat_d_bar, vec_b_bar=vec_b_bar)
    # return (obj, vec_y, vec_x, vec_gamma, vec_theta1, vec_theta2)
end

function main()
    # Question 1: MILP with Box Uncertainty and Budget of Uncertainty
    println("#########################################################\n",
            "   MILP with Box Uncertainty and Budget of Uncertainty\n",
            "#########################################################\n",)
    vec_gammaCap = [0.3, 0, 1]
    for i = 1: 3
        println("*********************************************************\n",
                "   When gammaCap = $(vec_gammaCap[i])\n",
                "*********************************************************\n",)
        sectionBudgetMILP(vec_gammaCap[i])
    end
    # Question 2: Two-Stage Stochastic LP with Box Uncertainty
    println("#########################################################\n",
            "   Two-Stage Stochastic LP with Box Uncertainty\n",
            "#########################################################\n",)
    vec_penalty = [1, 2, 10, 20, 50, 100, 200, 300, 500]
    for i = 1: length(vec_penalty)
        println("*********************************************************\n",
                "   When penalty = $(vec_penalty[i])\n",
                "*********************************************************\n")
        sectionStochasticLP(vec_penalty[i])
    end
end

main()

ext: jl


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the Jacchia-Robert 1971 model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/earth/atmospheric_models/jr1971
# ===========================================

# Functions: jr1971
# -----------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
#   Values obtained from GMAT R2018a using the following inputs:
#
#       Date:      2017-01-01 00:00:00 UTC
#       Latitude:  45 deg
#       Longitude: 0 deg
#       F10.7:     100
#       F10.7ₐ:    100
#       Kp:        4
#
#   Result:
#
#       | Altitude [km]  | Density [g/cm³] |
#       |----------------|------------------|
#       |             92 | 3.87506e-09      |
#       |            100 | 3.96585e-09      |
#       |          100.1 | 6.8354e-10       |
#       |          110.5 | 1.2124e-10       |
#       |            125 | 1.60849e-11      |
#       |          125.1 | 1.58997e-11      |
#       |            300 | 1.30609e-14      |
#       |            700 | 1.34785e-17      |
#       |           1500 | 4.00464e-19      |
#
# Scenario 02
# ===========
#
#   Values obtained from GMAT R2018a using the following inputs:
#
#       Date:      2017-01-01 00:00:00 UTC
#       Latitude:  45 deg
#       Longitude: 0 deg
#       F10.7:     100
#       F10.7ₐ:    100
#       Kp:        1
#
#   Result:
#
#       | Altitude [km]  | Density [g/cm³] |
#       |----------------|------------------|
#       |             92 | 3.56187e-09      |
#       |            100 | 3.65299e-09      |
#       |          100.1 | 6.28941e-10      |
#       |          110.5 | 1.11456e-10      |
#       |            125 | 1.46126e-11      |
#       |          125.1 | 1.44428e-11      |
#       |            300 | 9.08858e-15      |
#       |            700 | 8.51674e-18      |
#       |           1500 | 2.86915e-19      |
#
# Scenario 03
# ===========
#
#   Values obtained from GMAT R2018a using the following inputs:
#
#       Date:      2017-01-01 00:00:00 UTC
#       Latitude:  45 deg
#       Longitude: 0 deg
#       F10.7:     100
#       F10.7ₐ:    100
#       Kp:        9
#
#   Result:
#
#       | Altitude [km]  | Density [g/cm³] |
#       |----------------|------------------|
#       |             92 | 5.55597e-09      |
#       |            100 | 5.63386e-09      |
#       |          100.1 | 9.75634e-10      |
#       |          110.5 | 1.73699e-10      |
#       |            125 | 2.41828e-11      |
#       |          125.1 | 2.39097e-11      |
#       |            300 | 3.52129e-14      |
#       |            700 | 1.28622e-16      |
#       |           1500 | 1.97775e-18      |
#
################################################################################

@testset "Function jr1971" begin
    # Common inputs to all scenarios.
    JD   = DatetoJD(2017,1,1,0,0,0)
    glat = 45*pi/180
    glon = 0.0
    F10  = 100
    F10ₐ = 100
    h    = [92; 100; 100.1; 110.5; 125; 125.1; 300; 700; 1500]*1000

    # Scenario 01
    # ===========

    Kp = 4

    # Results in [kg/m³]
    results = [3.87506e-09;
               3.96585e-09;
               6.83540e-10;
               1.21240e-10;
               1.60849e-11;
               1.58997e-11;
               1.30609e-14;
               1.34785e-17;
               4.00464e-19]*1000

    for i = 1:length(h)
        ret = jr1971(JD, glat, glon, h[i], F10, F10ₐ, Kp)
        @test ret.rho ≈ results[i] rtol = 5e-4
    end

    # Scenario 02
    # ===========

    Kp = 1

    # Results in [kg/m³]
    results = [3.56187e-09;
               3.65299e-09;
               6.28941e-10;
               1.11456e-10;
               1.46126e-11;
               1.44428e-11;
               9.08858e-15;
               8.51674e-18;
               2.86915e-19]*1000

    for i = 1:length(h)
        ret = jr1971(JD, glat, glon, h[i], F10, F10ₐ, Kp)
        @test ret.rho ≈ results[i] rtol = 5e-4
    end

    # Scenario 03
    # ===========

    Kp = 9

    # Results in [kg/m³]
    results = [5.55597e-09
               5.63386e-09
               9.75634e-10
               1.73699e-10
               2.41828e-11
               2.39097e-11
               3.52129e-14
               1.28622e-16
               1.97775e-18]*1000

    for i = 1:length(h)
        ret = jr1971(JD, glat, glon, h[i], F10, F10ₐ, Kp)
        @test ret.rho ≈ results[i] rtol = 5e-4
    end
end


text:
module Calibration

using StaticArrays

export
    CalibrationData,
    fit

struct CalibrationData{T<:AbstractFloat}
    cnt::Int   # number of samples
    avg::T     # average value
    sqr::T     # average squared value
    fct::T     # brightness factor, e.g. 1 for flat field, 0 for dark or bias
    Δt::T      # exposure time
end

"""

The parameters are `[a, b, c, g, σ]` and the expectation and variance of
the data are given by:

    E(d|x) = (s*a + c)*Δt/g + b
    Var(d|x) = ((s*a + c)*Δt + σ²)/g²

with `d` the raw pixel data, `s = 1` for flat images and `s = 0` for bias and
dark images, `a` the flat flux times the quantum efficiency, `c` the dark
current, `Δt` the exposure time, `b` the bias, `g` the conversion gain and `σ`
the standard deviation of the readout noise.  Here, `a` and `c` are in
photo-electrons per pixel per unit of time, the gain `g` is en photo-electrons
per ADU (analog-digital units), `b` is in ADU per pixel per frame, and `σ` is
in photo-electrons per pixel per frame.

A more convenient set of parameters for the fit `[a',b,c',g,v']` where `a' =
a/g`, `c' = c/g` and `v' = σ²/g`

"""
function fit(::Type{Val{:detector_model_abc}},
             dat::AbstractVector{CalibrationData{T}},
             wgt::AbstractVector{T} = ones(T, length(dat));
             nonnegative::Bool = false) where {T<:AbstractFloat}

    n = length(dat)
    @assert length(wgt) == n
    A11 = A21 = A31 = zero(T)
    A22 = A32 = zero(T)
    A33 = zero(T)
    b1 = b2 = b3 = zero(T)
    γ = zero(T)
    @inbounds @simd for i in eachindex(dat, wgt)
        d = dat[i]
        w = T(d.cnt)*wgt[i] # total weight

        # Basis functions.
        h1 = d.fct*d.Δt
        h2 = one(T)
        h3 = d.Δt

        # Update coefficients of LHS matrix.
        h1_w, h2_w, h3_w = h1*w, h2*w, h3*w
        A11 += h1_w*h1
        A21 += h2_w*h1
        A31 += h3_w*h1
        A22 += h2_w*h2
        A32 += h3_w*h2
        A33 += h3_w*h3

        # Update coefficients of RHS vector.
        b1 += h1_w*d.avg
        b2 += h2_w*d.avg
        b3 += h3_w*d.avg

        # Update constant term.
        γ += w*d.sqr
    end

    # Solve the normal equations (without constraints and using fast "static
    # arrays").
    q0, s0 = solve3(A11, A21, A31, A22, A32, A33, b1, b2, b3)
    qmax, x1, x2, x3 = q0, s0[1], s0[2], s0[3]

    # Check whether solving the unsconstrained problem yields a feasible
    # solution.
    if nonnegative && (x1 < zero(T) || x3 < zero(T))
        # Assuming inactive constraints yields an unfeasible model.  We
        # have to determine which constraints are active.  We first assume
        # that both constraints are active.
        q1, s1 = solve1(A22, b2)
        qmax, x1, x2, x3 = q1, zero(T), s1[1], zero(T)
        # Second, assume only x1 is 0.
        q2, s2 = solve2(A22, A32, A33, b2, b3)
        if q2 > qmax && s2[2] ≥ zero(T)
            qmax, x1, x2, x3 = q2, zero(T), s2[1], s2[2]
        end
        # Third, assume only x3 is 0.
        q3, s3 = solve2(A11, A21, A22, b1, b2)
        if q3 > qmax && s3[1] ≥ zero(T)
            qmax, x1, x2, x3 = q2, s3[1], s3[2], zero(T)
        end
    end
    return γ - qmax/2, (x1, x2, x3)
end

# Idem but assume dark-current is zero.
function fit(::Type{Val{:detector_model_ab}},
             dat::AbstractVector{CalibrationData{T}},
             wgt::AbstractVector{T} = ones(T, length(dat));
             nonnegative::Bool = false) where {T<:AbstractFloat}

    n = length(dat)
    @assert length(wgt) == n
    A11 = A21 = zero(T)
    A22 = zero(T)
    b1 = b2 = zero(T)
    γ = zero(T)
    @inbounds @simd for i in eachindex(dat, wgt)
        d = dat[i]
        w = T(d.cnt)*wgt[i] # total weight

        # Basis functions.
        h1 = d.fct*d.Δt
        h2 = one(T)

        # Update coefficients of LHS matrix.
        h1_w, h2_w = h1*w, h2*w
        A11 += h1_w*h1
        A21 += h2_w*h1
        A22 += h2_w*h2

        # Update coefficients of RHS vector.
        b1 += h1_w*d.avg
        b2 += h2_w*d.avg

        # Update constant term.
        γ += w*d.sqr
    end

    # Solve the normal equations (without constraints and using fast "static
    # arrays").
    q0, s0 = solve2(A11, A21, A22, b1, b2)
    qmax, x1, x2, x3 = q0, s0[1], s0[2], zero(T)

    # Check whether solving the unsconstrained problem yields a feasible
    # solution.
    if nonnegative && x1 < zero(T)
        # Assuming inactive constraints yields an unfeasible model.  The only
        # other possibility is to have x1 = 0.
        q1, s1 = solve1(A22, b2)
        qmax, x1, x2, x3 = q1, zero(T), s1[1], zero(T)
    end
    return γ - qmax/2, (x1, x2, x3)
end

#------------------------------------------------------------------------------
# UTILITIES

"""
```julia
solve1(A11, b1) -> q, x
solve2(A11, A21, A22, b1, b2) -> q, x
solve3(A11, A21, A31, A22, A32, A33, b1, b2, b3) -> q, x
```

yields the maximal score `q = 2b'⋅x - x'⋅A⋅x` and the corresponding parameters
`x` where `A` is a symmetric positive definite matrix of small size `n×n` and
`b` a vector of length `n`.  For instance, `A` and `b` are the left hand-side
matrix and right hand-side vector of a linear least squares problem.

Only the lower-triangular part of `A` is specified.

"""
function solve1(A11::T, b1::T) where {T <: AbstractFloat}
    x1 = b1/A11
    return x1*b1, (x1,)
end

function solve2(A11::T, A21::T, A22::T,
                b1::T, b2::T) where {T <: AbstractFloat}
    A = SMatrix{2,2,T,4}(A11, A21,
                         A21, A22)
    b = SVector{2,T}(b1, b2)
    x = A\b
    return (x[1]*b1 + x[2]*b2), x
end

@doc @doc(solve1) solve2

function solve3(A11::T, A21::T, A31::T, A22::T, A32::T, A33::T,
                b1::T, b2::T, b3::T) where {T <: AbstractFloat}
    A = SMatrix{3,3,T,9}(A11, A21, A31,
                         A21, A22, A32,
                         A31, A32, A33)
    b = SVector{3,T}(b1, b2, b3)
    x = A\b
    return (x[1]*b1 + x[2]*b2 + x[3]*b3), x
end

@doc @doc(solve1) solve3

end # module


text:
#=
Created on 09/01/2021 14:50:19
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Implement Stirling's approximation for the factorial.
=#

using Plots

stirling(n) = n * log(n / exp(1)) + 0.5log(2π*n)

# here we use `m=:circle` to use the :circle symbols as a marker on our plots.

plot(1:10, n->log(factorial(big(n))), m=:circle, label="log-factorial")
xlabel!("n")
plot!(stirling, 1:10, m=:circle, label = "Stirling's approximation")


text:
using Indicomb
using Test, JSON3

@testset "Indico main function" begin
    t = get_events_catnum_name("https://indico.cern.ch", 6725, r"."; from="2021-06-01", to="2021-06-15")
    @test length(t) == 14
    t = get_events_catnum_name("https://indico.cern.ch", 6725, "Edition"; from="2021-06-01", to="2021-06-15")
    @test length(t) == 1
    @test contains(t[1][:title], "Edition of the Large")
end

@testset "Indico utility functions" begin
    re = Indicomb.indico_request("cern.ch/"; p = 3)
    @test re == "cern.ch/?p=3"

    t = Indicomb.get_indico_page(
        "https://indico.cern.ch",
        "/export/categ/6725.json";
        from="2021-06-01",
        to="2021-06-15",
    )
    @test t[:count] == 14
    @test t[:results][1][:_type] == "Conference"
    @test contains(t[:results][1][:title], "Edition of the Large")
end


text:
# nuclear Norm (times a constant)

export NuclearNorm

"""
**Nuclear norm**

    NuclearNorm(λ=1.0)

Returns the function
```math
f(X) = \\|X\\|_* = λ ∑_i σ_i(X),
```
where `λ` is a positive parameter and ``σ_i(X)`` is ``i``-th singular value of matrix ``X``.
"""
struct NuclearNorm{R <: Real} <: ProximableFunction
  lambda::R
  function NuclearNorm{R}(lambda::R) where {R <: Real}
    if lambda < 0
      error("parameter λ must be nonnegative")
    else
      new(lambda)
    end
  end
end

is_convex(f::NuclearNorm) = true

NuclearNorm(lambda::R=1.0) where {R <: Real} = NuclearNorm{R}(lambda)

function (f::NuclearNorm{R})(X::AbstractMatrix{T}) where {R <: Real, T <: Union{R, Complex{R}}}
  F = svd(X);
  return f.lambda * sum(F.S);
end

function prox!(Y::AbstractMatrix{T}, f::NuclearNorm{R}, X::AbstractMatrix{T}, gamma::R=one(R)) where {R <: Real, T <: Union{R, Complex{R}}}
  F = svd(X)
  S_thresh = max.(zero(R), F.S .- f.lambda*gamma)
  rankY = findfirst(S_thresh .== zero(R))
  if rankY === nothing
    rankY = minimum(size(X))
  end
  Vt_thresh = view(F.Vt, 1:rankY, :)
  U_thresh = view(F.U, :, 1:rankY)
  # TODO: the order of the following matrix products should depend on the shape of x
  M = S_thresh[1:rankY] .* Vt_thresh
  mul!(Y, U_thresh, M)
  return f.lambda * sum(S_thresh);
end

fun_name(f::NuclearNorm) = "nuclear norm"
fun_dom(f::NuclearNorm) = "AbstractArray{Real,2}, AbstractArray{Complex,2}"
fun_expr(f::NuclearNorm) = "X ↦ λ∑σ_i(X)"
fun_params(f::NuclearNorm) = "λ = $(f.lambda)"

function prox_naive(f::NuclearNorm, X::AbstractMatrix{T}, gamma=1.0) where T
  F = svd(X)
  S = max.(0, F.S .- f.lambda*gamma)
  Y = F.U * (Diagonal(S) * F.Vt)
  return Y, f.lambda * sum(S)
end


text:
using Meshes
using Tables
using Query
using Distances
using Statistics
using LinearAlgebra
using CategoricalArrays
using CircularArrays
using StaticArrays
using SparseArrays
using PlyIO
using Test, Random, Plots
using ReferenceTests, ImageIO

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

# helper function to read *.line files containing polygons
# generated with RPG (https://github.com/cgalab/genpoly-rpg)
function readpoly(T, fname)
  open(fname, "r") do f
    # read outer chain
    n = parse(Int, readline(f))
    outer = map(1:n) do _
      coords = readline(f)
      x, y = parse.(T, split(coords))
      Point(x, y)
    end

    # read inner chains
    inners = []
    while !eof(f)
      n = parse(Int, readline(f))
      inner = map(1:n) do _
        coords = readline(f)
        x, y = parse.(T, split(coords))
        Point(x, y)
      end
      push!(inners, inner)
    end

    # return polygonal area
    PolyArea(outer, inners)
  end
end

# helper function to read *.ply files containing meshes
function readply(T, fname)
  ply = load_ply(fname)
  x = ply["vertex"]["x"]
  y = ply["vertex"]["y"]
  z = ply["vertex"]["z"]
  points = Point{3,T}.(x, y, z)
  connec = [connect(Tuple(c.+1)) for c in ply["face"]["vertex_indices"]]
  SimpleMesh(points, connec)
end

include("dummy.jl")

# list of tests
testfiles = [
  "points.jl",
  "angles.jl",
  "polytopes.jl",
  "primitives.jl",
  "multi.jl",
  "collections.jl",
  "connectivities.jl",
  "topologies.jl",
  "toporelations.jl",
  "mesh.jl",
  "meshdata.jl",
  "traits.jl",
  "paths.jl",
  "distances.jl",
  "neighborhoods.jl",
  "neighborsearch.jl",
  "supportfun.jl",
  "diffops.jl",
  "views.jl",
  "viewing.jl",
  "sampling.jl",
  "partitioning.jl",
  "intersections.jl",
  "discretization.jl",
  "simplification.jl",
  "refinement.jl",
  "smoothing.jl",
  "boundingboxes.jl",
  "hulls.jl",
  "utils.jl"
]

# --------------------------------
# RUN TESTS WITH SINGLE PRECISION
# --------------------------------
T = Float32
P1, P2, P3 = Point{1,T}, Point{2,T}, Point{3,T}
V1, V2, V3 = Vec{1,T}, Vec{2,T}, Vec{3,T}
@testset "Meshes.jl ($T)" begin
  for testfile in testfiles
    println("Testing $testfile...")
    include(testfile)
  end
end

# --------------------------------
# RUN TESTS WITH DOUBLE PRECISION
# --------------------------------
T = Float64
P1, P2, P3 = Point{1,T}, Point{2,T}, Point{3,T}
V1, V2, V3 = Vec{1,T}, Vec{2,T}, Vec{3,T}
@testset "Meshes.jl ($T)" begin
  for testfile in testfiles
    println("Testing $testfile...")
    include(testfile)
  end
end


text:
# This file is part of the StateMachines.jl Julia package
#
# Paul Bayer
# MIT License
# Github: https://github.com/pbayer/StateMachines.jl
#

"A simple state has no internal vertices or transitions."
abstract type SimpleState <: AbstractState end

"A composite State contains at least one region."
abstract type CompositeState <: AbstractState end

"An orthogonal state is a composite state containing multiple regions."
abstract type OrthogonalState <: CompositeState end

"A submachine State refers to an entire state machine nested within."
abstract type SubmachineState <: AbstractState end

"""
    PseudoState

A Pseudostate is an abstraction for different transient states. A state machine
enters a pseudo state temporarily within a transition to the next stable state.
"""
abstract type PseudoState <: AbstractState end

# define typical simple states

struct Undefined <: SimpleState end
struct Idle <: SimpleState end
struct Setup <: SimpleState end
struct Busy <: SimpleState end
struct Blocked <: SimpleState end
struct Halted <: SimpleState end
struct Empty <: SimpleState end
struct Ready <: SimpleState end
struct Full <: SimpleState end
struct Failed <: SimpleState end
struct Waiting <: SimpleState end
struct InProcess <: SimpleState end

"""
    enter!(A::StateMachine, q::AbstractState)

A default function for entering a state if no more specific `enter!`-function
is defined. If entering a new state it calls the `exit`-action for the current state,
then it calls the `entry`- and `doActivity`-actions for the new state. If the new
state equals the current state, only the `doActivity`-action is executed.
"""
function enter!(A::StateMachine, q::AbstractState)
    if A.state == q
        doActivity(A, A.state)
    else
        exit(A, A.state)
        A.state = q
        entry(A, A.state)
        doActivity(A, A.state)
    end
end

"""
    entry(A::StateMachine, ::AbstractState)

A fallback entry function, doing nothing. The user can implement more specific
`entry`-actions, which are called by a concrete state machine upon entering a
concrete state.
"""
entry(A::StateMachine, ::AbstractState) = nothing

"""
    exit(A::StateMachine, ::AbstractState)

A fallback exit function, doing nothing. The user can implement more specific
`exit`-actions, which are called by a concrete state machine upon exiting a
concrete state.

"""
exit(A::StateMachine, ::AbstractState) = nothing

"""
    doActivity(A::StateMachine, ::AbstractState)

A fallback doActivity function, doing nothing. The user can implement more specific
`doActivity`-actions, which are called by a concrete state machine upon calling a
concrete state.

"""
doActivity(A::StateMachine, ::AbstractState) = nothing


text:
using Distributions

# Definitions of simulated environments
# and interaction of learners with the environment

# For the Categorical task (Section 2.3.2)
mutable struct EnvMultinomial
    ns::Int
    state::Int
    stochasticity::Float64
    changeprobability::Float64
    trans_probs::Array{Float64, 1}
    switchflag::Bool
    seed::Any
    rng::MersenneTwister # Used only for switches!
    nsteps::Int
end
function EnvMultinomial(; ns = 10, state = 1, stochasticity = 0.1,
                        changeprobability = .001, seed = 3,
                        nsteps = 10^5)
    rng = MersenneTwister(seed)
    trans_probs = rand(ENV_RNG, Dirichlet(ns, stochasticity))
    EnvMultinomial(ns, state, stochasticity, changeprobability,
                trans_probs, false, seed, rng, nsteps)
end
export EnvMultinomial

function interact!(env::EnvMultinomial)
    env.switchflag = false
    r = rand(env.rng)
     # @show r
    if r < env.changeprobability # Switch or not!
        #println("Switch!")
        env.trans_probs = rand(env.rng, Dirichlet(env.ns, env.stochasticity))
        env.switchflag = true
    end
    env.state = wsample(ENV_RNG, env.trans_probs)
end

# For the Gaussian task (Section 2.3.1)
mutable struct EnvGaussian
    mu_0::Float64
    sigma_0::Float64 # std of prior
    sigma::Float64 # std of observations
    mu_n::Float64
    state::Float64
    changeprobability::Float64
    switchflag::Bool
    seed::Any
    rng::MersenneTwister # Used only for switches!
    nsteps::Int
end
function EnvGaussian(; mu_0 = 0., sigma_0 = 1., sigma = 1.,
                        changeprobability = .001,
                        seed = 3,
                        nsteps = 10^5)
    rng = MersenneTwister(seed)
    mu_n = rand(ENV_RNG, Normal(mu_0, sigma_0))
    state = rand(ENV_RNG, Normal(mu_n, sigma))
    EnvGaussian(mu_0, sigma_0, sigma, mu_n, state, changeprobability,
                false, seed, rng, nsteps)
end
export EnvGaussian

function interact!(env::EnvGaussian)
    env.switchflag = false
    r = rand(env.rng)
     # @show r
    if r < env.changeprobability # Switch or not!
        # println("Switch!")
        env.mu_n = rand(env.rng, Normal(env.mu_0, env.sigma_0))
        env.switchflag = true
    end
    env.state = rand(ENV_RNG, Normal(env.mu_n, env.sigma))
end


text:
type Shot3
     isz :: Int64
     isx :: Int64
     irz :: Array{Int64, 1}
     irx :: Array{Int64, 1}
     ot  :: Float64
     nt  :: Int64
     dt  :: Float64
     vz  :: Array{Float64, 2}
     vx  :: Array{Float64, 2}
     p   :: Array{Float64, 2}
end

function InitShot3(isz::Int64, isx::Int64, pos::Array{Int64,2}, ot::Float64, nt::Int64, dt::Float64)
    nrec = size(pos, 1)
    irz  = pos[:, 1]
    irx  = pos[:, 2]
    d    = zeros(nt, nrec)
    shot = Shot3(isz, isx, irecz, irecx, ot, nt, dt, d, d, d)
    return shot
end

"""
copy(s3)

copy three components shot gather

"""
function Base.copy(s::Shot3)
    s1 = Shot3(s.isz, s.isx, s.irz, s.irx, s.ot, s.nt, s.dt, s.vz, s.vx, s.p)
    return s1
end

type Shot
     isz :: Int64
     isx :: Int64
     irz :: Array{Int64, 1}
     irx :: Array{Int64, 1}
     ot  :: Float64
     nt  :: Int64
     dt  :: Float64
     d   :: Array{Float64, 2}
end

function InitShot(isz::Int64, isx::Int64, pos::Array{Int64,2}, ot::Float64, nt::Int64, dt::Float64)
    nrec = size(pos, 1)
    irecz= pos[:, 1]
    irecx= pos[:, 2]
    d    = zeros(nt, nrec)
    shot = Shot(isz, isx, irecz, irecx, ot, nt, dt, d)
    return shot
end

function spt2shot3!(shot::Shot, spt::SnapShot, Nz::Int64, Nx::Int64, ext::Int64, zupper::Int64)
    nrec = length(shot.irx)
    it = spt.it
    for irec = 1 : nrec
        iz = shot.irz[irec] + zupper
        ix = shot.irx[irec] + ext
        ind= (ix-1)*Nz + iz
        shot.d[it, irec] = spt.pz[ind] + spt.px[ind]
    end
    return nothing
end

function AddShot32SnapShot!(spt::SnapShot, shot::Shot, Nz::Int64, Nx::Int64, ext::Int64, zupper::Int64)
    nrec = length(shot.irx)
    it = spt.it
    for irec = 1 : nrec
        iz = shot.irz[irec] + zupper
        ix = shot.irx[irec] + ext
        ind= (ix-1)*Nz + iz
        spt.pz[ind] = spt.pz[ind] + shot.d[it, irec]
        spt.px[ind] = spt.px[ind] + shot.d[it, irec]
    end
end

import Base.copy
function copy(s::Shot)
    s1 = Shot(s.isz, s.isx, s.irz, s.irx, s.ot, s.nt, s.dt, s.d)
    return s1
end

function spt2shot!(shot::Shot, spt::SnapShot, Nz::Int64, Nx::Int64, ext::Int64, zupper::Int64)
    nrec = length(shot.irx)
    it = spt.it
    for irec = 1 : nrec
        iz = shot.irz[irec] + zupper
        ix = shot.irx[irec] + ext
        ind= (ix-1)*Nz + iz
        shot.d[it, irec] = spt.pz[ind] + spt.px[ind]
    end
    return nothing
end

function AddShot2SnapShot!(spt::SnapShot, shot::Shot, Nz::Int64, Nx::Int64, ext::Int64, zupper::Int64)
    nrec = length(shot.irx)
    it = spt.it
    for irec = 1 : nrec
        iz = shot.irz[irec] + zupper
        ix = shot.irx[irec] + ext
        ind= (ix-1)*Nz + iz
        spt.pz[ind] = spt.pz[ind] + shot.d[it, irec]
        spt.px[ind] = spt.px[ind] + shot.d[it, irec]
    end
end

function wfd2shot(pos::Array{Int64,2}, path::ASCIIString)
    (isz, isx, nz, nx, ext, iflag, dt, ot, nt) = InfoWfd(path)
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    shot= InitShot(isz, isx, pos, ot, nt, dt)
    nrec= size(pos, 1)
    fid = open(path, "r"); seek(fid, sizeof(Int64)*8);
    for it = 1 : nt
        d = reshape(read(fid, Float64, Nz*Nx), Nz, Nx)
        for irec = 1 : nrec
            iz = shot.irz[irec] + zupper;
            ix = shot.irx[irec] + ext;
            shot.d[it, irec] = d[iz, ix]
        end
    end
    return shot
end

function fullWfd2Shot(pos::Array{Int64,2}, path::ASCIIString)
    (isz, isx, nz, nx, ext, iflag, ot, dt, nt) = InfoFullWfd(path)
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    shot= InitShot(isz, isx, pos, ot, nt, dt)
    nrec= size(pos, 1)
    fid = open(path, "r")
    head = sizeof(Float64) * (8+Nz*Nx*2)
    sizeSpt = sizeof(Float64)*Nz*Nx*4
    for it = 1 : nt
        position = head + sizeSpt*(it-1); seek(fid, position);
        pz = reshape(read(fid, Float64, Nz*Nx), Nz, Nx)
        px = reshape(read(fid, Float64, Nz*Nx), Nz, Nx)
        p  = pz + px
        for irec = 1 : nrec
            iz = shot.irz[irec] + zupper;
            ix = shot.irx[irec] + ext;
            shot.d[it, irec] = p[iz, ix]
        end
    end
    close(fid)
    return shot
end

function writeShot!(path::ASCIIString, s::Shot)
    fid = open(path, "w")
    write(fid, s.isz); write(fid, s.isx);
    nrec = length(s.irz); write(fid, nrec);
    write(fid, s.irz); write(fid, s.irx);
    write(fid, s.ot); write(fid, s.nt); write(fid, s.dt);
    write(fid, vec(s.d))
    close(fid)
end

function readShot(path::ASCIIString)
    fid = open(path, "r")
    isz = read(fid, Int64); isx = read(fid, Int64);
    nr  = read(fid, Int64)
    irz = read(fid, Int64, nr); irx = read(fid, Int64, nr);
    ot  = read(fid, Float64); nt = read(fid, Int64); dt = read(fid, Float64);
    d   = reshape(read(fid, Float64, nt*nr), nt, nr)
    s   = Shot(isz, isx, irz, irx, ot, nt, dt, d)
    return s
end

import Base.-
function -(shot1::Shot, shot2::Shot)
    if shot1.irz != shot2.irz || shot1.irx != shot2.irx || shot1.nt != shot2.nt
       error("size dismatch")
    end
    pos_rec = hcat(shot1.irz, shot1.irx)
    shot = InitShot(0, 0, pos_rec, 0.0, shot1.nt, shot1.dt)
    shot.d = shot1.d - shot2.d
    return shot
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
function optimize(obj::BaseOptimizer, arg0::Vector{OptimizationData})
    return jcall(obj, "optimize", Object, (Vector{OptimizationData},), arg0)
end


text:
using BerkeleyABC
using Test, Glob

abc_cmd = start_abc()

# for now use Nangate45_typ.lib for testing with a library
abc_cmd("read_lib libs/Nangate45_typ.lib")

@testset "ABC.jl" begin
    @testset "readers" begin
        @testset "aiger" begin
            for fn in Glob.glob("./designs/*.aiger")
                @test abc_cmd("read_aiger " * fn)[1] == 0
                @test abc_cmd("&get; &st; &dch; &nf; &put; ps; stime -p")[1] == 0
            end
        end
        @testset "blif" begin
            for fn in Glob.glob("./designs/*.blif")
                @test abc_cmd("read_blif " * fn)[1] == 0
                @test abc_cmd("&get; &st; &dch; &nf; &put; ps; stime -p")[1] == 0
            end
        end
    end
end


text:
module Backends # begin submodule PredictMD.Backends

import ..AbstractBackend

struct ClassImbalance <: AbstractBackend
end

struct DecisionTree <: AbstractBackend
end

struct Documenter <: AbstractBackend
end

struct Flux <: AbstractBackend
end

struct GLM <: AbstractBackend
end

struct Knet <: AbstractBackend
end

struct LIBSVM <: AbstractBackend
end

struct PGFPlots <: AbstractBackend
end

struct PGFPlotsX <: AbstractBackend
end

end # end submodule PredictMD.Backends


text:
stabilities(sweep::Codim1Sweep) = isstable.((sweep.timekind,), sweep.eigvals)

regions(xs::AbstractVector{<: Bool}) = regions(identity, xs)

function regions(f, xs::AbstractVector)
    ranges = UnitRange{Int}[]
    if isempty(xs)
        return ranges
    end

    s = 1
    b = f(xs[1])
    @inbounds for i in 2:length(xs)
        c = f(xs[i])
        if c != b
            push!(ranges, s:i-1)
            s = i
            b = c
        end
    end
    if s !== length(xs)
        push!(ranges, s:length(xs))
    end
    return ranges
end

# TODO: Merge curves_by_stability to curves(...; by=:stability)
function curves_by_stability(sweep::Codim1Sweep, vars)
    ss = stabilities(sweep)
    ranges = regions(ss)
    info = [ss[first(r)] for r in ranges]

    super = as(sweep, ContinuationSweep)
    TE = eltype(eltype(super))
    curves = Tuple(Vector{TE}[] for _ in vars)
    for r in ranges
        if first(r) == 1
            ran = r
        else
            ran = first(r) - 1:last(r)
        end
        for (i, v) in enumerate(vars)
            push!(curves[i], [u[v] for u in @view super.u[ran]])
        end
    end

    return info, curves
end

function curves_by_stability(sol::Codim1Solution)
    if isempty(sol.sweeps)
        return ([], [])
    end

    info, curves = curves_by_stability(sol.sweeps[1])
    for sweep in sol.sweeps[2:end]
        info2, curves2 = curves_by_stability(sweep)
        append!(info, info2)
        append!(curves, curves2)
    end
    return info, curves
end


text:
using Pkg
Pkg.activate(".")

using Genie, HTTP
import Genie.Router: route, POST, @params
import Genie.Requests: jsonpayload

route("/jsonpayload", method = POST) do
  @show @params(:JSON_PAYLOAD)
end

route("/jsontest", method = POST) do
    json_request = jsonpayload()
    json_request["test"]
end

Genie.AppServer.startup(; open_browser = false)

HTTP.request("POST", "http://localhost:8000/jsonpayload", [("Content-Type", "application/json; charset=utf-8")], """{"greeting":"hello"}""")
HTTP.request("POST", "http://localhost:8000/jsonpayload", [("Content-Type", "application/json")], """{"greeting":"hello"}""")
HTTP.request("POST", "http://localhost:8000/jsontest", [("Content-Type", "application/json")], """{"test": [1,2,3]}""")

exit(0)


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ b129ba7c-953a-11ea-3379-17adae34924c
md"# _Welcome to Pluto!_

Pluto is a programming environment for _Julia_, designed to be **interactive** and **helpful**.

In this introduction, we will go through the basics of using Pluto. To make it interesting, this notebook does something special: it **changes while you work on it**. Computer magic ✨"

# ╔═╡ 4d88b926-9543-11ea-293a-1379b1b5ae64
md"## Cats
Let's say you're like my grandma, and you have a lot of cats. Our story will be about them."

# ╔═╡ 3e8e381e-953f-11ea-3d3e-71d0fea52560
cat = "Kitty"

# ╔═╡ aeb3a6bc-9540-11ea-0b8f-6d37412bfe68
if cat == "Ks*;lj"
	md"Oh no! Someone messed with my pretty introduction. Change the code above to give our cat a proper name!"
else
	HTML("""
Well done, your cat is called $cat now. This text gets updated every time you change the name. To see how the magic works, click on the  to the left of this text.
""")
end

# ╔═╡ 611c28fa-9542-11ea-1751-fbdedcfb7690
html"""
To edit any code, just click on it. When you're done programming, press the  in the lower-right corner of a cell to run the code. You can also use `Shift+Enter` if you are in a hurry.
"""

# ╔═╡ d6a563f1-8676-4274-86d0-ad950abb2438
friend = "Meow"

# ╔═╡ 6f7eecec-9543-11ea-1284-dd52fce3ecca
md"I feel like our cat needs a friend. Let's call them $(friend)."

# ╔═╡ a1a20314-9543-11ea-17de-0b658da18992
if !@isdefined friend
	md"Uh oh, what is this? I forgot to add a cell defining our friend. Can you do it for me?"
else
	md"**Well done!** $cat and $friend are both happy with your performance."
end

# ╔═╡ f112b662-9543-11ea-3dcb-2906a99b2188
html"""
A cell is a container for code & output. To add one, click on the  above or below another cell. You can do it wherever you like. After you're done writing code in your cell, remember to run it!
"""

# ╔═╡ e0642f42-9545-11ea-14ee-fde52cb54ccc
md"## Feeding neighbors
Our cats have some neighbors. Let's involve them in the story too!"

# ╔═╡ 3653b1ac-9546-11ea-2a44-ddd3054636fe
neighbors = [cat, friend, "Smerfetka", "Latte"]

# ╔═╡ 19ff8d36-9547-11ea-0e08-e5cdd8338673
md"Now, if you're like my grandma, you're feeding the entire neighborhood by yourself. Let's see how many cans of cat food you need to prepare."

# ╔═╡ 270ac49e-9549-11ea-3ffd-71ddaee9f134
md"But what does `confusing_function` do? If you ever need help, click on 📚 **Live docs** in the lower right, and then place your cursor on the code you need help with.

If you don't see it, then your screen is too small! Maybe you need to zoom out?"

# ╔═╡ 745a4584-954a-11ea-028e-59011f268ec6
cans_in_stock = "🥫🥫🥫🥫"

# ╔═╡ 55ad7422-954e-11ea-0a33-a3b03febb56e
if @isdefined cans_in_stock
	md"Actually, I have a hunch there will be another cat coming. Uncomment the code below (remove the #) to add one more can. Remember to run it after making the change!"
else
	md"**Whoopsie!** Because Pluto figures out execution order for you, it doesn't really make sense to assign to the same variable twice. A smarter way to plan ahead is to write `cans_in_stock = consumption` — Pluto will take care of updating everything."
end

# ╔═╡ eac62fea-954e-11ea-2768-39ce6f4059ab
# cans_in_stock = "🥫🥫🥫🥫🥫"

# ╔═╡ 6c8e2108-9550-11ea-014d-235770ed4771
md"## Saving cats and notebooks
Alright, we have a neighborhood full of well-fed cats. But oh no, here comes..."

# ╔═╡ 9e89fc9a-9550-11ea-14b4-7f0e96225ec0
scary_dog = "Piesio"

# ╔═╡ bdd5d268-9550-11ea-1046-31efedc36872
if @isdefined scary_dog
	md" $scary_dog is terrorizing the neighborhood! We must do something about it!"
else
	md"You saved the neighborhood! Referencing `scary_dog` leads to an `UndefVarError`, as if it never even existed."
end

# ╔═╡ 36cd006a-9551-11ea-3c0c-df8b7f2843c1
HTML("""
To delete a cell like the one defining $scary_dog, click on the  on the right of its code.
""")

# ╔═╡ fb4e471c-9551-11ea-1ab5-41bbd5de76b8
md"""
Speaking of saving, this notebook is autosaved whenever you change something. The default location for new notebooks is $(joinpath(first(DEPOT_PATH), "pluto_notebooks")), you can find it using your file explorer. To change the name or the directory of a notebook, scroll to the top - you enter the notebook's path next to the Pluto logo.
"""

# ╔═╡ 9d3af596-9554-11ea-21bd-bf427c91c424
md"## ⚡ Pluto power ⚡
Oof, that dog situation was intense. Let's give our cats superpowers to make sure it never happens again!"

# ╔═╡ 3150bf1a-9555-11ea-306f-0fd4d9229a51
md"Remember learning HTML in junior high? Now you can use it for something! Pluto lets you `@bind` variables to HTML elements. As always, every time you change something, Pluto knows what to update!"

# ╔═╡ f2c79746-9554-11ea-39ca-298fd09248ad
@bind power_level html"<input type='range'>"

# ╔═╡ 0b094cea-9556-11ea-268e-0d270fd04d56
md"The power level is $power_level, but we should do more than just say it - let's equip our cats with $power_level emoji!"

# ╔═╡ 1908f9f2-9557-11ea-2abd-dd52f8d776f4
power_emoji = "⚡"
power = repeat(power_emoji, power_level)

# ╔═╡ 784b1774-9557-11ea-315e-d1ea277ce0fd
if !@isdefined power
	md"Uh oh! Pluto doesn't support multiple expressions per cell. This is a conscious choice - this restriction helps you write less buggy code once you get used to it. To fix the code, you can either split the above cell, or wrap it in a `begin ... end` block. Both work."
else
	md"**Well done!** Your cats have powers now."
end

# ╔═╡ 5edadcd2-9554-11ea-1714-b5b7692c4797
html"""
We're almost done! It's time to share your amazing story. Scroll to the top of the notebook, and click on  to see the export options - or you can always share this notebook's save file. (The file is pure Julia, by the way, and it's runnable! You'll learn more about this in the advanced introduction.)
"""

# ╔═╡ 4634c856-9553-11ea-008d-3539195970ea
md"## Final notes"

# ╔═╡ 4d0ebb46-9553-11ea-3431-2d203f594815
md"If anything about this introduction left you confused, something doesn't work, or you have a cool new idea - don't hesitate to contact us! You can do it right from this screen, using the `Instant feedback` form in the bottom right."

# ╔═╡ d736e096-9553-11ea-3ba5-277afde1afe7
md"Also, if you were wondering where `confusing_function` came from, here you go! Remember that you, too, can place code wherever you like."

# ╔═╡ 7366f1b6-954c-11ea-3f68-b12444c902c3
"""
confusing_function(text::String, array::Array)

Repeats the `text` as many times as there are elements in `array`.
"""
confusing_function(text::String, array::Array) = repeat(text, length(array))

# ╔═╡ a4a60262-9547-11ea-3a81-5bf7f9ee5d16
consumption = confusing_function("🥫", neighbors)

# ╔═╡ e11e1660-9549-11ea-22f6-8bb53dc045fe
md"Now we know to prepare $(length(consumption)) cans. Let's stock up!"

# ╔═╡ f27f90c2-954f-11ea-3f93-17acb2ce4280
md"We have $(length(cans_in_stock)) cans of cat food, and our cats need $(length(consumption)). Try adding another cat to the neighborhood to see what changes!"

# ╔═╡ 1deaaf36-9554-11ea-3dae-85851f73dbc6
md"**Have fun using Pluto!**

_~ Fons van der Plas & Nicholas Bochenski_"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.9.0-DEV"
manifest_format = "2.0"
project_hash = "da39a3ee5e6b4b0d3255bfef95601890afd80709"

[deps]
"""

# ╔═╡ Cell order:
# ╟─b129ba7c-953a-11ea-3379-17adae34924c
# ╟─4d88b926-9543-11ea-293a-1379b1b5ae64
# ╠═3e8e381e-953f-11ea-3d3e-71d0fea52560
# ╟─aeb3a6bc-9540-11ea-0b8f-6d37412bfe68
# ╟─611c28fa-9542-11ea-1751-fbdedcfb7690
# ╠═d6a563f1-8676-4274-86d0-ad950abb2438
# ╠═6f7eecec-9543-11ea-1284-dd52fce3ecca
# ╟─a1a20314-9543-11ea-17de-0b658da18992
# ╟─f112b662-9543-11ea-3dcb-2906a99b2188
# ╟─e0642f42-9545-11ea-14ee-fde52cb54ccc
# ╠═3653b1ac-9546-11ea-2a44-ddd3054636fe
# ╟─19ff8d36-9547-11ea-0e08-e5cdd8338673
# ╠═a4a60262-9547-11ea-3a81-5bf7f9ee5d16
# ╟─270ac49e-9549-11ea-3ffd-71ddaee9f134
# ╟─e11e1660-9549-11ea-22f6-8bb53dc045fe
# ╠═745a4584-954a-11ea-028e-59011f268ec6
# ╟─55ad7422-954e-11ea-0a33-a3b03febb56e
# ╠═eac62fea-954e-11ea-2768-39ce6f4059ab
# ╟─f27f90c2-954f-11ea-3f93-17acb2ce4280
# ╟─6c8e2108-9550-11ea-014d-235770ed4771
# ╠═9e89fc9a-9550-11ea-14b4-7f0e96225ec0
# ╟─bdd5d268-9550-11ea-1046-31efedc36872
# ╟─36cd006a-9551-11ea-3c0c-df8b7f2843c1
# ╟─fb4e471c-9551-11ea-1ab5-41bbd5de76b8
# ╟─9d3af596-9554-11ea-21bd-bf427c91c424
# ╟─3150bf1a-9555-11ea-306f-0fd4d9229a51
# ╠═f2c79746-9554-11ea-39ca-298fd09248ad
# ╟─0b094cea-9556-11ea-268e-0d270fd04d56
# ╠═1908f9f2-9557-11ea-2abd-dd52f8d776f4
# ╟─784b1774-9557-11ea-315e-d1ea277ce0fd
# ╟─5edadcd2-9554-11ea-1714-b5b7692c4797
# ╟─4634c856-9553-11ea-008d-3539195970ea
# ╟─4d0ebb46-9553-11ea-3431-2d203f594815
# ╟─d736e096-9553-11ea-3ba5-277afde1afe7
# ╟─7366f1b6-954c-11ea-3f68-b12444c902c3
# ╟─1deaaf36-9554-11ea-3dae-85851f73dbc6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
using DashHtmlComponents

function Example(example::String)
    return html_div(
        example,
        className = "example-container"
    )
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

const _default_PT = Float32 # Default type/eltype for variate/parameter values

abstract type AbstractVarBounds end

function Base.intersect(a::AbstractVarBounds, b::AbstractVarBounds)
    totalndof(a) != totalndof(b) && throw(ArgumentError("Can't intersect bounds with different number of DOF"))
    unsafe_intersect(a, b)
end

@enum BoundsType hard_bounds=1 cyclic_bounds=2 reflective_bounds=3

function Base.intersect(a::BoundsType, b::BoundsType)
    if a == hard_bounds || b == hard_bounds
        hard_bounds
    elseif a == reflective_bounds || b == reflective_bounds
        reflective_bounds
    else
        cyclic_bounds
    end
end

@inline float_iseven(n::T) where {T<:AbstractFloat} = (n - T(2) * floor((n + T(0.5)) * T(0.5))) < T(0.5)

@doc doc"""
    renormalize_variate!(v_renorm::Any, bounds::AbstractVarBounds, v::Any)

*BAT-internal, not part of stable public API.*

Apply renormalization (if any) inherent `bounds`, to variate/parameters `x`.
The renormalization must preserve phase-space volume locally (determinant
of Jacobian must be one).

Modifies and returns `v_renorm`.
"""
function renormalize_variate! end

@doc doc"""
    renormalize_variate_impl(x::<:Real, lo::<:Real, hi::<:Real, boundary_type::BoundsType)

*BAT-internal, not part of stable public API.*

Apply lower/upper bound `lo`/`hi` to value `x`. `boundary_type` may be
`hard_bounds`, `cyclic_bounds` or `reflective_bounds`.
"""
@inline function renormalize_variate_impl(x::X, lo::L, hi::H, boundary_type::BoundsType) where {X<:Real,L<:Real,H<:Real}
    T = float(promote_type(X, L, H))

    offs = ifelse(x < lo, lo - x, x - hi)
    hi_lo = hi - lo
    nwrapped = floor(offs / (hi_lo))
    even_nwrapped = float_iseven(nwrapped)
    wrapped_offs = muladd(-nwrapped, (hi_lo), offs)

    hb = (boundary_type == hard_bounds)
    rb = (boundary_type == reflective_bounds)

    ifelse(
        lo <= x <= hi,
        convert(T, x),
        ifelse(
            hb,
            convert(T, x),
            ifelse(
                (x < lo && (!rb || rb && !even_nwrapped)) || (x > lo && rb && even_nwrapped),
                convert(T, hi - wrapped_offs),
                convert(T, lo + wrapped_offs)
            )
        )
    )
end

@doc doc"""
    renormalize_variate_impl(x::Real, interval::ClosedInterval, boundary_type::BoundsType)

*BAT-internal, not part of stable public API.*

Specify lower and upper bound via `interval`.
"""
@inline renormalize_variate_impl(x::Real, interval::ClosedInterval, boundary_type::BoundsType) =
    renormalize_variate_impl(x, minimum(interval), maximum(interval), boundary_type)

struct NoVarBounds <: AbstractVarBounds
    ndims::Int
end

# TODO: Find better solution to determine value type if no bounds:
Base.eltype(bounds::NoVarBounds) = _default_PT

Base.in(x::Any, bounds::NoVarBounds) = true

Base.isinf(bounds::NoVarBounds) = true

ValueShapes.totalndof(b::NoVarBounds) = b.ndims

function renormalize_variate!(v_renorm::AbstractVector{<:Real}, bounds::NoVarBounds, v::AbstractVector{<:Real})
    if v_renorm !== v
        v_renorm .= v
    end
    v_renorm
end

unsafe_intersect(a::NoVarBounds, b::NoVarBounds) = a
unsafe_intersect(a::AbstractVarBounds, b::NoVarBounds) = a
unsafe_intersect(a::NoVarBounds, b::AbstractVarBounds) = b

abstract type VarVolumeBounds{T<:Real, V<:SpatialVolume{T}} <: AbstractVarBounds end

Base.in(x::Any, bounds::VarVolumeBounds) = in(x, bounds.vol)

# Random.rand(rng::AbstractRNG, bounds::VarVolumeBounds) =
#     rand!(rng, bounds, Vector{float(eltype(bounds))}(totalndof(bounds)))

# Random.rand(rng::AbstractRNG, bounds::VarVolumeBounds, n::Integer) =
#     rand!(rng, bounds, Matrix{float(eltype(bounds))}(npartotalndofams(bounds), n))
#
# Random.rand!(rng::AbstractRNG, bounds::VarVolumeBounds, x::StridedVecOrMat{<:Real}) = rand!(rng, spatialvolume(bounds), x)

ValueShapes.totalndof(b::VarVolumeBounds) = ndims(b.vol)

@doc doc"""
    spatialvolume(b::VarVolumeBounds)::SpatialVolume

*BAT-internal, not part of stable public API.*

Returns the spatial volume that defines the variate/parameter bounds.
"""
function spatialvolume end

struct HyperRectBounds{T<:Real} <: VarVolumeBounds{T, HyperRectVolume{T}}
    vol::HyperRectVolume{T}
    bt::Vector{BoundsType}

    function HyperRectBounds{T}(vol::HyperRectVolume{T}, bt::AbstractVector{BoundsType}) where {T<:Real}
        nd = ndims(vol)
        axes(bt) != (1:nd,) && throw(ArgumentError("bt must have indices (1:ndims(vol),)"))
        (nd > 0) && isempty(vol) && throw(ArgumentError("Cannot create bounds with emtpy volume of $nd dimensions"))
        new{T}(vol, bt)
    end
end

HyperRectBounds(vol::HyperRectVolume{T}, bt::AbstractVector{BoundsType}) where {T<:Real} = HyperRectBounds{T}(vol, bt)
HyperRectBounds(lo::AbstractVector{T}, hi::AbstractVector{T}, bt::AbstractVector{BoundsType}) where {T<:Real} = HyperRectBounds(HyperRectVolume(lo, hi), bt)
HyperRectBounds(lo::AbstractVector{T}, hi::AbstractVector{T}, bt::BoundsType) where {T<:Real} = HyperRectBounds(lo, hi, fill(bt, size(lo, 1)))
HyperRectBounds(intervals::AbstractVector{<:ClosedInterval{<:Real}}, bt) = HyperRectBounds(minimum.(intervals), maximum.(intervals), bt)
HyperRectBounds(bounds::AbstractInterval, bt::BoundsType) = HyperRectBounds([bounds.left], [bounds.right], bt)

function HyperRectBounds(bounds::Vector{<:AbstractInterval}, bt::BoundsType)
    lo = [b.left for b in bounds]
    hi = [b.right for b in bounds]

    return HyperRectBounds(lo, hi, bt)
end

Base.similar(bounds::HyperRectBounds) = HyperRectBounds(
        HyperRectVolume(
            fill!(similar(bounds.vol.lo), zero(eltype(bounds.vol.lo))),
            fill!(similar(bounds.vol.hi), one(eltype(bounds.vol.hi)))
        ),
        similar(bounds.bt)
    )

Base.eltype(bounds::HyperRectBounds{T}) where {T} = T

function Base.intersect(a::HyperRectBounds, b::HyperRectBounds)
    c = similar(a)
    for i in eachindex(a.vol.lo, a.vol.hi, a.bt, b.vol.lo, b.vol.hi, b.bt)
        iv_a = a.vol.lo[i]..a.vol.hi[i]
        iv_b = b.vol.lo[i]..b.vol.hi[i]
        if issubset(iv_a, iv_b)
            c.bt[i] = a.bt[i]
        elseif issubset(iv_b, iv_a)
            c.bt[i] = b.bt[i]
        else
            c.bt[i] = a.bt[i] ∩ b.bt[i]
        end

        iv_c = iv_a ∩ iv_b
        c.vol.lo[i] = minimum(iv_c)
        c.vol.hi[i] = maximum(iv_c)
    end
    c
end

Base.isinf(bounds::HyperRectBounds) = isinf(bounds.vol)

function Base.vcat(xs::HyperRectBounds...)
    lo = vcat(map(x -> x.vol.lo, xs)...)
    hi = vcat(map(x -> x.vol.hi, xs)...)
    bt = vcat(map(x -> x.bt, xs)...)
    HyperRectBounds(lo, hi, bt)
end

spatialvolume(bounds::HyperRectBounds) = bounds.vol

function renormalize_variate!(v_renorm::AbstractVector{<:Real}, bounds::HyperRectBounds, v::AbstractVector{<:Real})
    VT = typeof(v)
    v_renorm .= renormalize_variate_impl.(v, bounds.vol.lo, bounds.vol.hi, bounds.bt)
end


text:
"""
    House

Everything that makes a house a `House`.
"""
struct House
    nbedrooms::Int
    nbathrooms::Int
    rent::Float64
    area::Float64
end

"""
Describes the range of values we are dealing with.
"""
struct MarketData
    bedrooms_range::Vector{Int}
    bedrooms_frequency::Weights
    bathrooms_range::Vector{Int}
    bathrooms_frequency::Weights
    prices_range_pp::Vector{Float64}
    area_range::Vector{Float64}
end

"""
    maintenance(h::House)

Monthly cost we take on having invested in house `h`.
"""
maintenance(h::House) = h.rent / 10

"""
Average square footage of a house with a given number of bedrooms.
"""
function avg_area(nbedrooms::Int, nbathrooms::Int)
    nbedrooms * 300 + nbathrooms * 15
end
"""
Average rent we expect for a house with a given number of bedrooms, bathrooms,
    and area.
"""
function avg_rent(nbedrooms::Int, nbathrooms::Int, area::Float64)
    area * 1.4 + nbathrooms * 50 - nbedrooms * 10
end

"""
    Characteristic

Everything that anyone ever cared about when choosing where to live.
"""
struct Characteristic
    nbedrooms::Int
    nbathrooms::Int
    price_at_most::Int
    area_at_least::Int
end

"""
    house_fits_characteristic(house::House, characteristic::Characteristic)

Returns a boolean depending on whether `house` fits a description given by
`characteristic`.
"""
function house_fits_characteristic(house::House, characteristic::Characteristic, md::MarketData)
    house.nbedrooms == characteristic.nbedrooms &&
    house.nbathrooms == characteristic.nbathrooms &&
    house.rent <= md.prices_range_pp[characteristic.price_at_most] &&
    house.area >= md.area_range[characteristic.area_at_least]
end

"""
This will become our column generation part later.
"""
function get_all_characteristics(md::MarketData)
    all_possible_characteristics = Characteristic[]
    for nbed in md.bedrooms_range
        for nbath in md.bathrooms_range
            nbath > nbed && continue
            for p in 1:length(md.prices_range_pp)
                for a in 1:length(md.area_range)
                    push!(all_possible_characteristics, Characteristic(nbed, nbath, p, a))
                end
            end
        end
    end
    all_possible_characteristics
end

"""
    gethouses(nhouses::Int)

Get the houses that we expect will be available on the market. Since we have no
real data, we generate some synthetic ones here.
"""
function gethouses(nhouses::Int, md::MarketData)
    @assert minimum(md.bathrooms_range) <= minimum(md.bedrooms_range)
    houses = House[]
    for h = 1:nhouses
        nbedrooms = sample(md.bedrooms_range, md.bedrooms_frequency)
        nbathrooms = 1
        while true
            nbathrooms = sample(md.bathrooms_range, md.bathrooms_frequency)
            (nbathrooms <= nbedrooms) && break
        end
        area = avg_area(nbedrooms, nbathrooms) + rand(Normal(0.0, 10.0))
        rent = avg_rent(nbedrooms, nbathrooms, area) + rand(Normal(0.0, 100.0))
        push!(houses, House(nbedrooms, nbathrooms, rent, area))
    end
    houses
end

"""
    getdemand(npatterns::Int, nnoises::Int, nstages::Int, d::ContinuousDistribution)

Generates radnom demands for a problem where the number of periods is
`nstages` and we want `nnoises` possible realizations of demand each period.
All demands are sampled from `d`.
"""
function getdemand(npatterns::Int, nnoises::Int, nstages::Int, d::ContinuousDistribution)
    demands = Array{Float64,3}(npatterns, nnoises, nstages)
    for n in 1:nnoises
        demands[:, n, :] .= round.(rand(d, npatterns, nstages))
    end
    demands
end

"""
    function dominates(c1::Characteristic, c2::Characteristic)

Returns true if `c1` dominates `c2`.
"""
function dominates(c1::Characteristic, c2::Characteristic)
    (c1.nbedrooms <= c2.nbedrooms) && # fewer roommates
    (c1.nbathrooms >= c2.nbathrooms) &&
    (c1.price_at_most <= c2.price_at_most) &&
    (c1.area_at_least >= c2.area_at_least)
end

"""
    pattern_is_legal(p::Vector{Int}, all_characteristics::Vector{Characteristic})

Checks whether there are logical conflicts in the characteristics referenced
within `p`.
"""
function pattern_is_legal(p::Vector{Int}, all_characteristics::Vector{Characteristic})
    for i = 1:length(p)
        if p[i] == 1
            # check that it is not dominated by anything with a 0
            for j = 1:length(p)
                i == j && continue
                if p[j] == 0 && dominates(all_characteristics[j], all_characteristics[i])
                    return false
                end
            end
        end
    end
    # Didn't find any conflicts, so if not empty, then pattern is legal
    true
end

# function make_all_patterns(all_characteristics::Vector{Characteristic})
#     C = length(all_characteristics)
#     C > 63 && error("Too many characteristics.")
#     max_patterns = 2^C
#     # Divide patterns to process evenly among processors
#     all_patterns = SharedArray{Vector{Int}}(max_patterns)
#     procs = collect(1:nprocs())
#     breakpoints = [linspace(1, max_patterns, i) for i = 1:nprocs()]
#     all_patterns = Vector{Int}[]
#     p = zeros(Int, length(all_characteristics))
#     while true
#         addone!(p)
#         if pattern_is_legal(p, all_characteristics)
#             push!(all_patterns, copy(p))
#         end
#         if all(p .== 1)
#             break
#         end
#     end
#     all_patterns
# end

function make_all_patterns(all_characteristics::Vector{Characteristic})
    all_patterns = Vector{Int}[]
    p = zeros(Int, length(all_characteristics))
    while true
        addone!(p)
        if pattern_is_legal(p, all_characteristics)
            push!(all_patterns, copy(p))
        end
        if all(p .== 1)
            break
        end
    end
    all_patterns
end

function beds_needed(p::Int, easy::Bool=true)
    if easy
        ones(p)
    else
        round.(rand(Uniform(1.0, 3.0), p))
    end
end

struct StudentHousingData
    nstages::Int
    budget::Float64
    all_characteristics::Vector{Characteristic}
    houses::Vector{House}
    market_data::MarketData
    demands::Array{Float64,3}
    patterns_allow::Vector{Vector{Int}}
    beds_needed::Vector{Int}
end

function check_data(nhouses::Int, nnoises::Int)
    (nhouses > 0) && (nnoises > 0)
end

"""
    StudentHousingData(market_data; nstages::Int=1,
                            nhouses::Int=1,
                            nnoises::Int=1,
                            demand_distribution::ContinuousDistribution=Normal(100,10),
                            budget::Float64=1e6

Holds all data for the housing problem.
"""
function StudentHousingData(market_data; nstages::Int=1,
                            nhouses::Int=1,
                            nnoises::Int=1,
                            demand_distribution::ContinuousDistribution=Normal(100,10),
                            budget::Float64=1e6,
                            easy::Bool=true
                        )

    check_data(nhouses, nnoises) || error("Specified non-positive values for positive parameters.")
    # Get all possible combinations of market characteristics
    all_characteristics = get_all_characteristics(market_data)
    # Generate some data describing houses
    houses = gethouses(nhouses, market_data)

    # Isolate only the patterns that make logical sense without ever computing
    # all patterns
    patterns = make_all_patterns(all_characteristics)
    # Cache characteristics each pattern will allow
    npatterns = length(patterns)
    patterns_sparse = Vector{Vector{Int}}(npatterns)
    for i = 1:npatterns
        patterns_sparse[i] = find(patterns[i] .> 0)
    end
    bedsneeded = beds_needed(npatterns, easy)

    # The number of patterns is the number of logically sound patterns
    npatterns = length(patterns_sparse)
    # Generate demand data
    demands = getdemand(npatterns, nnoises, nstages, demand_distribution)
    StudentHousingData(nstages, budget, all_characteristics, houses,
        market_data, demands, patterns_sparse, bedsneeded)
end

# Some helper functions

beds_avail(i::Int, houses::Vector{House}) = houses[i].nbedrooms
maintenance(i::Int, houses::Vector{House}) = maintenance(houses[i])
get_ncharacteristics(d::StudentHousingData) = length(d.all_characteristics)
get_npatterns(d::StudentHousingData) = length(d.patterns_allow)
beds_needed(p::Int, d::StudentHousingData) = d.beds_needed[p]


text:
using Combinatorics
using SparseArrays
using TensorOperations
using LinearAlgebra
using ArnoldiMethod

function CASCI{T}(Alg::ACI) where T <: AbstractFloat
    # we need this implementation
    @output "Getting molecule...\n"
    #molecule = Molecule()
    @output "Computing AO Integrals...\n"
    #aoint = ConventionalAOIntegrals()

    @output "Calling RHF module...\n"
    refwfn = Fermi.HartreeFock.RHF()
    CASCI{T}(refwfn, Alg)
end

function CASCI{T}(refwfn::Fermi.HartreeFock.RHF, Alg::ACI; ci = nothing) where T <: AbstractFloat
    @output "Generating Integrals for CAS computation...\n"
    #aoint = ConventionalAOIntegrals()
    ints = refwfn.ints
    @output "Transforming Integrals for CAS computation...\n"
    frozen = Fermi.CurrentOptions["cas_frozen"]

    nmo = refwfn.ndocc + refwfn.nvir

    act_elec = 2*(refwfn.ndocc - frozen)

    if act_elec < 0
        error("\nInvalid number of frozen orbitals ($frozen) for $(2*refwfn.ndocc) electrons.")
    end

    # Active = -1 means FCI, with frozen
    if Fermi.CurrentOptions["cas_active"] == -1
        active = nmo - frozen
    else
        active = Fermi.CurrentOptions["cas_active"]
    end

    if active ≤ act_elec/2
        error("\nNumber of active orbitals ($active) too small for $(act_elec) active electrons")
    end

    if active+frozen > nmo
        error("\nNumber of active ($active) and frozen orbitals ($frozen) greater than number of orbitals ($nmo)")
    end

    s = 1:(frozen+active)

    h = T.(Fermi.Integrals.transform_fock(ints["T"] + ints["V"], ints.orbs["FU"][s], ints.orbs["FU"][s]))
    V = T.(Fermi.Integrals.transform_eri(ints["μ"], ints.orbs["FU"][s], ints.orbs["FU"][s], ints.orbs["FU"][s], ints.orbs["FU"][s]))

    aoint = nothing

    aoint = nothing
    CASCI{T}(refwfn, h, V, frozen, act_elec, active, Alg, ci=ci)
end

function CASCI{T}(refwfn::Fermi.HartreeFock.RHF, h::Array{T,2}, V::Array{T,4}, frozen::Int, act_elec::Int, active::Int, Alg::ACI; ci=nothing) where T <: AbstractFloat

    # Print intro
    Fermi.ConfigurationInteraction.print_header()
    ttotal = @elapsed begin
    @output "\n    • Computing FCI with the ACI algorithm.\n\n"
    act_range = (frozen+1):(active+frozen)
    σ = Fermi.CurrentOptions["σ"]
    γ = Fermi.CurrentOptions["γ"]
    pe = Fermi.CurrentOptions["aci_print_screen"]

    @output "\n →  ACTIVE SPACE\n"
    @output "Frozen Orbitals:  {:3d}\n" frozen
    @output "Active Electrons: {:3d}\n" act_elec
    @output "Active Orbitals:  {:3d}\n" active

    # Start reference space as HF
    zeroth = repeat('1', frozen)*repeat('1', Int(act_elec/2))
    if ci == nothing
        P = [Determinant(zeroth, zeroth)]
        Pcoef = [1.0]
    else
        P = deepcopy(ci.dets)
        Pcoef = deepcopy(ci.coef)
    end
    E = refwfn.energy - refwfn.molecule.Vnuc
    ΔE = 1.0
    ite = 1

    @output repeat("=",50)*"\n"
    Nα = Int(act_elec/2)
    Nβ = Int(act_elec/2)
    Lenny = length(P)
    M = nothing
    ϵsum = nothing
    ϵest = nothing
    oldP = nothing
    cflag = false
    while true
        @output " → Iteration {}\n\n" ite
        @output "   • P\n"
        @output "Initial model space size: {}\n\n" length(P)

        @output "   • P ⇒ F\n"
        @output "Generating First Order Interacting Space...\n"
        t = @elapsed F = get_fois(P, Nα, Nβ, act_range)
        @output "FOIS size:                {}\n" length(F)
        @output "FOIS contructed in {:5.5f} s.\n\n" t

        @output "   • F ⇒ Q\n"
        @output "Screening FOIS using 2-D Hamiltonian\n" σ
        t = @elapsed Fe = ϵI(F, P, Pcoef, E, h, V, act_elec, active)
        @output "Screen complete in {:5.5} s.\n" t
        @output "Sorting F space...\n"
        _t = @elapsed begin
            Fperm = zeros(Int, length(Fe))
            sortperm!(Fperm, Fe, by=abs)
            reverse!(Fperm)
            Fe = Fe[Fperm]
            F = F[Fperm]
        end
        @output "Sorted in {:5.5f} s.\n" _t
        @output "Filtering F..."
        t = @elapsed begin
        ϵest = 0.0
        ϵsum = 0.0
        while true
            if length(Fe) == 0
                #then no determinants were deemed important - exit ACI
                break
            end
            if ϵsum ≤ σ
                ϵest += Fe[end]
                ϵsum += abs(Fe[end])
                pop!(Fe)
                pop!(F)
            else
                ϵest -= Fe[end]
                ϵsum -= abs(Fe[end])
                break
            end
        end
        end
        @output " Secondary space (Q) built in {:5.5f}\n\n" t
        @output "Size of Q {}\n" length(Fe)

        @output "   • M = P ∪ Q\n"
        ΔE = -E
        M = vcat(P, F)
        @output "Model space size: {}\n" length(M)
        @output "Updating model space...\n" length(M)
        t = @elapsed E, Pcoef, P = update_model_space(M, h, V)
        @output " Model space updated in {:5.5f}\n" t
        ΔE += E
        @output "Model Space Energy           {:15.10f}\n" E + refwfn.molecule.Vnuc
        @output "Energy Change                {:15.10f}\n" ΔE

        if oldP == Set(P)
            cflag = true
            break
        end
        ite += 1
        if ite > 100
            break
        end
        oldP = Set(deepcopy(P))
        Lenny = length(P)
        @output "Coarse graining model space for next iteration\n"
        # Coarse grain
        Cperm = zeros(Int, length(P))
        sortperm!(Cperm, Pcoef, by=i->i^2)
        reverse!(Cperm)

        Pcoef = Pcoef[Cperm]
        P = P[Cperm]

        while true
            if sum(Pcoef[1:end-1].^2) >= 1-γ*σ
                pop!(Pcoef)
                pop!(P)
            else
                break
            end
        end
        @output "Final coarse grained model space size is {}\n" length(P)
        @output repeat("=",50)*"\n"
    end
    end #@elapsed

    @output repeat("=",50)*"\n"
    if cflag
        @output "🔥🔥🔥🔥🔥 ACI procedure has converged. 🔥🔥🔥🔥🔥\n"
    else
        @output "😲😲😲😲😲 ACI procedure has failed!!!! 😲😲😲😲😲\n"
    end
    @output "Computation finished in {:5.5} seconds.\n" ttotal
    @output "Model space size: {}\n" length(M)
    @output "E[ACI:{}]     = {:15.10f}\n" σ E + refwfn.molecule.Vnuc
    @output "E[ACI:{}]+PT2 = {:15.10f}\n" σ E + refwfn.molecule.Vnuc + ϵest
    @output repeat("=",51)*"\n\n"
    E = (E+refwfn.molecule.Vnuc)

    @output "\n • Most important determinants:\n\n"

    for i in 1:(min(20,length(P)))
        @output "{:15.5f}      {}\n" Pcoef[i]  detstring(P[i], frozen+active)
    end

    CASCI{T}(refwfn, E, P, Pcoef)
end

@fastmath @inbounds function get_fois(dets::Array{Determinant,1}, Nα::Int, Nβ::Int, act_range::UnitRange{Int64})::Array{Determinant,1}

    # Ns must be > 1
    αoccs = [zeros(Int,Nα) for i=1:Threads.nthreads()]
    βoccs = [zeros(Int,Nβ) for i=1:Threads.nthreads()]
    αunos = [zeros(Int,length(act_range)-Nα) for i=1:Threads.nthreads()]
    βunos = [zeros(Int,length(act_range)-Nβ) for i=1:Threads.nthreads()]

    lf_per_det = (length(αoccs[1])^2*length(αunos[1])^2 + length(αoccs[1])*length(αunos[1])
                       + length(βoccs[1])^2*length(βunos[1])^2 + length(βoccs[1])*length(βunos[1])
                       + length(αoccs[1])*length(αunos[1])*length(βoccs[1])*length(βunos[1]))
    lf_crit = Int(round(length(dets)*lf_per_det))
    fois = [Determinant(0,0) for i=1:lf_crit]
    @sync for _DI in eachindex(dets)
        Threads.@spawn begin
            d = dets[_DI]
            DI = (_DI-1)*lf_per_det + 1
            ct = 0
            id = Threads.threadid()
            αocc = αoccs[id]
            βocc = βoccs[id]
            αuno = αunos[id]
            βuno = βunos[id]
            αocc!(d, act_range, αocc)
            βocc!(d, act_range, βocc)
            αvir!(d, act_range, αuno)
            βvir!(d, act_range, βuno)
            # Get αα -> αα excitations
            for i in αocc
                for a in αuno
                    newα = (d.α ⊻ (1<<(i-1))) | (1<<(a-1))
                    _det = Determinant(newα, d.β)
                    fois[DI+ct] = _det
                    ct += 1
                    for j in αocc
                        if j ≥ i
                            break
                        end
                        for b in αuno
                            if b ≥ a
                                break
                            end
                            newestα = (newα ⊻ (1<<(j-1))) | (1<<(b-1))
                            _det = Determinant(newestα, d.β)
                            fois[DI+ct] = _det
                            ct += 1
                        end
                    end
                end
            end
            # Get ββ -> ββ excitations
            for i in βocc
                for a in βuno
                    newβ = (d.β ⊻ (1<<(i-1))) | (1<<(a-1))
                    _det = Determinant(d.α, newβ)
                    fois[DI+ct] = _det
                    ct += 1
                    for j in βocc
                        if j ≥ i
                            break
                        end
                        for b in βuno
                            if b ≥ a
                                break
                            end
                            newestβ = (newβ ⊻ (1<<(j-1))) | (1<<(b-1))
                            _det = Determinant(d.α, newestβ)
                            fois[DI+ct] = _det
                            ct += 1
                        end
                    end
                end
            end
            # Get αβ -> αβ excitations
            for i in αocc
                for a in αuno
                    for j in βocc
                        for b in βuno
                            newα = (d.α ⊻ (1<<(i-1))) | (1<<(a-1))
                            newβ = (d.β ⊻ (1<<(j-1))) | (1<<(b-1))
                            _det = Determinant(newα, newβ)
                            fois[DI+ct] = _det
                            ct += 1
                        end
                    end
                end
            end
        end #Threads.@spawn
    end
    fois = filter((x)->x != Determinant(0,0), fois)
    fois = Set(fois)
    setdiff!(fois, dets)
    fois = collect(fois)
    return fois end

#lexicographic premutations generation, By Donald Knuth
function lpermutations(a::BitArray)
  b=BitArray[]
  sort!(a)
  n=length(a)
  while(true)
    push!(b,copy(a))
    j=n-1
    while(a[j]>=a[j+1])
      j-=1
      j==0 && return(b)
    end
    l=n
    while(a[j]>=a[l])
      l-=1
    end
    tmp=a[l]
    a[l]=a[j]
    a[j]=tmp
    k=j+1
    l=n
    while(k<l)
      tmp=a[k]
      a[k]=a[l]
      a[l]=tmp
      k+=1
      l-=1
    end
  end
end

function ϵI(Fdets, P::Array{Determinant,1}, Pcoef::Array{Float64,1}, Ep::T, h::Array{T,2}, V::Array{T,4},act_elec,active) where T <: AbstractFloat
    Fe = zeros(length(Fdets))
    if false
        @output "Using Residue Array algorithm\n"
        N = count_ones(P[1].α) #number of electrons (assumes RHF)
        Ne = N*2
        Vints = zeros(length(Fdets))
        αind = zeros(Int64,N)
        βind = zeros(Int64,N)
        t = @elapsed for i in eachindex(Fdets)
            αindex!(Fdets[i], αind)
            βindex!(Fdets[i], βind)
            Ei = Hd0(αind, βind, h, V)
            Fe[i] = Ei - Ep #Δ
        end
        @output "time making Δ {}\n" t
        t = @elapsed begin
        mask = BitArray(undef,active)
        mask[1:2] .= 1
        masks_ss = lpermutations(mask)
        mask[2:2] .= 0
        masks_os = lpermutations(mask)
        for i in eachindex(masks_os)
            masks_os[i] = .~masks_os[i]
        end
        for i in eachindex(masks_ss)
            masks_ss[i] = .~masks_ss[i]
        end

        A = BitArray(undef,active)
        ae2 = Int(act_elec/2)
        A[1:ae2] .= 1
        oo = lpermutations(A)
        A[ae2-1:ae2] .= 0
        ss = lpermutations(A)
        aa = reshape(collect(Base.product(ss,oo)),(1,:))
        bb = reshape(collect(Base.product(oo,ss)),(1,:))
        A[ae2-1:ae2-1] .= 1
        os = (lpermutations(A))
        ab = reshape(collect(Base.product(os,os)),(1,:))
        ba = reshape(collect(Base.product(os,os)),(1,:))
        residues = hcat(aa,bb,ab,ba)
        residues = collect(Set(residues))
        end
        @output "time making perms {}\n" t

        # bit arrays for Fdet info
        DFa = BitArray(undef,active)
        DFb = BitArray(undef,active)
        DF1 = BitArray(undef,active)
        DF2 = BitArray(undef,active)

        ncore = Int((Ne - act_elec)/2)
        core = BitArray(undef,ncore)
        core[:] .= 1
        ct = 1
        E2 = Int(Ne/2 - 2)
        E1 = Int(Ne/2 - 1)

        t = @elapsed begin
        #resP = Dict(residues .=> [Set(Int64[]) for i in (residues)])
        resP = Array{Tuple}(undef,0)
        for i in eachindex(P)
            dF = P[i]
            DFa.chunks[1] = dF.α >> (ncore )
            DFb.chunks[1] = dF.β >> (ncore )
            for mask in masks_ss
                nua = DFa .& mask
                nub = DFb .& mask
                if (sum(nua) == E2 )#|| sum(nua) == E1)
                    push!(resP,((nua,DFb),i))
                    push!(resP,((DFb,nua),i))
                end
                if (sum(nub) == E2 )#|| sum(nub) == E1)
                    push!(resP,(((DFa),nub),i))
                    push!(resP,((nub,(DFa)),i))
                end
            end
            for mask1 in masks_os
                for mask2 in masks_os
                    nua = DFa .& mask1
                    nub = DFb .& mask2
                    if sum(nua) == E1 || sum(nub) == E1
                        push!(resP,((nua,nub),i))
                        push!(resP,((nub,nua),i))
                        push!(resP,((nua,nua),i))
                        push!(resP,((nub,nub),i))
                    end
                end
            end
        end
        end
        @output "time in P loop {}\n" t
        t = @elapsed begin
        resF = Array{Tuple}(undef,0)
        for i in eachindex(Fdets)
            dF = Fdets[i]
            DFa.chunks[1] = (dF.α >> (ncore))
            DFb.chunks[1] = (dF.β >> (ncore))
            for mask in masks_ss
                nua = DFa .& mask
                nub = DFb .& mask
                if (sum(nua) == E2 )#|| sum(nua) == E)
                    push!(resF,((nua,(DFb)),i))
                    push!(resF,(((DFb),nua),i))
                end
                if (sum(nub) == E2 )#|| sum(nub) == E2)
                    #push!(resP[(DFa,nub)],i)
                    push!(resF,(((DFa),nub),i))
                    push!(resF,((nub,(DFa)),i))
                end
            end
            for mask1 in masks_os
                for mask2 in masks_os
                    nua = DFa .& mask1
                    nub = DFb .& mask2
                    if sum(nua) == E1 || sum(nub) == E1
                        push!(resF,((nua,nub),i))
                        push!(resF,((nub,nua),i))
                        push!(resF,((nua,nua),i))
                        push!(resF,((nub,nub),i))
                    end
                end
            end
        end
        end
        @output "time in F loop {}\n" t
        t = @elapsed begin
        resP = collect(Set(resP))
        resF = collect(Set(resF))
        sort!(resP,by=first)
        sort!(resF,by=first)
        lastP,lastF = resP[1][1],resF[1][1]
        Fstart = 1
        Pstart = 1
        Fend = 1
        Fs = nothing
        detpairs = []
        end
        @output "setup time {}\n" t
        t = @elapsed begin
        pd = Dict()
        fd = Dict()
        for respair in resP
            try
                push!(pd[respair[1]],respair[2])
            catch KeyError
                pd[respair[1]] = []
                push!(pd[respair[1]],respair[2])
            end
        end
        for respair in resF
            try
                push!(fd[respair[1]],respair[2])
            catch KeyError
                fd[respair[1]] = []
                push!(fd[respair[1]],respair[2])
            end
        end

        #println(pd)
        #println(fd)

        for key in keys(pd)
            try
                push!(detpairs,collect(Base.product(fd[key],pd[key])))
            catch KeyError
                continue
            end
        end

        for idxF in eachindex(resF)
            if resF[idxF][1] != lastF
                Fs = last.(resF[Fstart:idxF-1])
                Fstart = idxF
                for idxP in Pstart:length(resP)
                    if resP[idxP][1] != lastP
                        Ps = last.(resP[Pstart:idxP-1])
                        Pstart = idxP
                        prod = collect(Base.product(Fs,Ps))
                        if length(prod) != 0
                            push!(detpairs,collect(Base.product(Fs,Ps)))
                        end
                        lastP = resP[idxP][1]
                        Ps = nothing
                        break
                    end
                end
                lastF = resF[idxF][1]
            end
            Fs = nothing
        end
        end
        @output "time spent pairing {}\n" t

        t = @elapsed begin
        detpairs = [reshape(dp,1,:) for dp in detpairs]
        detpairs = hcat(detpairs...)
        detpairs = collect(Set(detpairs))
        println(length(detpairs))
        #perm = sortperm(detpairs,by=first)
        #detpairs = detpairs[perm]
        sort!(detpairs,by=first)
        end
        @output "cleanup time {}\n" t
        ct = 0
        t = @elapsed for detpair in detpairs
            i = detpair[1]
            j = detpair[2]
            D1 = Fdets[i]
            αindex!(D1, αind)
            βindex!(D1, βind)
            D2 = P[j]
            αexc = αexcitation_level(D1,D2)
            βexc = βexcitation_level(D1,D2)
            el = αexc + βexc
            if el == 2
                Vints[i] += Pcoef[j]*Hd2(D1, D2, V, αexc)
            elseif el == 1
                Vints[i] += Pcoef[j]*Hd1(αind, βind, D1, D2, h, V, αexc)
            else
                ct += 1
            end
        end
        @output "time in P space {}\n" t
        @output "There were {} false positives\n" ct
        for i in eachindex(Fe)
            Fe[i] = Fe[i]/2 - √((Fe[i]^2)/4 + Vints[i]^2)
        end
    else
        @output "Using double loop algorithm\n"
        N = sum(αlist(P[1]))
        αinds = [Array{Int64,1}(undef,N) for i=1:Threads.nthreads()]
        βinds = [Array{Int64,1}(undef,N) for i=1:Threads.nthreads()]
        @sync for i in eachindex(Fdets)
            begin
            D1 = Fdets[i]
            id = Threads.threadid()
            αind = αinds[id]
            βind = βinds[id]
            αindex!(D1, αind)
            βindex!(D1, βind)
            Ei = Hd0(αind, βind, h, V)
            Δ = Ei - Ep
            Vint = 0.0
            for j in eachindex(P)
                D2 = P[j]
                αexc = αexcitation_level(D1,D2)
                βexc = βexcitation_level(D1,D2)
                el = αexc + βexc
                if el > 2
                    continue
                elseif el == 2
                    Vint += Pcoef[j]*Hd2(D1, D2, V, αexc)
                elseif el == 1
                    Vint += Pcoef[j]*Hd1(αind, βind, D1, D2, h, V, αexc)
                end
            end

            @fastmath Fe[i] = Δ/2 - √((Δ^2)/4 + Vint^2)
            end #Threads.@spawn
        end
    end
    return Fe
end

function update_model_space(M::Array{Determinant,1}, h::Array{T,2}, V::Array{T,4}) where T <: AbstractFloat

    M = complete_set(M)
    H = get_sparse_hamiltonian_matrix(M, h, V, Fermi.CurrentOptions["cas_cutoff"])

    @output "Diagonalizing Hamiltonian...\n"
    decomp, history = partialschur(H, nev=1, tol=10^-12, which=LM())
    λ, ϕ = partialeigen(decomp)
    #λ,ϕ = eigen(Array(H))

    return λ[1], ϕ[:,1], deepcopy(M)
end

function incomp(dets::Array{Determinant,1})

    ming = 0
    for d in dets

        _det = Determinant(d.β, d.α)
        if !(_det in dets)
            ming += 1
        end
    end

    return ming
end

function complete_set(dets::Array{Determinant,1})

    newdets = [Determinant[] for i = 1:Threads.nthreads()]
    @Threads.threads for d in dets

        asym = d.α ⊻ d.β
        if asym == 0
            continue
        end

        sym = d.α & d.β

        n = count_ones(asym)
        e = Int(n/2)
        idx = Int[]

        str = vcat(repeat([1],e), repeat([0],e))
        perms = multiset_permutations(str, n)

        i = 1
        while i ≤ asym
            if 1<<(i-1) & asym ≠ 0
                push!(idx, i)
            end
            i += 1
        end

        for p in perms
            newα = sym
            newβ = sym
            for (x,i) in zip(p,idx)
                if x == 1
                    newα = newα | (1<<(i-1))
                elseif x == 0
                    newβ = newβ | (1<<(i-1))
                end
            end
            push!(newdets[Threads.threadid()], Determinant(newα, newβ))
        end
    end
    newdets = vcat(newdets...)
    return unique(vcat(dets,newdets))
end

function woz(A::Vector{T}, sentinel::T) where T
         p = sortperm(A)
         q = A[p]
         res = Vector{Vector{T}}()
         grp = Vector{T}()
         first = true
         last = sentinel
         for (i,v) in enumerate(q)
           if !first && last != v
             push!(res, grp)
             grp = [p[i]]
           else
             push!(grp, p[i])
           end
           last = v
           first = false
         end
         push!(res, grp)
         res
       end
groupby(f, list::Array) = begin
  foldl(list; init = Dict()) do dict, v
    push!(get!(dict, f(v), []), v)
    dict
  end
end


text:
using RICoverAlg
using Test

@testset "Testing Preprocessing" begin include("preprocessinput_test.jl") end
@testset "Testing algorithm" begin include("cover_algorithm_test.jl") end
@testset "Testing whole package" begin include("RICoverAlg_test.jl") end
# @testset "Is the package still fast" begin include("benchmarking.jl") end


text:
module RATTLE

export RattleIntegrator
include("./types.jl")
include("./main.jl")

end


text:
# ------------------------------------
# time integration with explicit methods
# ------------------------------------
function time_integration_explicit(dt, Qcon_hat, RHS, cellxmax, cellymax, cellzmax, nval, icell)
    Threads.@threads for l in 1:nval
        for k in 1+icell:cellzmax-icell
            for j in 1+icell:cellymax-icell
                for i in 1+icell:cellxmax-icell
                    Qcon_hat[i,j,k,l] = Qcon_hat[i,j,k,l] + dt*RHS[i,j,k,l]
                end
            end
        end
    end
    return Qcon_hat
end

# ------------------------------------
# Calculation of the advection Jacobian matrix by one-wave approximation for lusgs
# ------------------------------------
function one_wave(A_adv_hat_m, A_adv_hat_p, B_adv_hat_m, B_adv_hat_p, C_adv_hat_m, C_adv_hat_p, A_beta_shig, B_beta_shig, C_beta_shig, I,
                    Qbase, Qcon, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, specific_heat_ratio, volume, nval)
    beta = 1.1
    kx_av = 0.0
    ky_av = 0.0
    kz_av = 0.0

    for k in 2:cellzmax
        for j in 2:cellymax
            for i in 2:cellxmax
                for nn in 1:2   #A, B, C
                    jacob_temp = zeros(nval, nval)
                    if nn == 1
                        kx_av = 0.5*(vecAx[i,j,k,1]+vecAx[i+1,j,k,1]) / volume[i,j,k]
                        ky_av = 0.5*(vecAx[i,j,k,2]+vecAx[i+1,j,k,2]) / volume[i,j,k]
                        kz_av = 0.5*(vecAx[i,j,k,3]+vecAx[i+1,j,k,3]) / volume[i,j,k]
                    elseif nn == 2
                        kx_av = 0.5*(vecAy[i,j,k,1]+vecAy[i,j+1,k,1]) / volume[i,j,k]
                        ky_av = 0.5*(vecAy[i,j,k,2]+vecAy[i,j+1,k,2]) / volume[i,j,k]
                        kz_av = 0.5*(vecAy[i,j,k,3]+vecAy[i,j+1,k,3]) / volume[i,j,k]
                    elseif nn == 3
                        kx_av = 0.5*(vecAz[i,j,k,1]+vecAz[i,j,k+1,1]) / volume[i,j,k]
                        ky_av = 0.5*(vecAz[i,j,k,2]+vecAz[i,j,k+1,2]) / volume[i,j,k]
                        kz_av = 0.5*(vecAz[i,j,k,3]+vecAz[i,j,k+1,3]) / volume[i,j,k]
                    end

                    rho = Qbase[i,j,k,1]
                    u = Qbase[i,j,k,2]
                    v = Qbase[i,j,k,3]
                    w = Qbase[i,j,k,4]
                    e = Qcon[i,j,k,5]
                    p = Qbase[i,j,k,5]

                    g = specific_heat_ratio

                    Z = kx_av*u + ky_av*v + kz_av*w
                    q2   = 0.5*(u^2 + v^2 + w^2)
                    b1c2 = 0.5*q2*(g-1)
                    gebyrho = g*e/rho

                    # advection Jacobian matrix in the general coordinate system
                    jacob_temp[1,1] = 0.0
                    jacob_temp[1,2] = kx_av
                    jacob_temp[1,3] = ky_av
                    jacob_temp[1,4] = kz_av
                    jacob_temp[1,5] = 0.0

                    jacob_temp[2,1] = -u*Z + kx_av*b1c2
                    jacob_temp[2,2] = Z - (g-2)*kx_av*u
                    jacob_temp[2,3] = ky_av*u - kx_av*(g-1)*v
                    jacob_temp[2,4] = kz_av*u - kx_av*(g-1)*w
                    jacob_temp[2,5] = (g-1)*kx_av

                    jacob_temp[3,1] = -v*Z + ky_av*b1c2
                    jacob_temp[3,2] = kx_av*v - ky_av*(g-1)*u
                    jacob_temp[3,3] = Z - ky_av*(g-2)*v
                    jacob_temp[3,4] = kz_av*v - ky_av*(g-1)*w
                    jacob_temp[3,5] = (g-1)*ky_av

                    jacob_temp[4,1] = -w*Z + kz_av*b1c2
                    jacob_temp[4,2] = kx_av*w - kz_av*(g-1)*u
                    jacob_temp[4,3] = ky_av*w - kz_av*(g-1)*v
                    jacob_temp[4,4] = Z - kz_av*(g-2)*w
                    jacob_temp[4,5] = (g-1)*kz_av

                    jacob_temp[5,1] = Z*(-gebyrho + 2*b1c2)
                    jacob_temp[5,2] = kx_av*(gebyrho-b1c2) - (g-1)*u*Z
                    jacob_temp[5,3] = ky_av*(gebyrho-b1c2) - (g-1)*v*Z
                    jacob_temp[5,4] = kz_av*(gebyrho-b1c2) - (g-1)*w*Z
                    jacob_temp[5,5] = g*Z

                    c = (g*rho/p)^0.5
                    shigma = abs(Z) + c*(kx_av^2 + ky_av^2 + kz_av^2)^0.5

                    for l in 1:nval
                        I[l,l] = beta * shigma
                    end

                    if nn == 1
                        for m in 1:nval
                            for l in 1:nval
                                A_adv_hat_p[i,j,k,l,m] = 0.5*(jacob_temp[l,m] + I[l,m])
                                A_adv_hat_m[i,j,k,l,m] = 0.5*(jacob_temp[l,m] - I[l,m])
                            end
                        end
                        A_beta_shig[i,j,k] = beta * shigma
                    elseif nn == 2
                        for m in 1:nval
                            for l in 1:nval
                                B_adv_hat_p[i,j,k,l,m] = 0.5*(jacob_temp[l,m] + I[l,m])
                                B_adv_hat_m[i,j,k,l,m] = 0.5*(jacob_temp[l,m] - I[l,m])
                            end
                        end
                        B_beta_shig[i,j,k] = beta * shigma
                    elseif nn == 3
                        for m in 1:nval
                            for l in 1:nval
                                C_adv_hat_p[i,j,k,l,m] = 0.5*(jacob_temp[l,m] + I[l,m])
                                C_adv_hat_m[i,j,k,l,m] = 0.5*(jacob_temp[l,m] - I[l,m])
                            end
                        end
                        C_beta_shig[i,j,k] = beta * shigma
                    end

                end
            end
        end
    end

    # reset
    for l in 1:nval
        I[l,l] = 1.0
    end

    return A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, C_adv_hat_p, C_adv_hat_m, A_beta_shig, B_beta_shig, C_beta_shig
end

# ------------------------------------
# Calculation of the viscosity Jacobian matrix by center differential for lusgs
# ------------------------------------
function central_diff_jacobian(jalphaP, jbetaP, jgammaP, Qbase, Qcon, cellxmax, cellymax, cellzmax, mu, lambda,
                            vecAx, vecAy, vecAz, specific_heat_ratio, volume, nval)

    for k in 1:cellzmax
        for j in 1:cellymax
            for i in 1:cellxmax
                rho = Qbase[i,j,k,1]

                # alpha
                xi_x_av = 0.5*(vecAx[i,j,k,1] + vecAx[i+1,j,k,1]) / volume[i,j,k]
                xi_y_av = 0.5*(vecAx[i,j,k,2] + vecAx[i+1,j,k,2]) / volume[i,j,k]
                xi_z_av = 0.5*(vecAx[i,j,k,3] + vecAx[i+1,j,k,3]) / volume[i,j,k]
                alpha = (xi_x_av^2 + xi_y_av^2 + xi_z_av^2)^0.5
                jalphaP[i,j,k] = alpha / volume[i,j,k] * (2*mu[i,j,k]/rho)

                # beta
                eta_x_av = 0.5*(vecAy[i,j,k,1] + vecAy[i,j+1,k,1]) / volume[i,j,k]
                eta_y_av = 0.5*(vecAy[i,j,k,2] + vecAy[i,j+1,k,2]) / volume[i,j,k]
                eta_z_av = 0.5*(vecAy[i,j,k,3] + vecAy[i,j+1,k,3]) / volume[i,j,k]
                beta = (eta_x_av^2 + eta_y_av^2 + eta_z_av^2)^0.5
                jbetaP[i,j,k] = beta / volume[i,j,k] * (2*mu[i,j,k]/rho)

                # gamma
                zeta_x_av = 0.5*(vecAz[i,j,k,1] + vecAz[i,j,k+1,1]) / volume[i,j,k]
                zeta_y_av = 0.5*(vecAz[i,j,k,2] + vecAz[i,j,k+1,2]) / volume[i,j,k]
                zeta_z_av = 0.5*(vecAz[i,j,k,3] + vecAz[i,j,k+1,3]) / volume[i,j,k]
                gamma = (zeta_x_av^2 + zeta_y_av^2 + zeta_z_av^2)^0.5
                jgammaP[i,j,k] = gamma / volume[i,j,k] * (2*mu[i,j,k]/rho)
            end
        end
    end
    return jalphaP, jbetaP, jgammaP
end

# ------------------------------------
# lusgs
# ------------------------------------
function lusgs(D, Lx, Ly, Lz, Ux, Uy, Uz, LdQ, UdQ, RHS_temp, I, dt, dtau, Qcon_hat, Qconn_hat, delta_Q,
                A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, C_adv_hat_p, C_adv_hat_m, A_beta_shig, B_beta_shig, C_beta_shig,
                jalphaP, jbetaP, jgammaP, RHS, cellxmax, cellymax, cellzmax, volume, nval, icell)

    # calculate L and U
    for m in 1:nval
        for l in 1:nval
            for k in 1:cellzmax
                for j in 1:cellymax
                    for i in 1:cellxmax
                        Lx[i,j,k,l,m] = dt*(A_adv_hat_p[i,j,k,l,m] + jalphaP[i,j,k]*I[l,m])
                        Ly[i,j,k,l,m] = dt*(B_adv_hat_p[i,j,k,l,m] + jbetaP[i,j,k]*I[l,m])
                        Lz[i,j,k,l,m] = dt*(C_adv_hat_p[i,j,k,l,m] + jgammaP[i,j,k]*I[l,m])
                        Ux[i,j,k,l,m] = dt*(A_adv_hat_m[i,j,k,l,m] - jalphaP[i,j,k]*I[l,m])
                        Uy[i,j,k,l,m] = dt*(B_adv_hat_m[i,j,k,l,m] - jbetaP[i,j,k]*I[l,m])
                        Uz[i,j,k,l,m] = dt*(C_adv_hat_m[i,j,k,l,m] - jgammaP[i,j,k]*I[l,m])
                    end
                end
            end
        end
    end

    for k in 1+icell:cellzmax-icell
        for j in 1+icell:cellymax-icell
            for i in 1+icell:cellxmax-icell
                for l in 1:nval
                    for m in 1:nval
                        LdQ[i,j,k,l] += Lx[i-1,j,k,l,m]*delta_Q[i-1,j,k,m]
                                      + Ly[i,j-1,k,l,m]*delta_Q[i,j-1,k,m]
                                      + Lz[i,j,k-1,l,m]*delta_Q[i,j,k-1,m]
                        UdQ[i,j,k,l] += Ux[i+1,j,k,l,m]*delta_Q[i+1,j,k,m]
                                      + Uy[i,j+1,k,l,m]*delta_Q[i,j+1,k,m]
                                      + Uz[i,j,k+1,l,m]*delta_Q[i,j,k+1,m]
                    end
                end
            end
        end
    end

    # diagonal
    for k in 1+icell:cellzmax-icell
        for j in 1+icell:cellymax-icell
            for i in 1+icell:cellxmax-icell
                D[i,j,k] = dt/dtau[i,j,k] + 1.0 + dt*(A_beta_shig[i,j,k] + 2*jalphaP[i,j,k] +
                                                      B_beta_shig[i,j,k] + 2*jbetaP[i,j,k] +
                                                      C_beta_shig[i,j,k] + 2*jgammaP[i,j,k])
            end
        end
    end

    # RHS
    for l in 1:nval
        for k in 1+icell:cellzmax-icell
            for j in 1+icell:cellymax-icell
                for i in 1+icell:cellxmax-icell
                    RHS_temp[i,j,k,l] = - (Qcon_hat[i,j,k,l]-Qconn_hat[i,j,k,l])*1.0 + dt*RHS[i,j,k,l]
                end
            end
        end
    end

    # lower sweep
    for l in 1:nval
        for k in 1+icell:cellzmax-icell
            for j in 1+icell:cellymax-icell
                for i in 1+icell:cellxmax-icell
                    delta_Q[i,j,k,l] = D[i,j,k]^(-1) * (RHS_temp[i,j,k,l]+LdQ[i,j,k,l])
                end
            end
        end
    end

    # upepr sweep
    for l in 1:nval
        for k in 1+icell:cellzmax-icell
            for j in 1+icell:cellymax-icell
                for i in 1+icell:cellxmax-icell
                    delta_Q[i,j,k,l] = delta_Q[i,j,k,l] - D[i,j,k]^(-1) * UdQ[i,j,k,l]
                end
            end
        end
    end

    # reset
    for l in 1:nval
        for k in 1+icell:cellzmax-icell
            for j in 1+icell:cellymax-icell
                for i in 1+icell:cellxmax-icell
                    LdQ[i,j,k,l] = 0.0
                    UdQ[i,j,k,l] = 0.0
                end
            end
        end
    end

    return delta_Q
end


text:
#-------------------------------------------------------------------------------
#transformed coordinates are handy

#upward calculations
P2ω(P) = -sqrt(P)
ω2P(ω) = ω*ω
dωfac(P) = 2*sqrt(P)
P2ω(Pₛ, Pₜ) = P2ω(Pₛ), P2ω(Pₜ)

#downward calculations
P2ι(P) = sqrt(P)
ι2P(ι) = ι*ι
dιfac(P) = 2*sqrt(P)
P2ι(Pₜ, Pₛ) = P2ι(Pₜ), P2ι(Pₛ)

#-------------------------------------------------------------------------------

export pressuregrid
function pressuregrid(Pₜ, Pₛ, n)
    @assert Pₛ > Pₜ
    @assert n >= 3
    exp.(chebygrid(log(Pₜ), log(Pₛ), n))
end

export trapz
function trapz(x::AbstractVector, y::AbstractVector)
    @assert length(x) == length(y) "vectors must be equal length"
    s = zero(eltype(y))
    @inbounds for i ∈ 1:length(x) - 1
        s += (x[i+1] - x[i])*(y[i] + y[i+1])/2
    end
    return s
end

export meshgrid
function meshgrid(x::AbstractVector, y::AbstractVector)
    X = x' .* ones(length(y))
    Y = y .* ones(length(x))'
    return X, Y
end

export logrange
function logrange(a, b, N::Int=101, γ::Real=1)
    ((10 .^ LinRange(0, γ, N)) .- 1)*(b - a)/(10^γ - 1) .+ a
end

#evaluates the 0th, 1st, and 2nd derivative of a parabola through three points, anywhere
function quaddiff(x, x₁, x₂, x₃, y₁, y₂, y₃)
    z₁ = x₁^2 - x₂^2
    z₂ = x₂^2 - x₃^2
    w = z₁/z₂
    b = (y₁ - y₂ - (y₂ - y₃)*w)/(x₁ - x₂ - (x₂ - x₃)*w)
    a = (y₂ - y₃ - b*(x₂ - x₃))/z₂
    c = y₁ - (a*x₁^2 + b*x₁)
    f = a*x^2 + b*x + c
    f′ = 2*a*x + b
    f′′ = 2*a
    return f, f′, f′′
end

function quaddiff(x, xₚ, yₚ)
    @assert length(xₚ) == length(yₚ) == 3
    @inbounds quaddiff(x, xₚ[1], xₚ[2], xₚ[3], yₚ[1], yₚ[2], yₚ[3])[2]
end

export deriv
function deriv(x, y, T)
    d = similar(x)
    d[1] = (y[2] - y[1])/(x[2] - x[1])
    for i in 2:length(d)-1
        #d₁ = (y[i] - y[i-1])/(x[i] - x[i-1])
        #d₂ = (y[i+1] - y[i])/(x[i+1] - x[i])
        #if i == 25
        #    d[i] = (d₁ + d₂)/2
        #elseif i > 25
        #    d[i] = d₁
        #else
        #    d[i] = d₂
        #end
        #if d₁*d₂ > 0
        #    d[i] = (d₁ + d₂)/2
        #else
        #    if T[i-1] < T[i+1]
        #        d[i] = d₁
        #    else
        #        d[i] = d₂
        #    end
        #end
        d[i] = quaddiff(x[i], view(x,i-1:i+1), view(y,i-1:i+1))
    end
    d[end] = (y[end] - y[end-1])/(x[end] - x[end-1])
    return d
end

#------------------------------------------------------------------------------
# a couple of root finding methods

terminate(a, b, tol)::Bool = abs(a - b) < (tol + tol*abs(b)) ? true : false

function terminate(x₁, x₂, y₁, y₂, tol)::Bool
    terminate(x₁, x₂, tol) && terminate(y₁, y₂, tol) && return true
    return false
end

export regulafalsi
function regulafalsi(F, x₁, x₂, p=nothing; tol=1e-6)
    @assert x₁ != x₂ "starting points must not be identical"
    y₁ = F(x₁, p)
    y₁ == 0 && return x₁
    y₂ = F(x₂, p)
    y₂ == 0 && return x₂
    @assert sign(y₁) != sign(y₂) "regula falsi non-bracketing"
    yₘ = floatmax(y₂)
    yₚ = zero(y₂)
    n = 0
    while !terminate(x₁, x₂, yₚ, yₘ, tol) || (n < 2)
        #store previous evaluation
        yₚ = yₘ
        #approximate zero
        xₘ = x₁ - y₁*(x₂ - x₁)/(y₂ - y₁)
        yₘ = F(xₘ, p)
        #reduce the interval
        if y₁*yₘ > 0
            x₁ = xₘ
        else
            x₂ = xₘ
        end
        #count
        n += 1
    end
    return (x₁ + x₂)/2
end

export secant
function secant(F, x₁, x₂, p=nothing; tol=1e-6)
    @assert x₁ != x₂ "starting points must not be identical"
    y₁ = F(x₁, p)
    y₁ == 0 && return x₁
    y₂ = F(x₂, p)
    y₂ == 0 && return x₂
    y₃ = floatmax(y₂)
    x₃ = zero(x₂)
    n = 0
    while !terminate(x₁, x₂, y₁, y₂, tol) || (n < 2)
        #approximate zero
        x₃ = x₁ - y₁*(x₂ - x₁)/(y₂ - y₁)
        y₃ = F(x₃, p)
        #swap values
        x₁, x₂ = x₂, x₃
        y₁, y₂ = y₂, y₃
        #count
        n += 1
    end
    return x₃
end