text:
"""
    Ormsby(; <keyword arguments>)

Create a Ormsby wavelet sampled every dt seconds with corner frequencies
defined by the vector f = [f1, f2, f3, f4]. The final wavelet is multiplied by
a Hamming window.

# Arguments

**Keyword arguments**

* `dt::Real=0.002`: sampling interval in secs.
* `f::Vector{Real}=[2.0, 10.0, 40.0, 60.0]`: corner frequencies in Hz.

      ^
    1 |     ***************
      |    *               *
      |   *                 *
      |  *                   *
      | *                     *
      -----------------------------> f
        f1  f2           f3  f4

# Example
```julia
julia> w = Ormsby(); plot(w);
```
"""

function Ormsby{T<:Real}(; dt::Real=0.002, f::Vector{T}=[2.0, 10.0, 40.0, 60.0])

    f1 = f[1]
    f2 = f[2]
    f3 = f[3]
    f4 = f[4]

    fc = (f2+f3)/2.0
    nw = 2.2/(fc*dt)
    nc = floor(Int, nw/2)
    t = dt*collect(-nc:1:nc)
    nw = 2*nc + 1
    a4 = (pi*f4)^2/(pi*(f4-f3))
    a3 = (pi*f3)^2/(pi*(f4-f3))
    a2 = (pi*f2)^2/(pi*(f2-f1))
    a1 = (pi*f1)^2/(pi*(f2-f1))

    u = a4*(sinc(f4*t)).^2 - a3*(sinc(f3*t)).^2
    v = a2*(sinc(f2*t)).^2 - a1*(sinc(f1*t)).^2

    w = u - v
    w = w.*Hamming(nw)/maximum(w)

end


text:
@assert ends_inside(23, -10, day17_test)
@assert ends_inside(25, -9, day17_test)
@assert ends_inside(27, -5, day17_test)
@assert ends_inside(29, -6, day17_test)
@assert ends_inside(22, -6, day17_test)
@assert ends_inside(21, -7, day17_test)
@assert ends_inside(9, 0, day17_test)
@assert ends_inside(27, -7, day17_test)
@assert ends_inside(24, -5, day17_test)
@assert ends_inside(25, -7, day17_test)
@assert ends_inside(26, -6, day17_test)
@assert ends_inside(25, -5, day17_test)
@assert ends_inside(6, 8, day17_test)
@assert ends_inside(11, -2, day17_test)
@assert ends_inside(20, -5, day17_test)
@assert ends_inside(29, -10, day17_test)
@assert ends_inside(6, 3, day17_test)
@assert ends_inside(28, -7, day17_test)
@assert ends_inside(8, 0, day17_test)
@assert ends_inside(30, -6, day17_test)
@assert ends_inside(29, -8, day17_test)
@assert ends_inside(20, -10, day17_test)
@assert ends_inside(6, 7, day17_test)
@assert ends_inside(6, 4, day17_test)
@assert ends_inside(6, 1, day17_test)
@assert ends_inside(14, -4, day17_test)
@assert ends_inside(21, -6, day17_test)
@assert ends_inside(26, -10, day17_test)
@assert ends_inside(7, -1, day17_test)
@assert ends_inside(7, 7, day17_test)
@assert ends_inside(8, -1, day17_test)
@assert ends_inside(21, -9, day17_test)
@assert ends_inside(6, 2, day17_test)
@assert ends_inside(20, -7, day17_test)
@assert ends_inside(30, -10, day17_test)
@assert ends_inside(14, -3, day17_test)
@assert ends_inside(20, -8, day17_test)
@assert ends_inside(13, -2, day17_test)
@assert ends_inside(7, 3, day17_test)
@assert ends_inside(28, -8, day17_test)
@assert ends_inside(29, -9, day17_test)
@assert ends_inside(15, -3, day17_test)
@assert ends_inside(22, -5, day17_test)
@assert ends_inside(26, -8, day17_test)
@assert ends_inside(25, -8, day17_test)
@assert ends_inside(25, -6, day17_test)
@assert ends_inside(15, -4, day17_test)
@assert ends_inside(9, -2, day17_test)
@assert ends_inside(15, -2, day17_test)
@assert ends_inside(12, -2, day17_test)
@assert ends_inside(28, -9, day17_test)
@assert ends_inside(12, -3, day17_test)
@assert ends_inside(24, -6, day17_test)
@assert ends_inside(23, -7, day17_test)
@assert ends_inside(25, -10, day17_test)
@assert ends_inside(7, 8, day17_test)
@assert ends_inside(11, -3, day17_test)
@assert ends_inside(26, -7, day17_test)
@assert ends_inside(7, 1, day17_test)
@assert ends_inside(23, -9, day17_test)
@assert ends_inside(6, 0, day17_test)
@assert ends_inside(22, -10, day17_test)
@assert ends_inside(27, -6, day17_test)
@assert ends_inside(8, 1, day17_test)
@assert ends_inside(22, -8, day17_test)
@assert ends_inside(13, -4, day17_test)
@assert ends_inside(7, 6, day17_test)
@assert ends_inside(28, -6, day17_test)
@assert ends_inside(11, -4, day17_test)
@assert ends_inside(12, -4, day17_test)
@assert ends_inside(26, -9, day17_test)
@assert ends_inside(7, 4, day17_test)
@assert ends_inside(24, -10, day17_test)
@assert ends_inside(23, -8, day17_test)
@assert ends_inside(30, -8, day17_test)
@assert ends_inside(7, 0, day17_test)
@assert ends_inside(9, -1, day17_test)
@assert ends_inside(10, -1, day17_test)
@assert ends_inside(26, -5, day17_test)
@assert ends_inside(22, -9, day17_test)
@assert ends_inside(6, 5, day17_test)
@assert ends_inside(7, 5, day17_test)
@assert ends_inside(23, -6, day17_test)
@assert ends_inside(28, -10, day17_test)
@assert ends_inside(10, -2, day17_test)
@assert ends_inside(11, -1, day17_test)
@assert ends_inside(20, -9, day17_test)
@assert ends_inside(14, -2, day17_test)
@assert ends_inside(29, -7, day17_test)
@assert ends_inside(13, -3, day17_test)
@assert ends_inside(23, -5, day17_test)
@assert ends_inside(24, -8, day17_test)
@assert ends_inside(27, -9, day17_test)
@assert ends_inside(30, -7, day17_test)
@assert ends_inside(28, -5, day17_test)
@assert ends_inside(21, -10, day17_test)
@assert ends_inside(7, 9, day17_test)
@assert ends_inside(6, 6, day17_test)
@assert ends_inside(21, -5, day17_test)
@assert ends_inside(27, -10, day17_test)
@assert ends_inside(7, 2, day17_test)
@assert ends_inside(30, -9, day17_test)
@assert ends_inside(21, -8, day17_test)
@assert ends_inside(22, -7, day17_test)
@assert ends_inside(24, -9, day17_test)
@assert ends_inside(20, -6, day17_test)
@assert ends_inside(6, 9, day17_test)
@assert ends_inside(29, -5, day17_test)
@assert ends_inside(8, -2, day17_test)
@assert ends_inside(27, -8, day17_test)
@assert ends_inside(30, -5, day17_test)
@assert ends_inside(24, -7, day17_test)

ext: jl


text:
using Test

@time @testset "Assemblers" begin include("test_assemblers.jl") end
@time @testset "Reference shapes" begin include("test_refshapes.jl") end
@time @testset "Finite elements" begin include("test_felements.jl") end
@time @testset "Fields" begin include("test_fefields.jl") end
@time @testset "Spaces" begin include("test_fespaces.jl") end
@time @testset "Quadrature iterators" begin include("test_qpiterators.jl") end
@time @testset "FE iterators" begin include("test_feiterators.jl") end
@time @testset "Heat" begin include("test_heat.jl") end
@time @testset "Elasticity" begin include("test_elasticity.jl") end
@time @testset "Stokes" begin include("test_stokes.jl") end


text:
module Hertz

using JuLIP
using JuLIP.Potentials: pairs
import JuLIP: energy, forces, cutoff
using NeighbourLists
using StaticArrays

const SVecG{T} = SVector{3, T}
const SMatG{T} = SMatrix{3, 3, T, 9}

export HHZPotential

@pot struct HP <: AbstractCalculator #HP means Hooke/Hertzian Potential
   V::PairPotential
   rcut::Float64
   rtol::Float64
end

# mutable struct mynlist
#     i::Vector{Int}
#     j::Vector{Int}
#     r::Vector{Float64}
#     R::Vector{JVecF}
#     s::Vector{Float64}
# end

cutoff(V::HP) = V.rcut
"""
Contact Model
k is the energy scale,
here, s = r - r0, r0 is the sum of radii of two grains
α = 2  => harmonic spring model
α = 2.5 => Hertzian contact model
"""
HHZ(k, α) = @analytic s -> k * (-s)^α # harmonic or Hertzian potential
@inline evaluate_d_HHZ(p::Shift{-1}, s) = r < p.rcut ? (@D p.V(s)) : 0.0
@inline grad_HHZ(p::PairPotential, s::T, r::T, R::JVec{T}) where {T <: Real} = (evaluate_d_HHZ(p, s) / r) * R

@inline function sum_energy_HHZ(V::PairPotential, s::Vector{T}) where {T <: Real}
    E = 0.0
    @simd for n = 1:length(s)
        @inbounds E += V(s[n])
    end
    return E
end

function HHZPotential(k, α, r0::Float64; rcutfact = 2.5)
    rcut = rcutfact * r0 # rcut is rcutfact*2*R_{max}, used when searching possible pairs
    rtol = 0.5 * rcut    # rtol is the neighbourlist recomputing tolerance
    V = HHZ(k, α) * HS(0.0)
    return HP(V, rcut, rtol)
end

function energy(V::HP, at::AbstractAtoms)
    mynlist = update_neighbourlist!(V, at)
    return 0.5*sum_energy_HHZ(V.V, get_data(at, :nlist_s))
end
# compute the maximum displacement of all particles
@inline max_dr2(X::Vector{JVecF}, X_Old::Vector{JVecF}) =
                    maximum( [ sum( (X[i] - X_Old[i]) .* (X[i] - X_Old[i]) ) for i = 1:length(X) ] )
@inline s_r_r0(i::Vector{Int}, j::Vector{Int}, Rad::Vector{Float64}, r::Vector{Float64}) =
                    [ r[k] - Rad[i[k]] - Rad[j[k]] for k = 1:length(i) ]

function update_neighbourlist!(V::HP, at::AbstractAtoms)
    if has_data(at, :nlist_X)
        X_Old = get_data(at, :nlist_X)::Vector{JVecF}

        X = at.X::Vector{JVecF}
        max_d = max_dr2(X, X_Old)
        if max_d < V.rtol*V.rtol
            nlist_i = get_data(at, :nlist_i)::Vector{Int}
            nlist_j = get_data(at, :nlist_j)::Vector{Int}
            update_distance!(nlist_i, nlist_j, at)
            return nothing
        end
    end
    nlist = neighbourlist(at, V.rcut)
    Rad = get_data(at, :Rad)
    nlist_s = s_r_r0(nlist.i, nlist.j, Rad, nlist.r)
    set_data!(at, :nlist_X, copy(at.X))
    set_data!(at, :nlist_i, nlist.i)
    set_data!(at, :nlist_j, nlist.j)
    set_data!(at, :nlist_r, nlist.r)
    set_data!(at, :nlist_R, nlist.R)
    set_data!(at, :nlist_s, nlist_s)
#    return mynlist(nlist.i, nlist.j, nlist.r, nlist.R, nlist_s)
end

@inline wrap_to_unit(inv_bins::SMatG{T}, x::SVecG{T}) where {T <: Real} = mod.(inv_bins * x, 1)
@inline adjust_dx(bins::SMatG{T}, dx::SVecG{T}) where {T <: Real} = bins * (dx - (dx .> 0.5))
# @inline x2dx(xi::SVecG{T}, xj::SVecG{T}, inv_bins::SMatG{T}, bins::SMatG{T}) where {T <: Real} =
#                                     adjust_dx(bins, wrap_to_unit())
"""
Note: The following function is correct only when periodic boundary condition
"""
function update_distance!(nlist_i::Vector{TI}, nlist_j::Vector{TI}, at::AbstractAtoms) where TI <: Integer
    nn = length(nlist_i)
    nlist_r = Vector{Float64}()
    nlist_R = Vector{JVecF}()
    nlist_s = Vector{Float64}()
    X = at.X
    Rad = get_data(at, :Rad)
    bins = at.cell'
    inv_bins = inv(bins)
    for k = 1:nn
 #       @inbounds dx = X[nlist_j[k]] - X[nlist_i[k]]
        i = nlist_i[k]
        j = nlist_j[k]
        xi = X[i]
        xj = X[j]
        dx = xj - xi
        # Xi = wrap_to_unit(inv_bins, xi)
        # Xj = wrap_to_unit(inv_bins, xj)
        r2 = sum(abs2, dx)
        if r2 > 0.25
            dX = wrap_to_unit(inv_bins, dx)
            dx = adjust_dx(bins, dX)
            r2 = sum(abs2, dx)
        end
        r = sqrt(r2)
        push!(nlist_r, r)
        push!(nlist_R, dx)
    end
    nlist_s = s_r_r0(nlist_i, nlist_j, Rad, nlist_r)
    set_data!(at, :nlist_r, nlist_r)
    set_data!(at, :nlist_s, nlist_s)
    set_data!(at, :nlist_R, nlist_R)
#    return mynlist(nlist_i, nlist_j, nlist_r, nlist_R, nlist_s)
end

end


text:
"""
A type for positive semi-definite matrices.

Relates to PSDMatrices.jl, and I might move this to PSDMatrices.jl in the
future, but having this here allowed for easier development.
"""

struct SquarerootMatrix{T<:Real,S<:AbstractMatrix{T},M<:AbstractMatrix{T}} <:
       AbstractMatrix{T}
    squareroot::S
    mat::M
end
SquarerootMatrix(S::AbstractMatrix{T}, mat::AbstractMatrix{T}) where {T} =
    SquarerootMatrix{T,typeof(S),typeof(mat)}(S, mat)
SquarerootMatrix(S) = SquarerootMatrix(S, S * S')

Base.Matrix(M::SquarerootMatrix) = M.mat
# TODO Maybe cache the above, into M.mat or something. But do so lazily!

Base.size(M::SquarerootMatrix) = size(M.mat)
# getindex is expensive, since each call performs a matrix multiplication
Base.getindex(M::SquarerootMatrix, I::Vararg{Int,N}) where {N} = getindex(M.mat, I...)
Base.copy(M::SquarerootMatrix) = SquarerootMatrix(
    # I want the copy to look exactly as the original one
    M.squareroot isa LinearAlgebra.Adjoint ? copy(M.squareroot')' : copy(M.squareroot),
    copy(M.mat),
)
Base.copy!(dst::SquarerootMatrix, src::SquarerootMatrix) =
    (Base.copy!(dst.squareroot, src.squareroot);
    Base.copy!(dst.mat, src.mat);
    dst)
Base.similar(M::SquarerootMatrix) = SquarerootMatrix(similar(M.squareroot), similar(M.mat))

X_A_Xt(M::SquarerootMatrix, X::AbstractMatrix) = SquarerootMatrix(X * M.squareroot)
X_A_Xt!(out::SquarerootMatrix, M::SquarerootMatrix, X::AbstractMatrix) = begin
    _matmul!(out.squareroot, X, M.squareroot)
    _matmul!(out.mat, out.squareroot, out.squareroot')
    return out
end
X_A_Xt!(out::SquarerootMatrix, M::SquarerootMatrix, D::Diagonal) = begin
    # Basically just to optimize PI*Q*PI
    out.squareroot .= D.diag .* M.squareroot
    out.mat .= D.diag .* M.mat .* D.diag'
    return out
end

Base.inv(M::SquarerootMatrix) = Base.inv(M.mat)
LinearAlgebra.diag(M::SquarerootMatrix) = LinearAlgebra.diag(M.mat)
LinearAlgebra.cholesky(M::SquarerootMatrix) = LinearAlgebra.cholesky(M.mat)
LinearAlgebra.cholesky!(M::SquarerootMatrix) = LinearAlgebra.cholesky!(M.mat)
LinearAlgebra.logdet(M::SquarerootMatrix) = LinearAlgebra.logdet(M.mat)


text:
using Test, VTKDataTypes, VTKDataIO

function write_read_3d_image()
    origin = (0.0, 0.0, 0.0)
    spacing = (1.0, 2.0, 1.0)
    _extents = (2, 3, 5)

    image = VTKImageData(origin, spacing, _extents)
    image.point_data["Point scalar data"] = rand(extents(image)...)
    image.cell_data["Cell scalar data"] = rand(cell_extents(image)...)
    image.point_data["Point vector data"] = rand(3, extents(image)...)
    image.cell_data["Cell vector data"] = rand(3, cell_extents(image)...)

    write_vtk(image, "test_3d_image") # Without extension

    _image = read_vtk("test_3d_image.vti")
    @test image == _image

    return image
end

function write_read_3d_rectilinear()
    image = write_read_3d_image()
    rectilinear = VTKRectilinearData(image)
    write_vtk(rectilinear, "test_3d_rect") # Without extension
    _rectilinear = read_vtk("test_3d_rect.vtr")
    @test rectilinear == _rectilinear
    return rectilinear
end

function write_read_3d_structured()
    rectilinear = write_read_3d_rectilinear()
    structured = VTKStructuredData(rectilinear)
    write_vtk(structured, "test_3d_struct") # Without extension
    _structured = read_vtk("test_3d_struct.vts")
    @test structured == _structured
    return structured
end

function write_read_3d_unstructured()
    structured = write_read_3d_structured()
    unstruct = VTKUnstructuredData(structured)
    write_vtk(unstruct, "test_3d_unstruct") # Without extension
    _unstruct = read_vtk("test_3d_unstruct.vtu")
    @test unstruct == _unstruct
    return unstruct
end

function write_read_3d_polydata()
    rectilinear = write_read_3d_rectilinear()
    polydata = VTKPolyData(rectilinear)
    write_vtk(polydata, "test_3d_poly") # Without extension
    #When writing 2D polydata, it is automatically changed to 3D
    _polydata = read_vtk("test_3d_poly.vtp")
    @test polydata == _polydata
    return polydata
end

function write_read_3d_multiblock()
    mb = VTKMultiblockData((write_read_3d_image(), write_read_3d_rectilinear(),
        write_read_3d_structured(), write_read_3d_unstructured(),
        write_read_3d_polydata()));
    write_vtk(mb, "test_3d_multiblock") # Without extension
    _mb = read_vtk("test_3d_multiblock.vtm")
    @test mb == _mb
    return mb
end

function write_read_3d_timeseries()
    #Making a timeseries data
    x = y = z = [0., 1., 2.];
    point_coords = (x, y, z);
    point_data = Dict{String, Array{Float64}}();
    cell_data = Dict{String, Array{Float64}}();
    a = VTKRectilinearData(point_coords, point_data, cell_data);
    a.cell_data["Color"] = reshape(rand(num_of_cells(a)), cell_extents(a));
    b = VTKPolyData(a);

    x = y = z = [3., 4., 5.];
    point_coords = (x, y, z);
    c = VTKPolyData(VTKRectilinearData(point_coords, point_data, cell_data));
    c.cell_data["Color"] = rand(num_of_cells(c));

    #Multiblock data
    m = VTKMultiblockData((b, c));
    timeseries = VTKTimeSeriesData([0.], [m]);

    # Adding new time steps
    no_of_timesteps = 5;
    timestep = 0.5;
    speed = 1;
    i = 0
    while i < no_of_timesteps
        i += 1;
        new_data = deepcopy(timeseries.data[i]);

        #Random walk of first dataset
        new_data[1].point_coords = new_data[1].point_coords +
            vcat(fill(speed*rand()-speed/2, (1,num_of_points(new_data[1]))),
            fill(speed*rand()-speed/2, (1,num_of_points(new_data[1]))),
            fill(speed*rand()-speed/2, (1,num_of_points(new_data[1]))));

        #Random walk of second dataset
        new_data[2].point_coords = new_data[2].point_coords +
            vcat(fill(speed*rand()-speed/2, (1,num_of_points(new_data[2]))),
            fill(speed*rand()-speed/2, (1,num_of_points(new_data[2]))),
            fill(speed*rand()-speed/2, (1,num_of_points(new_data[2]))));

        insert_timed_data!(timeseries, i*timestep, new_data);
    end

    #Writing and reading
    write_vtk(timeseries, "test_3d_timeseries") # Without extension
    _timeseries = read_vtk("test_3d_timeseries.pvd")
    @test timeseries == _timeseries
    return timeseries
end

write_read_3d_image();
write_read_3d_rectilinear();
write_read_3d_structured();
write_read_3d_unstructured();
write_read_3d_polydata();
write_read_3d_multiblock();
#write_read_3d_timeseries();


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

function string_dictkeys_tosymbols(d::Dict)
    d2 = Dict()
    for (k, v) in d
        d2[Symbol(k)] = v
    end
    return d2
end

function filter_dict_to_match_struct_field_names(d::Dict, s::DataType)
    f = fieldnames(s)
    d2 = Dict()
    for k in f
        if haskey(d, k)
            d2[k] = d[k]
        else
            @warn "utils.jl: dict is missing struct field $(k)!"
        end
    end
    return d2
end

"""
    array_of_array_to_2D_array(aa)

Convert Array of Array to a two dimensional array
```@example
array_of_array_to_2D_array([[1,2,3], [4, 5, 6]])
````
"""
function array_of_array_to_2D_array(aa)
    transpose([aa[i][j] for j in 1:length(aa[1]), i in 1:length(aa)])
end

Base.@kwdef struct Parameter
	 ###  Sets  ###
	 Storage::Array{String,1}      # Set B in math; new; B = "Elec","HotThermal","ColdThermal"
     TechClass::Array{String,1}    # Set C
	 DemandBin::UnitRange{Int64}   # Set E
	 FuelType::Array{String,1}     # Set F; new; F = {"NaturalGas"} for CHP and whatever fuel source is used for Boiler
     TimeStep::UnitRange{Int64}    # Set H
     TimeStepBat::UnitRange{Int64} # Set H union {0}
	 Subdivision::Array{String,1}	# Set K; new; elements = { "CapCost", "FuelBurn" }
	 Month::UnitRange{Int64} 	    # Set M
	 DemandMonthsBin::UnitRange{Int64}	# Set N
	 Ratchets::UnitRange{Int64}	   # Set R
     Seg::UnitRange{Int64}	       # Set S
	 Tech::Array{String,1}         # Set T in math
	 FuelBin::UnitRange{Int64}	   # To be removed
	 PricingTier::UnitRange{Int64}  # Set U: Pricing Tiers (proposed revision) #new
	 NMILRegime::Array{String,1}	# Set V: Net-metering Regimes
     CPPeriod::UnitRange{Int64}     # Set P: Coincident Peak Periods

	 ###  Subsets and Indexed Sets  ####
	 ElecStorage::Array{String,1}  # B^{e} \subset B: Electrical energy storage systems
	 TimeStepRatchetsMonth::Array{Array{Int64,1},1}   #  H_m: Time steps in month m
	 TimeStepRatchets::Array{Array{Int64,1},1}   #  H_r: Time steps in ratchet r
     TimeStepsWithGrid::Array{Int64,1}  # H_g: Time steps with grid connection
     TimeStepsWithoutGrid::Array{Int64,1}	 # H \setminus H_g: Time steps without grid connection
	 DemandLookbackMonths::Array{Any,1}   # M^{lb}: fixed Look back months considered for peak pricing
	 DemandLookbackRange::Int   # number of Look back months considered for peak pricing
	 SegByTechSubdivision::AxisArray # S_{kt}: System size segments from segmentation k applied to technology t
	 TechsInClass::AxisArray # T_c \subset T: Technologies that are in class c
	 TechsByFuelType::AxisArray # T_f \subset T: Technologies that burn fuel type f
	 TechsByNMILRegime::AxisArray # T_v \subset T: Technologies that may acess net-meterng regime v
	 ElectricTechs::Array{String,1}  # T^{e} \subset T: Electricity-producing technologies
	 FuelBurningTechs::Array{String,1}  # T^{f} \subset T: Fuel-burning technologies
	 TechsNoTurndown::Array{String,1}  # T^{ac} \subset T: Technologies that cannot turn down, i.e., PV and wind
     CoincidentPeakLoadTimeSteps::AxisArray # H^{cp}_m: Coincident peak time steps in month m

	 ###  Parameters and Tables supporting Indexed Sets ###
	 TechToNMILMapping::AxisArray  # Defines set T_v: Technologies that may be access net-metering regime v

	 ###  Scaling Parameters ###
	 TimeStepScaling::Float64  # \Delta: Time step scaling [h], eg. 30 minute resolution -> TimeStepScaling = 0.5

	 ###  Parameters for Costs and their Functional Forms ###
     AnnualMinCharge::Float64    # c^{amc}: Utility annual minimum charge
     MonthlyMinCharge::Float64    # c^{mmc}: Utility monthly minimum charge  (not in math; will use this in min charge calculation)
	 FixedMonthlyCharge::Float64  # c^{fmc}: Utility monthly fixed charge
	 FuelCost::AxisArray # c^{u}_{f}: Unit cost of fuel type f [$/kWh]  in math  (NEW)
	 ElecRate::Array{Float64, 2}  #   c^{g}_{uh}: Grid energy cost in energy demand tier u during time step h  (NEW)
	 OMperUnitSize::AxisArray # c^{om}_{t}: Operation and maintenance cost of technology t per unit of system size [$/kW]
     OMcostPerUnitProd::AxisArray
	 OMcostPerUnitHourPerSize::AxisArray

	 GridExportRates::AxisArray  # c^{e}_{uh}: Export rate for energy in energy pricing tier u in time step h   (NEW)
	 CapCostSlope::AxisArray   # c^{cm}_{ts}: Slope of capital cost curve for technology t in segment s
     CapCostYInt::AxisArray  # c^{cb}_{ts}: Y-Intercept of capital cost curve for technology t in segment s
     CapCostX::AxisArray    # X-value of inflection point (will be changed)
	 #For the replacement of CapCostX, see new parameters SegmentLB and SegmentUB in section "System size and fuel limit parameters"
	 DemandRates::Array{Float64, 2}  # c^{r}_{re}: Cost per unit peak demand in tier e during ratchet r
	 DemandRatesMonth::Array{Float64, 2}   # c^{rm}_{mn}: Cost per unit peak demand in tier n during month m
     CoincidentPeakRates::AxisArray   # c^{cp}_p: Cost per unit peak demand during coincident peak hours of CP period p

	 ###  Demand Parameters ###
	 ElecLoad::Array{Float64,1}  # \delta^{d}_{h}: Electrical load in time step h   [kW]
     DemandLookbackPercent::Float64    # \delta^{lp}: Demand Lookback proportion [fraction]
     MaxDemandInTier::Array{Float64,1}  # \delta^{t}_{e}: Maximum power demand in ratchet e
     MaxDemandMonthsInTier::Array{Float64,1}   # \delta^{mt}_{n}: Maximum monthly power demand in tier n
     MaxUsageInTier::Array{Float64,1}   # \delta^{tu}_{u}: Maximum monthly energy demand in tier u

	 ###  Incentive Parameters ###
	 NMILLimits::AxisArray   # i^{n}_{v}: Net metering and interconnect limits in net metering regime v [kW]
     MaxProdIncent::AxisArray      # \bar{i}_t: Upper incentive limit for technology t [$]
	 ProductionIncentiveRate::AxisArray  # i^{r}_{t}: Incentive rate for technology t [$/kWh] (NEW)
	 MaxSizeForProdIncent::AxisArray  # \bar{i}^{\sigma}_t: Maximum system size to obtain production incentive for technology t [kW]

	 ###  Technology-specific Time-series Factor Parameters ###
	 ProductionFactor::AxisArray    #f^{p}_{th}  Production factor of technology t and time step h  [unitless]  (NEW)

	 ###  Technology-specific Factor Parameters ###
	 TurbineDerate::AxisArray  # f^{d}_{t}: Derate factor for turbine technologyt [unitless]
     MinTurndown::AxisArray    # f^{td}_{t}:  Minimum turn down for technology t [unitless]
     pwf_prod_incent::AxisArray   # f^{pi}_t: Present worth factor for incentives for technology t [unitless]
	 LevelizationFactor::AxisArray    # f^{l}_{t}: Levelization factor of technology t [unitless]

	 ###  Generic Factor Parameters ###
	 pwf_om::Float64  # f^{om}: Operations and maintenance present worth factor [unitless]
     pwf_e::Float64   # f^{e}: Energy present worth factor [unitless]
	 pwf_fuel::AxisArray
	 r_tax_owner::Float64      # f^{tow}: Tax rate factor for owner [fraction]
     r_tax_offtaker::Float64   # f^{tot}: Tax rate factor for offtaker [fraction]
     pwf_owner::Float64    # Annuity with zero escalation and owner's discount rate
     pwf_offtaker::Float64 # Annuity with zero escalation and offtaker's discount rate
     pwfs_emissions_cost::Dict{String, Any} # Cost of emissions present worth factors for grid and onsite fuelburn emissions [unitless]
     pwfs_grid_emissions_lbs::Dict{String, Any} # Emissions [lbs] present worth factors for grid emissions [unitless]

	 ###  System Size and Fuel Limit Parameters ###
	 TechClassMinSize::AxisArray   #  \ubar{b}^{\sigma}_{c}: Minimum system size for technology class c [kW]
	 MaxSize::AxisArray    #  \bar{b}^{\sigma}_{t}: Maximum system size for technology t [kW]
	 SegmentMinSize::AxisArray # \ubar{b}^{\sigma s}_{tks}: Minimum system size for technology t, subdivision k, segments
	 SegmentMaxSize::AxisArray  # \bar{b}^{\sigma s}_{tks}: Maximum system size for technology t, subdivision k, segments
	 FuelLimit::AxisArray # b^{fa}_{f}: Amount of available fuel for type f [kWh]   (NEW)

	 ###  Efficiency Parameters ###
	 ChargeEfficiency::AxisArray  # \eta^{esi}_{bt}: Efficiency of charging storage system b using technology t  [fraction] (NEW)
	 GridChargeEfficiency::Float64   # \eta^{esig}: Efficiency of charging electrical storage using grid power [fraction] (NEW)
     DischargeEfficiency::AxisArray  # \eta^{eso}_{b}: Efficiency of discharging storage system b [fraction] (NEW)

	 ###  Storage Parameters ###   # \ubar{w}^{bkW}_{b}: Minimum power capacity of storage system b (needs to be indexed on b )
     StorageMinChargePcent::Float64     #  \ubar{w}^{mcp}_{b}: Minimum state of charge of storage system b
     InitSOC::Float64    # w^{i}_{b} Initial percent state of charge for storage system b
     StorageMinSizeEnergy::AxisArray     # \bar{w}^{bkWh}_{b}: Maximum energy capacity of storage system b [kWh]
     StorageMaxSizeEnergy::AxisArray     # \ubar{w}^{bkWh}_{b}: Minimum energy capacity of storage system b [kWh]
     StorageMinSizePower::AxisArray     # \bar{w}^{bkW}_{b}: Maximum power capacity of storage system b [kW]
     StorageMaxSizePower::AxisArray     # \ubar{w}^{bkW}_{b}: Minimum power capacity of storage system b [kW]
     StorageMinSOC::AxisArray     #  \ubar{w}^{mcp}_{b}: Minimum state of charge of storage system b [fraction]
     StorageInitSOC::AxisArray  #Initial state of charge of storage system b [fraction]
     StorageCanGridCharge::Bool  # Boolean for storage system [fraction]

	 ###  Fuel Burn Parameters ###
	 FuelBurnSlope::AxisArray # m^\text{fm}_{t}: Fuel burn rate slope parameter for technology t
	 FuelBurnYInt::AxisArray # m^\text{fb}_{t}: Fuel burn rate slope parameter for technology t

	 ### Not used or used for calculation of other parameters ###
	 two_party_factor::Float64
     analysis_years::Int64     # Used to calculate present worth factors maybe?
     AnnualElecLoadkWh::Float64
     CapCostSegCount::Int64    # Size of set S
     FuelBinCount::Int64       # Size of set F
     DemandBinCount ::Int64    # Size of set E
     DemandMonthsBinCount::Int64   # Size of set N
     TimeStepCount::Int64          # Size of set H
     Points::UnitRange{Int64}      # CapCostSegCount+1; this is going to be the size of set S^{c} now
     PricingTierCount::Int64    # Size of set U

     # new parameters for reformulation
     StorageCostPerKW::AxisArray
     StorageCostPerKWH::AxisArray
	 ElectricDerate::AxisArray

     # New Sets
     ExportTiers::Array{String,1}
	 ExportTiersByTech::AxisArray
	 TechsByExportTier::AxisArray
     ExportTiersBeyondSiteLoad::Array{String, 1}
     TechsCannotCurtail::Array{String, 1}

    # Feature Additions
     TechToLocation::AxisArray
     MaxSizesLocation::Array{Float64, 1}
     Location::UnitRange
	 AddSOCIncentive::Int64

    # Annual RE parameters
     TechPercentRE::AxisArray
     MinAnnualPercentREElec::Union{Float64,Nothing}
     MaxAnnualPercentREElec::Union{Float64,Nothing}
     IncludeExportedREElecinTotal::Bool

    # Emissions parameters
     IncludeExportedElecEmissionsInTotal::Bool
     MinPercentCO2EmissionsReduction::Union{Float64,Nothing}
     MaxPercentCO2EmissionsReduction::Union{Float64,Nothing}
     BAUYr1Emissions_CO2::Float64
     BAUYr1Emissions_NOx::Float64
     BAUYr1Emissions_SO2::Float64
     BAUYr1Emissions_PM25::Float64
     BAUYr1Emissions_grid_CO2::Float64
     BAUYr1Emissions_grid_NOx::Float64
     BAUYr1Emissions_grid_SO2::Float64
     BAUYr1Emissions_grid_PM25::Float64
     GridEmissionsFactor_CO2::Array{Float64,1}
     GridEmissionsFactor_NOx::Array{Float64,1}
     GridEmissionsFactor_SO2::Array{Float64,1}
     GridEmissionsFactor_PM25::Array{Float64,1}
     TechEmissionsFactors_CO2::AxisArray
     TechEmissionsFactors_NOx::AxisArray
     TechEmissionsFactors_SO2::AxisArray
     TechEmissionsFactors_PM25::AxisArray
     CO2_dollars_tonne::Float64
     NOx_dollars_tonne_grid::Float64
     SO2_dollars_tonne_grid::Float64
     PM25_dollars_tonne_grid::Float64
     NOx_dollars_tonne_onsite_fuelburn::Float64
     SO2_dollars_tonne_onsite_fuelburn::Float64
     PM25_dollars_tonne_onsite_fuelburn::Float64
     Include_climate_in_objective::Bool
     Include_health_in_objective::Bool
     Lbs_per_tonne::Float64

	# Added for CHP
	HotTES::Array{String,1}
	ColdTES::Array{String,1}
	ThermalStorage::Array{String,1}
	CHPTechs::Array{String,1}
	ElectricChillers::Array{String,1}
	AbsorptionChillers::Array{String,1}
	CoolingTechs::Array{String,1}
	HeatingTechs::Array{String,1}
	BoilerTechs::Array{String,1}
	HeatingLoad::Array{Float64,1}
	CoolingLoad::Array{Float64,1}
	BoilerEfficiency::AxisArray
	ElectricChillerCOP::Float64
    AbsorptionChillerCOP::Float64
    AbsorptionChillerElecCOP::Float64
	CHPThermalProdSlope::AxisArray
	CHPThermalProdIntercept::AxisArray
	FuelBurnYIntRate::AxisArray
	CHPThermalProdFactor::AxisArray
	CHPStandbyCharge::Float64
	CHPDoesNotReduceDemandCharges::Int64
	StorageDecayRate::AxisArray
    AllBoilerTechs::Array{String,1}
    AllTechsForSteamTurbine::Array{String,1}
    SteamTurbineTechs::Array{String,1}
    TechCanSupplySteamTurbine::Array{String,1}
    STElecOutToThermInRatio::Float64
    STThermOutToThermInRatio::Float64
    # GHP Arrays of different GHP options with index 1 being NO GHP
    GHPOptions::UnitRange{Int64}
    RequireGHPPurchase::Int64
    GHPHeatingThermalServed::Array{Float64,2}  # Array of heating load (thermal!) profiles served by GHP
    GHPCoolingThermalServed::Array{Float64,2}  # Array of cooling load profiles served by GHP
    HeatingThermalReductionWithGHP::Array{Float64,2}  # Array of heating load reduction (thermal!) profile from GHP retrofit
    CoolingThermalReductionWithGHP::Array{Float64,2}  # Array of cooling load reduction (thermal!) profile from GHP retrofit
    GHPElectricConsumed::Array{Float64,2}  # Array of electric load profiles consumed by GHP
    GHPInstalledCost::Array{Float64,1}  # Array of installed cost for GHP options
    GHPOMCost::Array{Float64,1}  # Array of O&M cost for GHP options
    CHPSupplementaryFireMaxRatio::Float64
    CHPSupplementaryFireEfficiency::Float64
    CapCostSupplementaryFiring::AxisArray # Array of capital cost for supplementary firing for CHP
	#Offgrid systems
	OffGridFlag::Bool
	TechsRequiringSR::Array{String,1}
	TechsProvidingSR::Array{String,1}
	MinLoadMetPct::Float64
	SRrequiredPctLoad::Float64
	SRrequiredPctTechs::AxisArray
    OtherCapitalCosts::Float64
    OtherAnnualCosts::Float64
end

function Parameter(d::Dict)
    can_be_empty = (
        "MaxSize",
        "OMperUnitSize",
        "OMcostPerUnitProd",
		"OMcostPerUnitHourPerSize",
        "MaxProdIncent",
        "MaxSizeForProdIncent",
        "TurbineDerate",
        "MinTurndown",
        "pwf_prod_incent",
        "LevelizationFactor",
        "NMILLimits",
        "TechClassMinSize",
		"TechsByExportTier",
		"ExportTiersByTech",
		"NMILRegime",
		"TechsByNMILRegime",
		"TechsByFuelType",
		"FuelCost"
     )
	for x in ["Tech","FuelType","CHPTechs"]
		if typeof(d[x]) === Array{Any, 1}  # came from Python as empty array
			d[x] = convert(Array{String, 1}, d[x])
		end
	end
    for x in can_be_empty
        if typeof(x) === Array{Any, 1}  # came from Python as empty array
            d[x] = convert(Array{Float64, 1}, d[x])
        end
    end
    d = convert(Dict{Union{String, Symbol}, Any}, d)

    d[:Seg] = 1:d["CapCostSegCount"]
    d[:Points] = 0:d["CapCostSegCount"]
    d[:Month] = 1:12
    d[:Ratchets] = 1:d["NumRatchets"]
    d[:FuelBin] = 1:d["FuelBinCount"]
	d[:PricingTier] = 1:d["PricingTierCount"]
    d[:DemandBin] = 1:d["DemandBinCount"]
    d[:DemandMonthsBin] = 1:d["DemandMonthsBinCount"]
    d[:TimeStep] = 1:d["TimeStepCount"]
    d[:TimeStepBat] = 0:d["TimeStepCount"]
    n_location = length(d["MaxSizesLocation"])
    d[:Location] = 1:n_location
    d[:CPPeriod] = 1:d["CoincidentPeakPeriodCount"]
    d[:GHPOptions] = 1:length(d["GHPInstalledCost"])

    # the following array manipulation may have to adapt once length(d["Subdivision"]) > 1
    seg_min_size_array = reshape(transpose(reshape(d["SegmentMinSize"], length(d[:Seg]), length(d["Tech"]))),
                                 length(d["Tech"]), length(d["Subdivision"]), length(d[:Seg]))
    seg_max_size_array = reshape(transpose(reshape(d["SegmentMaxSize"], length(d[:Seg]), length(d["Tech"]))),
                                length(d["Tech"]), length(d["Subdivision"]), length(d[:Seg]))

    # convert vectors to AxisArray's with axes for REopt JuMP model
    d["TurbineDerate"] = AxisArray(d["TurbineDerate"], d["Tech"])
    d["TechToLocation"] = vector_to_axisarray(d["TechToLocation"], d["Tech"], d[:Location])
    d["pwf_prod_incent"] = AxisArray(d["pwf_prod_incent"], d["Tech"])
    d["LevelizationFactor"] = AxisArray(d["LevelizationFactor"], d["Tech"])
    d["OMperUnitSize"] = AxisArray(d["OMperUnitSize"], d["Tech"])
    d["CapCostSlope"] = vector_to_axisarray(d["CapCostSlope"], d["Tech"], d[:Seg])
    d["CapCostYInt"] = vector_to_axisarray(d["CapCostYInt"], d["Tech"], d[:Seg])
    d["CapCostX"] = vector_to_axisarray(d["CapCostX"], d["Tech"], d[:Points])
    d["MaxProdIncent"] = AxisArray(d["MaxProdIncent"], d["Tech"])
    d["MaxSizeForProdIncent"] = AxisArray(d["MaxSizeForProdIncent"], d["Tech"])
    d["MaxSize"] = AxisArray(d["MaxSize"], d["Tech"])
    d["TechClassMinSize"] = AxisArray(d["TechClassMinSize"], d["TechClass"])
    d["MinTurndown"] = AxisArray(d["MinTurndown"], d["Tech"])
    d["DemandRates"] = transpose(reshape(d["DemandRates"], d["DemandBinCount"], d["NumRatchets"]))
    d["DemandRatesMonth"] = transpose(reshape(d["DemandRatesMonth"], d["DemandMonthsBinCount"], 12))
    d["NMILLimits"] = AxisArray(d["NMILLimits"], d["NMILRegime"])
    d["TechToNMILMapping"] = vector_to_axisarray(d["TechToNMILMapping"], d["Tech"], d["NMILRegime"])
    d["OMcostPerUnitProd"] = AxisArray(d["OMcostPerUnitProd"], d["Tech"])
    d["OMcostPerUnitHourPerSize"] = AxisArray(d["OMcostPerUnitHourPerSize"], d["Tech"])

    d["TechEmissionsFactors_CO2"] = AxisArray(d["TechEmissionsFactors_CO2"], d["Tech"])
    d["TechEmissionsFactors_NOx"] = AxisArray(d["TechEmissionsFactors_NOx"], d["Tech"])
    d["TechEmissionsFactors_SO2"] = AxisArray(d["TechEmissionsFactors_SO2"], d["Tech"])
    d["TechEmissionsFactors_PM25"] = AxisArray(d["TechEmissionsFactors_PM25"], d["Tech"])
    d["TechPercentRE"] = AxisArray(d["TechPercentRE"], d["Tech"])

    if !isempty(d["CoincidentPeakLoadTimeSteps"])
        d["CoincidentPeakRates"] = AxisArray(d["CoincidentPeakRates"], d[:CPPeriod])
        d["CoincidentPeakLoadTimeSteps"] = permutedims(hcat(d["CoincidentPeakLoadTimeSteps"]...), (2,1))
        d["CoincidentPeakLoadTimeSteps"] = AxisArray(d["CoincidentPeakLoadTimeSteps"], d[:CPPeriod], 1:size(d["CoincidentPeakLoadTimeSteps"],2))
    else
        d["CoincidentPeakRates"] = AxisArray([])
        d["CoincidentPeakLoadTimeSteps"] = AxisArray([])
    end

    # Reformulation additions
    d["StorageCostPerKW"] = AxisArray(d["StorageCostPerKW"], d["Storage"])
    d["StorageCostPerKWH"] = AxisArray(d["StorageCostPerKWH"], d["Storage"])
	d["FuelCost"] = vector_to_axisarray(d["FuelCost"], d["FuelType"], d[:TimeStep])
    d["ElecRate"] = transpose(reshape(d["ElecRate"], d["TimeStepCount"], d["PricingTierCount"]))

    if !isempty(d["GridExportRates"])
        d["GridExportRates"] = AxisArray(array_of_array_to_2D_array(d["GridExportRates"]),
                                         d["ExportTiers"], d[:TimeStep])
    else
        d["GridExportRates"] = AxisArray([])
    end
    d["FuelBurnSlope"] = AxisArray(d["FuelBurnSlope"], d["FuelBurningTechs"])
    d["FuelBurnYInt"] = AxisArray(d["FuelBurnYInt"], d["FuelBurningTechs"])
    d["ProductionFactor"] = vector_to_axisarray(d["ProductionFactor"], d["Tech"], d[:TimeStep])
    d["ProductionIncentiveRate"] = AxisArray(d["ProductionIncentiveRate"], d["Tech"])
    d["FuelLimit"] = AxisArray(d["FuelLimit"], d["FuelType"])
    d["ChargeEfficiency"] = vector_to_axisarray(d["ChargeEfficiency"], d["Tech"], d["Storage"])
    d["DischargeEfficiency"] = AxisArray(d["DischargeEfficiency"], d["Storage"])
    d["StorageMinSizeEnergy"] = AxisArray(d["StorageMinSizeEnergy"], d["Storage"])
    d["StorageMaxSizeEnergy"] = AxisArray(d["StorageMaxSizeEnergy"], d["Storage"])
    d["StorageMinSizePower"] = AxisArray(d["StorageMinSizePower"], d["Storage"])
    d["StorageMaxSizePower"] = AxisArray(d["StorageMaxSizePower"], d["Storage"])
    d["StorageMinSOC"] = AxisArray(d["StorageMinSOC"], d["Storage"])
    d["StorageInitSOC"] = AxisArray(d["StorageInitSOC"], d["Storage"])

    d["SegmentMinSize"] = AxisArray(seg_min_size_array, d["Tech"], d["Subdivision"], d[:Seg])
    d["SegmentMaxSize"] = AxisArray(seg_max_size_array, d["Tech"], d["Subdivision"], d[:Seg])
	d["ElectricDerate"] = vector_to_axisarray(d["ElectricDerate"], d["Tech"], d[:TimeStep])

	# CHP Additions
	d["CHPThermalProdSlope"] = AxisArray(d["CHPThermalProdSlope"],d["CHPTechs"])
	d["CHPThermalProdIntercept"] = AxisArray(d["CHPThermalProdIntercept"],d["CHPTechs"])
	d["FuelBurnYIntRate"] = AxisArray(d["FuelBurnYIntRate"],d["CHPTechs"])
	d["CHPThermalProdFactor"] = vector_to_axisarray(d["CHPThermalProdFactor"],d["CHPTechs"],d[:TimeStep])
	d["pwf_fuel"] = AxisArray(d["pwf_fuel"], d["Tech"])
	d["StorageDecayRate"] = AxisArray(d["StorageDecayRate"], d["Storage"])
    d["BoilerEfficiency"] = AxisArray(d["BoilerEfficiency"], d["AllBoilerTechs"])  # Always passes both values, even if partial/none
    d["CapCostSupplementaryFiring"] = AxisArray(d["CapCostSupplementaryFiring"],d["CHPTechs"])

	# Off-grid Modeling
	d["SRrequiredPctTechs"] = AxisArray(d["SRrequiredPctTechs"], d["TechsProvidingSR"])

    # Indexed Sets
    if isempty(d["FuelType"])
        d["TechsByFuelType"] = []  # array of arrays is not empty, but needs to be for AxisArray conversion
    end
    d["SegByTechSubdivision"] = vector_to_axisarray(d["SegByTechSubdivision"], d["Subdivision"], d["Tech"])
    d["TechsByFuelType"] = AxisArray(d["TechsByFuelType"], d["FuelType"])
    d["TechsInClass"] = AxisArray(d["TechsInClass"], d["TechClass"])
    d["ExportTiersByTech"] = AxisArray(d["ExportTiersByTech"], d["Tech"])
	d["TechsByExportTier"] = AxisArray(d["TechsByExportTier"], d["ExportTiers"])
    d["TechsByNMILRegime"] = AxisArray(d["TechsByNMILRegime"], d["NMILRegime"])

    d["GHPHeatingThermalServed"] = array_of_array_to_2D_array(d["GHPHeatingThermalServed"])
    d["GHPCoolingThermalServed"] = array_of_array_to_2D_array(d["GHPCoolingThermalServed"])
    d["HeatingThermalReductionWithGHP"] = array_of_array_to_2D_array(d["HeatingThermalReductionWithGHP"])
    d["CoolingThermalReductionWithGHP"] = array_of_array_to_2D_array(d["CoolingThermalReductionWithGHP"])
    d["GHPElectricConsumed"] = array_of_array_to_2D_array(d["GHPElectricConsumed"])

    d = string_dictkeys_tosymbols(d)
    d = filter_dict_to_match_struct_field_names(d, Parameter)

    param = Parameter(;d...)
end

# Additional dispatches to make things easier
function length(::Symbol)
    return 1
end

JuMP.value(::Val{false}) = 0

JuMP.value(x::Number) = x  # for @expression(REopt, ExportBenefitYr1, 0.0) and similar

JuMP.value(::MutableArithmetics.Zero) = 0

function vector_to_axisarray(v::Array{<:Any, 1}, ax1::Array{<:Any, 1}, ax2::Union{UnitRange, Array{<:Any, 1}})
    l1 = length(ax1)
    l2 = length(ax2)
    a = transpose(reshape(v, l2, l1))
    return AxisArray(a, ax1, ax2)
end


text:
"""
Native binary RData format IO stream wrapper.

TODO write readers
"""
struct NativeIO{T<:IO} <: RDAIO
    sub::T               # underlying IO stream

    NativeIO(io::T) where {T<:IO} = new{T}(io)
end


text:
include(joinpath(@__DIR__, "FLS.jl"))
using Base.Iterators: product as cartprod
using Catlab.CategoricalAlgebra.FinSets: preimage

"""
A key data structure in this file is a partial model. This is replacing a C-Set
with a C-Rel; we allow a morphism to map from one element to multiple targets.
The model is COMPLETE when each element has exactly one arrow per morphism.
It is INCONSISTENT if any element has zero arrows. A model we know nothing about
has a complete digraph for each morphism. Partial models are represented as
CSets with a particular schema (morphisms are replaced with relations, i.e.
objects a pair of morphisms: FK_src and FK_tgt).
"""

"""
For p1=p2, if p1 is completely known but p2 is only known up to penultimate
value, then we learn new info.
"""
function propagate_diagram!(fls::LabeledFLS, cset::CSet)::Nothing
  for sym in fls.labels[1]
    propagate_diagram!(fls, cset, sym)
  end
  return nothing
end

"""
Use diagrams to rule out certain possibilities
"""
function propagate_diagram!(fls::LabeledFLS, cset::CSet, sym::Symbol)::Nothing
  paths = all_paths(fls, sym)
  verts = sort(collect(Set(values(paths))))
  vertsyms = [fls.labels[1][fls.fls[:dvMap][v]] for v in verts]
  to_delete = Dict{Symbol, Set{Int}}([l=>Set{Int}() for l in fls.labels[2]])
  vals = Dict{Vector{Symbol},Vector{Set{Int}}}(
      [Symbol[] => Set{Int}[Set{Int}([i]) for i in 1:nparts(cset, sym)]])
  for (p, _) in filter(x->!isempty(x[1]), paths)
    tail, head = p[1:end-1], p[end]
    head_src, head_tgt = add_srctgt(head)
    valsp = Set{Int}[union([Set{Int}(
      cset[head_tgt][cset.indices[head_src][i]]) for i in iset]...)
           for iset in vals[tail]]
    vals[p] = valsp
  end
  # for each junction, check if any penultimate node is completely determined
  junctions = filter(x-> count(==(x[1]), values(paths)) > 1,
                     collect(zip(verts, vertsyms)))
  for (j, jsym) in junctions
    jpaths = [p for (p, v) in collect(paths) if v == j]
    for i in 1:nparts(cset, jsym)  # consider paths starting at each elem
      alljvals = [vals[jpath][i] for jpath in jpaths]
      jvals = union(Set{Int}(), [jv for jv in alljvals if length(jv)==1]...)
      @assert length(jvals) in [0, 1]
      if !isempty(jvals) # we know what value to set this path equal to
        jval = pop!(jvals)
        # can only use this info if the PENULTIMATE path is unambiguous
        for jpath in filter(!isempty, jpaths)
          tail, head = jpath[1:end-1], jpath[end]
          head_src, head_tgt = add_srctgt(head)
          if length(vals[tail][i])==1
            penult_val = collect(vals[tail][i])[1]
            # keep just the FK that goes to jval
            all_fks = cset.indices[head_src][penult_val]
            bad_fks = [fk for fk in all_fks if cset[head_tgt][fk] != jval]
            union!(to_delete[head], Set(bad_fks))
          end
        end
      end
    end
  end
  rem!(cset, to_delete)
end

"""Get the indices that satisfy equations starting at a particular node"""
function check_diagram(fls::LabeledFLS, cset::CSet, sym::Symbol;
                       cset_is_crel::Bool=false)::Set{Int}
  cset_ = cset_is_crel ? cset : cset_to_crels(fls, cset)
  dia = dia_to_crel(fls, get_diagram(fls, sym))
  hs = homomorphisms(dia, cset_)
  sat_inds = [h.components[sym](1) for h in hs]
  return Set(sat_inds)
end

"""
Coerce a LABELED GRAPH diagram into a C-Rel. The existence of a homomorphism
from this CRel into a partial model CRel indicates that the diagram is
satisfiable for some choice of FKs
"""
function dia_to_crel(fls::LabeledFLS, dia::LabeledGraph)::CSet
  res = crels(fls)()
  ref = Dict{Pair{Symbol, Int}, Int}()
  for (i, lab) in enumerate(dia[:vlabel])
    ind = add_part!(res, lab)
    ref[lab => i]=ind
  end
  for (e, ss, tt) in zip(dia[:elabel], dia[:src], dia[:tgt])
    slab, tlab = dia[:vlabel][[ss,tt]]
    s, t = add_srctgt(e)
    args = Dict([s => ref[slab => ss], t=> ref[tlab => tt]])
    add_part!(res, e; args...)
  end
  return res
end

"""
Take a real CSet and construct CRel with exactly one FK value per PK.
"""
function cset_to_crels(fls::LabeledFLS, cset::CSet)::CSet
  res = crels(fls)()
  for ob in fls.labels[1]
    add_parts!(res, ob, nparts(cset, ob))
  end
  for (i, ss) in enumerate(fls.fls[:src])
    isym = fls.labels[2][i]
    s, t = add_srctgt(isym)
    for src_ind in 1:nparts(cset, fls.labels[1][ss])
      sym = string(isym)
      tval = length(sym)>3 && sym[1:4] == "_id_" ? src_ind : cset[isym][src_ind]
      args = Dict([s=>src_ind, t=>tval])
      add_part!(res, isym; args...)
    end
  end
  return res
end

"""
Cone constraints can allow us to infer things about a partial model
"""
function propagate_cone!(fls::LabeledFLS, cset::CSet)::Nothing
  return nothing # to do
end

"""
Return false if unsat

Demand that every apex has at least one base (# of bases could
potentially decrease as we remove FKs).

If the legs have been finalized, then we demand every apex has
exactly one base

For apexes that have been finalized (have a determinate base)
then we check that these are unique.
"""
function check_cone(fls::LabeledFLS, cset::CSet)::Bool
  for sym in fls.labels[1]
    if !check_cone(fls, cset, sym)
      return false
    end
  end
  return true
end

"""
Check if every edge in a diagram is fully determined in a partial model
"""
function check_base_final(fls::LabeledFLS, crel::CSet, base::LabeledGraph)::Bool
  ns = n_fks(fls, crel)
  for e in Set(base[:elabel])
    for n in ns[e]
      if n != 1
        return false
      end
    end
  end
  return true
end

"""
Return true if a cone is satisfiable by a partial model
"""
function check_cone(fls::LabeledFLS, cset::CSet, apx::Symbol)::Bool
  cone = get_cone(fls, apx, false)
  base = get_cone(fls, apx, true)
  n = nparts(cset, apx)
  if base === nothing
    return true # no cone
  end

  basefinal = check_base_final(fls, cset, base)
  basedia = dia_to_crel(fls, base)
  matches = [h.components for h in homomorphisms(basedia, cset)]
  if basefinal && length(matches) != n
    # println("base is final but $(length(matches)) != # of apexes ($n)")
    return false
  end
  apexes = [get_apexes(cset, apx, cone, h) for h in matches]
  all_apexes = union(Set{Int}(), apexes...)
  cover = length(all_apexes) == n
  if !cover
    return false  # there are elements in the apex component w/ no base
  end
  seen_base = Set()
  for i in 1:n
    known_base = []
    known = true
    for e in cone.indices[:src][nv(cone)]
      esrc, etgt = add_srctgt(cone[:elabel][e])
      leg = cset[etgt][cset.indices[esrc][i]]
      if length(leg) == 1
        push!(known_base, leg[1])
      else
        known = false
        break
      end
    end
    if known
      if known_base in seen_base
        # println("Duplicate base")
        return false
      else
        push!(seen_base, known_base)
      end
    end
  end
  return true
end

"""
Return potential apexes for a given base
"""
function get_apexes(cset::CSet, apx::Symbol, cone::ACSet, morph)::Set{Int}
  poss_apex_inds = Set(1:nparts(cset, apx))
  for e in cone.indices[:src][nv(cone)]
    esrc, etgt = add_srctgt(cone[:elabel][e])
    tgt = cone[:tgt][e]
    tgt_sym = cone[:vlabel][tgt]
    tgt_ind = findfirst(==(tgt), [i for (i,s) in enumerate(cone[:vlabel])
                                  if s == tgt_sym])
    tgt_in_cset = morph[tgt_sym](tgt_ind)
    intersect!(poss_apex_inds, cset[esrc][cset.indices[etgt][tgt_in_cset]])
  end
  return poss_apex_inds
end

"""CSet type of C-Rel"""
function crels(fls::LabeledFLS)::Type
  pres = Presentation(FreeSchema)
  nv = length(fls.labels[1])
  alledge = Symbol[]
  xobs = [Ob(FreeSchema, s) for s in vcat(collect(fls.labels)...)]
  for x in xobs
    add_generator!(pres, x)
  end
  for (i,(src, tgt)) in enumerate(zip(fls.fls[:src], fls.fls[:tgt]))
    s, t = add_srctgt(fls.labels[2][i])
    add_generator!(pres, Hom(s, xobs[nv+i], xobs[src]))
    add_generator!(pres, Hom(t, xobs[nv+i], xobs[tgt]))
    append!(alledge, [s,t])
  end
  return CSetType(pres, index=alledge)
end

"""Cset type associated with FLS graph"""
function fls_csettype(fls::LabeledFLS)::Type
  pres = Presentation(FreeSchema)
  xobs = [Ob(FreeSchema, s) for s in fls.labels[1]]
  for x in xobs
    add_generator!(pres, x)
  end
  for (isym,s, t) in zip(fls.labels[2], fls.fls[:src], fls.fls[:tgt])
    add_generator!(pres, Hom(isym, xobs[s], xobs[t]))
  end
  return CSetType(pres, index=fls.labels[2])
end

"""
Convert from CRel to CSet by making an arbitrary choice for each relation.
Result is unique if the partial model is full determined, and it will fail if
the model is inconsistent.
"""
function crels_to_cset(fls::LabeledFLS, crels::CSet)::CSet
  res = fls_csettype(fls)()
  for v in fls.labels[1]
    add_parts!(res, v, nparts(crels, v))
  end
  for e in fls.labels[2]
    esrc, etgt = add_srctgt(e)
    evec = crels[etgt][[x[1] for x in crels.indices[esrc]]]
    set_subpart!(res, e, evec)
  end
  return res
end

"""Yes or no, does a CSet/CRel satisfy the diagrams of the FLS"""
function check_diagrams(fls::LabeledFLS, cset::CSet, error::Bool=true;
            cset_is_crel::Bool=false)::Bool
  for s in fls.labels[1]
    sat_inds = check_diagram(fls, cset, s,
                             cset_is_crel=cset_is_crel)
    if length(sat_inds) != nparts(cset, s)
      bad = sort(collect(setdiff(1:nparts(cset,s), sat_inds)))
      if error
        @assert false "Equations for $s not satisfied by indices $bad"
      else
        return false
      end
    end
  end
  return true
end

"""
Given a CSet from a FreeSchema presentation (e.g. Reflexive Graph) and a CSet
(e.g. an arbitrary instance of Reflexive Graph) throw an error if the CSet
does not conform to the equations in the presentation.

Once equational information is stored in the CatDesc of the CSet, we will no
longer need to pass a presentation as well.
"""
function check_diagrams(p::Presentation, cset::CSet)::Nothing
  fls = catpres_to_linear(p)
  crel = cset_to_crels(fls, cset)
  check_diagrams(fls, crel, true)
  return nothing
end

"""Populate crel with total relations: all possibilities"""
function allposs(fls::LabeledFLS, consts::Vector{Int})::CSet
  crel = crels(fls)()
  for (i, c) in enumerate(consts)
    add_parts!(crel, fls.labels[1][i], c)
  end
  for (i, (s, t)) in enumerate(zip(fls.fls[:src], fls.fls[:tgt]))
    fs, ft = add_srctgt(fls.labels[2][i])
    for j in 1:consts[s]
      for k in 1:consts[t]
        args = Dict([fs => j, ft => k])
        add_part!(crel, fls.labels[2][i]; args...)
      end
    end
  end
  return crel
end

"""
Naive search, needs to be improved.
We propagate info (or fail) based on diagram and cone constraints. Then branch
on an arbitrary choice.
"""
function search(fls::LabeledFLS,
                consts::Vector{Int},
                n::Int=-1,
                guess::Union{Nothing,CSet}=nothing;
               )::Vector{CSet}
  # initialize
  @assert all(x->x>=0, consts)
  @assert length(consts) == length(fls.labels[1])
  seen, res = Set(), Set()
  init = guess===nothing ? allposs(fls, consts) : guess

  function search_rec(crel::CSet, depth::Int)::Nothing
    if length(res) == n
      return nothing
    end

    if !check_diagrams(fls, crel, false; cset_is_crel=true)
      return nothing
    end
    if !check_cone(fls, crel)
      return nothing
    end
    old = deepcopy(crel)
    propagate_diagram!(fls, crel)
    if crel!=old
      hsh = canonical_hash(crel)
      if hsh in seen
        return nothing
      else
        push!(seen, hsh)
        old = deepcopy(crel)
      end
    end
    propagate_cone!(fls, crel)
    if crel!=old
      hsh = canonical_hash(crel)
      if hsh in seen
        return nothing
      else
        push!(seen, hsh)
      end
    end

    lens = filter(x->!isempty(x[2]), n_fks(fls, crel))
    maxes = [k=>findmax(v) for (k, v) in collect(lens)]
    mins = [minimum(v) for (_, v) in collect(lens)]
    nmax = maximum([x[1] for (_, x) in maxes])

    if minimum(mins) < 1
      return nothing
    elseif nmax == 1
      push!(res, crel)
    else
      emax, maxind = [(k, v[2]) for (k, v) in maxes if v[1] == nmax][1]
      for i in 1:nmax # split on the largest FK
        newer_crel = choose(crel, maxind, emax, i)
        newhsh = canonical_hash(newer_crel)
        if !(newhsh in seen)
          push!(seen, newhsh)
          search_rec(newer_crel, depth+1)
        end
      end
    end
  return nothing
  end

  search_rec(init, 0)
  return collect(res)
end

"""
For each PK+FK combo, how many FKs have that source
"""
function n_fks(fls::LabeledFLS, crel::CSet)::Dict{Symbol, Vector{Int}}
  res = Dict{Symbol, Vector{Int}}()
  for edge_sym in fls.labels[2]
    fks = map(length, crel.indices[add_srctgt(edge_sym, true)])
    res[edge_sym] = fks
  end
  return res
end

"""
Remove elements from tables - does not change any fk indices
so is only valid for tables which have no arrows to them
OR removing the LAST indices of a table
"""
function rem!(crel::CSet, rem::Dict{Symbol, Set{Int}})::Nothing
  for (comp, reminds) in collect(rem)
    rem_parts!(crel, comp, sort(collect(reminds)))
  end
  return nothing
end

"""Single table version of `rem`"""
function rem!(crel::CSet, tab::Symbol, rem::Set{Int})::Nothing
  rem_parts!(crel, tab, sort(collect(rem)))
  return nothing
end

"""
Eliminate options from a partial model by choosing a particular
FK combination for a specified PK.
"""
function choose(crel::CSet{CD}, srcind::Int, fk::Symbol, choice::Int
         )::CSet{CD} where {CD}
  fksrc = add_srctgt(fk, true)
  fks = crel.indices[fksrc][srcind]
  @assert 0 < choice <= length(fks)
  chosen = fks[choice]
  delinds = Set([fk for fk in fks if fk != chosen])
  res = deepcopy(crel)
  rem!(res, fk, delinds)
  return res
end

"""
1. Enumerate elements of ℕᵏ for an underlying graph with k nodes.
2. For each of these: (c₁, ..., cₖ) create a term model with that many constants

Do the first enumeration by incrementing n_nonzero and finding partitions so
that ∑(c₁,...) = n_nonzero

In the future, this function will write results to a database
that hashes the FLS as well as the set of constants that generated the model.

Also crucial is to decompose FLS into subparts that can be efficiently solved
and have solutions stitched together.
"""
function enum(fls::LabeledFLS, n::Int)::Dict{Tuple,Vector{CSet}}
    res = Dict{Tuple,CSet}()
    n_const = 0 # total number of constants across all sets
    n_v = nv(fls.G)
    while length(res) < n
        n_const += 1
        for n_nonzero in 1:n_v
            # values we'll assign to nodes
            c_parts = partitions(n_const, n_nonzero)
            # Which nodes we'll assign them to
            indices = permutations(1:n_v,n_nonzero)

            for c_part in c_parts
                for index_assign in indices
                    n_const = zeros(n_v)
                    n_const[index_assign] = c_part
                    n_c = tuple(n_const...)
                    if !(n_c in res)
                        res[n_c] =  term_models(flv, n_c)
                    end
                end
            end
        end
    end
end


text:
using Femshop, BenchmarkTools

SUITE = BenchmarkGroup()

include("rawcode.jl")


text:
# utils
export variables, operator, affine
import Base: convert
import AbstractOperators: displacement

convert(::Type{Expression},x::Variable{T,N,A}) where {T,N,A} =
Expression{1}((x,),Eye(T,size(x)))

"""
    variables(ex::Expression)

Returns a tuple containing the `Variable`s of expression `ex`.

# Example

```julia
julia> x,y = Variable(2),Variable(2);

julia> ex = x+y;

julia> variables(ex)
(Variable(Float64, (2,)), Variable(Float64, (2,)))

```

"""
variables(A::Expression)    = A.x
variables(x::Variable)    = x

"""
    operator(ex::Expression)

Returns the `AbstractOperator` of expression `ex`.

# Example

```julia
julia> x = Variable(3)
Variable(Float64, (3,))

julia> ex = fft(x);

julia> operator(ex)
ℱ  ℝ^3 -> ℂ^3

```
"""
operator(A::Expression) = remove_displacement(A.L)
operator(x::Variable)   = Eye(~x)

"""
    affine(ex::Expression)

Returns the `AbstractOperator` of expression `ex` keeping any affine addition.

"""
affine(A::Expression) = A.L
affine(x::Variable)   = Eye(~x)

"""
    displacement(ex::Expression)

Returns the displacement of expression `ex`.

# Example

```julia
julia> x = Variable(3)
Variable(Float64, (3,))

julia> ex = fft(x)+[1.+im*2.;0.;3.+im*4];

julia> displacement(ex)
3-element Array{Complex{Float64},1}:
1.0+2.0im
0.0+0.0im
3.0+4.0im

```
"""
displacement(A::Variable{T}) where {T} = zero(T)
displacement(E::Expression) = displacement(affine(E))


text:
function NewtonRaphsonSolver()
    return NewtonRaphsonSolver(())
end

function NewtonRaphsonSolver(arg0::jdouble)
    return NewtonRaphsonSolver((jdouble,), arg0)
end

function solve(obj::NewtonRaphsonSolver, arg0::jint, arg1::UnivariateDifferentiableFunction, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateDifferentiableFunction, jdouble, jdouble), arg0, arg1, arg2, arg3)
end


text:
using LinearAlgebra: det

using ..Inputs: Ibrav

"""
    Bravais(nml::SystemNamelist)

Return a `Bravais` from a `SystemNamelist`.
"""
Bravais(nml::SystemNamelist) = Bravais(Ibrav(nml.ibrav))

"""
    Lattice(nml::SystemNamelist)

Return a `Lattice` from a `SystemNamelist`.
"""
Lattice(nml::SystemNamelist) = Lattice(Bravais(nml), nml.celldm)

struct InformationNotEnough <: Exception
    msg::AbstractString
end

"""
    cellvolume(card)

Return the cell volume of a `CellParametersCard` or `RefCellParametersCard`, in atomic unit.

!!! warning
    It will throw an error if the option is `"alat"`.
"""
function cellvolume(card::AbstractCellParametersCard)
    option = optionof(card)
    if option == "bohr"
        return abs(det(card.data))
    elseif option == "angstrom"
        return ustrip(u"bohr^3", abs(det(card.data)) * u"angstrom^3")
    else  # option == "alat"
        throw(InformationNotEnough("parameter `celldm[1]` needed!"))
    end
end
"""
    cellvolume(nml::SystemNamelist)

Return the volume of the cell based on the information given in a `SystemNamelist`, in atomic unit.
"""
cellvolume(nml::SystemNamelist) = cellvolume(Lattice(nml))
"""
    cellvolume(input::PWInput)

Return the volume of the cell based on the information given in a `PWInput`, in atomic unit.
"""
function cellvolume(input::PWInput)
    if input.system.ibrav == 0
        if isnothing(input.cell_parameters)
            throw(InformationNotEnough("`ibrav` is 0, must read cell parameters!"))
        else
            if optionof(input.cell_parameters) == "alat"
                # If no value of `celldm` is changed...
                if isnothing(input.system.celldm[1])
                    throw(InformationNotEnough("parameter `celldm[1]` needed!"))
                else
                    return input.system.celldm[1]^3 * abs(det(input.cell_parameters.data))
                end
            else  # "bohr" or "angstrom"
                return cellvolume(input.cell_parameters)
            end
        end
    else
        return cellvolume(Lattice(input.system))
    end
end


text:
## NL
function nl(lx::Float64,ly::Float64,nx::Int,ny::Int,Ξ::Float64,β::Float64,τ::Float64=0.0,
    νn::Float64=0.0;jw::Float64=0.1,ic::Array{ComplexF64,2},dt::Float64=0.01,t_end::Float64=1000.0,savefreq::Int=20,kwargs...)
    A = acoeffs(ly,ny,Ξ,τ)
    B = bcoeffs(lx,ly,nx,ny,β,τ,νn)
    Cp,Cm = ccoeffs(lx,ly,nx,ny)
    p = [nx,ny,A,B,Cp,Cm]
    tspan = (0.0,t_end)
    prob = ODEProblem(nl_eqs!,ic,tspan,p)
    @info "Solving NL equations on $(nx-1)x$(ny-1) grid"
    @info "Parameters: Ξ = $Ξ, Δθ = $jw, β = $β, τ = $τ"
    solve(prob,RK4(),dt=dt,adaptive=false,progress=true,progress_steps=10000,save_start=true,saveat=savefreq,save_everystep=savefreq==1 ? true : false)
end

## GQL
function gql(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,Ξ::Float64,β::Float64,τ::Float64=0.0,
    νn::Float64=0.0;jw::Float64=0.1,ic::Array{ComplexF64,2},dt::Float64=0.01,t_end::Float64=1000.0,savefreq::Int=20,kwargs...)
    A = acoeffs(ly,ny,Ξ,τ)
    B = bcoeffs(lx,ly,nx,ny,β,τ,νn)
    Cp,Cm = ccoeffs(lx,ly,nx,ny,Λ)
    p = [nx,ny,Λ,A,B,Cp,Cm]
    tspan = (0.0,t_end)
    prob = ODEProblem(gql_eqs!,ic,tspan,p)
    @info "Solving GQL equations on $(nx-1)x$(ny-1) grid with Λ = $Λ"
    @info "Parameters: Ξ = $Ξ, Δθ = $jw, β = $β, τ = $τ"
    solve(prob,RK4(),dt=dt,adaptive=false,progress=true,progress_steps=10000,save_start=true,save_everystep=false,dense=false,saveat=savefreq)
end

## GCE2
function gce2(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,Ξ::Float64,β::Float64,τ::Float64=0.0,
    νn::Float64=0.0;jw::Float64=0.1,icnl::Bool=true,ic::Array{ComplexF64,2},dt::Float64=0.01,t_end::Float64=1000.0,savefreq::Int=20,saveinfo::Bool=false,saveinfofreq::Int=50,poscheck::Bool=false,poscheckfreq::Int=20,kwargs...)
    A = acoeffs(ly,ny,Ξ,τ)
    B = bcoeffs(lx,ly,nx,ny,β,τ,νn)
    Cp,Cm = ccoeffs(lx,ly,nx,ny,Λ)
    @info "Solving GCE2 equations on $(nx-1)x$(ny-1) grid with Λ = $Λ"
    @info "Parameters: Ξ = $Ξ, Δθ = $jw, β = $β, τ = $τ"
    tspan = (0.0,t_end)
    u0 = icnl == true ? ic_cumulants(nx,ny,Λ,ic) : ic_cumulants(nx,ny,Λ,1e-6,ic)
    p = [nx,ny,Λ,A,B,Cp,Cm,fill!(similar(u0.x[1]),0),fill!(similar(u0.x[2]),0),fill!(similar(u0.x[2]),0)]
    prob = ODEProblem(gce2_eqs!,u0,tspan,p)
    if poscheck && Λ < nx - 1
        poschecktimes = [tt for tt=1.0:poscheckfreq:t_end]
        condition(u,t,integrator) = t ∈ poschecktimes && !ispositive(u.x[2],nx,ny,Λ)
        affect!(integrator) = positivity!(integrator.u.x[2],nx,ny,Λ)
        cbp = DiscreteCallback(condition,affect!,save_positions=(false,false))
        # return solve(prob,RK4(),callback=cbp,tstops=poschecktimes,dt=dt,adaptive=false,progress=true,progress_steps=10000,save_start=true,save_everystep=false,dense=false,saveat=savefreq)
    end
    if saveinfo && Λ < nx - 1
        saved_values = SavedValues(Float64, Tuple{Float64,Array{Float64,1}})
        save_func(u,t,integrator) = rankis(integrator.u.x[2],nx,ny,Λ)
        saveat_array = [i for i=0.0:saveinfofreq:t_end]
        cbs = SavingCallback(save_func,saved_values,save_start=true,save_everystep=false,saveat=saveat_array)
    end
    if saveinfo && !poscheck
        @info "Saving twopoint correlation rank info"
        sol = solve(prob,RK4(),dt=dt,callback=cbs,adaptive=false,progress=true,progress_steps=10000,save_start=true,save_everystep=false,saveat=savefreq)
        return sol,saved_values
    elseif saveinfo && poscheck
        @info "Running with positivity condition and saving twopoint correlation rank info"
        cb = CallbackSet(cbp,cbs)
        sol = solve(prob,RK4(),dt=dt,callback=cb,tstops=poschecktimes,adaptive=false,progress=true,progress_steps=10000,save_start=true,save_everystep=false,saveat=savefreq)
        return sol,saved_values
    else
        return solve(prob,RK4(),dt=dt,adaptive=false,progress=true,progress_steps=10000,save_start=true,save_everystep=false,saveat=savefreq)
    end
end


text:
# Concrete definition
"""
    HeapKNNGraph{V, U, D, M}

A weighted, directed graph representing an approximate k-nearest neighbors graph
using binary max heaps to store each vertex's forward edges, allowing for
efficient updates of the candidate neighbors.
"""
struct HeapKNNGraph{V<:Integer,
                    U<:Real,
                    D<:AbstractVector,
                    M<:PreMetric} <: ApproximateKNNGraph{V, U, D, M}
    heaps::Vector{BinaryMaxHeap{HeapKNNGraphEdge{V, U}}}
    data::D
    n_neighbors::V
    metric::M
end

"""
    HeapKNNGraph(data, metric, indices, distances)

Create a HeapKNNGraph from `KxN` matrices of the indices and distances, where `indices[i, v]` and
`distances[i, v]` are the index and distance to node `v`s `i`th candidate neighbor. Note that each
column of `indices` cannot have duplicate entries, but they need not be sorted by distance.
"""
function HeapKNNGraph(data::D,
                      metric::M,
                      indices::AbstractMatrix{V},
                      distances::AbstractMatrix{U}) where {D <: AbstractVector,
                                                           M <: PreMetric,
                                                           V <: Integer,
                                                           U <: Real}
    n_neighbors = size(indices, 1)
    n_points = size(indices, 2)
    length(data) == size(indices, 2) || error("`indices` must have num columns equal to `length(data)`")
    size(indices) == size(distances) || error("`indices` and `distances` must have same shape")
    n_neighbors < n_points || error("`Must have more columns than rows`")
    # ...
    HeapType = BinaryMaxHeap{HeapKNNGraphEdge{V, U}}
    knn_heaps = HeapType[HeapType() for _ in 1:n_points]
    for v in 1:n_points, i in 1:n_neighbors
        _heappush!(knn_heaps[v], HeapKNNGraphEdge(v, indices[i, v], distances[i, v]), n_neighbors)
    end
    return HeapKNNGraph{V, U, D, M}(knn_heaps, data, n_neighbors, metric)
end

"""
    HeapKNNGraph(data, n_neighbors, metric)

Create a HeapKNNGraph by randomly sampling `n_neighbors` for each point and using `metric`
to calculate weights.
"""
function HeapKNNGraph(data::D, n_neighbors::Integer, metric::M) where {D <: AbstractVector,
                                                                       M <: PreMetric}
    # assert some invariants
    n_neighbors > 0 || error("HeapKNNGraph needs positive `n_neighbors`")
    length(data) > n_neighbors || error("HeapKNNGraph needs `length(data) > n_neighbors`")

    n_points = length(data)
    U = result_type(metric, data[1], data[1])
    HeapType = BinaryMaxHeap{HeapKNNGraphEdge{typeof(n_neighbors), U}}
    knn_heaps = HeapType[HeapType() for _ in 1:n_points]
    # create approximate knn heaps for each point
    for i in eachindex(data)
        k_idxs = sample_neighbors(n_points, n_neighbors, exclude=[i])
        for j in k_idxs
            dist = evaluate(metric, data[i], data[j])
            # calculate reverese distance if dist measure is not symmetric
            rev_dist = metric isa SemiMetric ? dist : evaluate(metric, data[j], data[i])
            _heappush!(knn_heaps[i], HeapKNNGraphEdge(i, j, dist), n_neighbors)
            _heappush!(knn_heaps[j], HeapKNNGraphEdge(j, i, rev_dist), n_neighbors)
        end
    end
    return HeapKNNGraph{typeof(n_neighbors), U, D, M}(knn_heaps, data, n_neighbors, metric)
end

# convert matrix data input to vector of vectors (views)
function HeapKNNGraph(data::D, args...; kwargs...) where D <: AbstractMatrix
    return HeapKNNGraph(collect(eachcol(data)), args...; kwargs...)
end

# lightgraphs interface
"""
    edges(g::HeapKNNGraph)

Return an iterator of the edges in this KNN graph. Note that mutating the edges
while iterating through `edges(g)` may result in undefined behavior.
"""
@inline function LightGraphs.edges(g::Union{HeapKNNGraph, LockHeapKNNGraph})
    return (e for heap in g.heaps for e in heap.valtree)
end

"""
    vertices(g::HeapKNNGraph)

Return an iterator of the vertices in the KNN graph.
"""
@inline function LightGraphs.vertices(g::Union{HeapKNNGraph{V}, LockHeapKNNGraph{V}}) where V
    return one(V):V(length(g.heaps))
end

"""
    weights(g::HeapKNNGraph{V, U})

Returns a SparseMatrixCSC{U, V} `w` where `w[j, i]` returns the weight of the
directed edge from vertex `i` to `j`.
"""
function LightGraphs.weights(g::Union{HeapKNNGraph{V, U}, LockHeapKNNGraph{V, U}}) where {V, U}
    # dests, srcs, ws where weights[dests[k], srcs[k]] = ws[k]
    srcs, dsts = V[], V[]
    vals = U[]
    for e in edges(g)
        push!(srcs, src(e))
        push!(dsts, dst(e))
        push!(vals, weight(e))
    end
    return sparse(dsts, srcs, vals, nv(g), nv(g))
end

"""
    edgetype(g::HeapKNNGraph)

Return the type of the edges in `g`.
"""
function LightGraphs.edgetype(g::Union{HeapKNNGraph{V, U}, LockHeapKNNGraph{V, U}}) where {V, U}
    return HeapKNNGraphEdge{V, U}
end

"""
    has_edge(g::HeapKNNGraph, s, d)

Return true if the graph g has an edge from `s` to `d`, else false.
"""
@inline function LightGraphs.has_edge(g::Union{HeapKNNGraph{V}, LockHeapKNNGraph{V}}, s, d) where V
    for e in g.heaps[s].valtree
        if dst(e) == d
            return true
        end
    end
    return false
end
"""
    has_edge(g::HeapKNNGraph, e::HeapKNNGraphEdge)

Return true if `g` contains the edge `e`.
"""
@inline function LightGraphs.has_edge(g::Union{HeapKNNGraph{V, U}, LockHeapKNNGraph{V, U}},
                                      e::HeapKNNGraphEdge{V, U}) where {V, U}
    return e in g.heaps[src(e)].valtree
end

"""
    has_vertex(g::HeapKNNGraph, v)

Return true if vertex `v` is in `g`.
"""
@inline LightGraphs.has_vertex(g::Union{HeapKNNGraph{V}, LockHeapKNNGraph{V}}, v::V) where V = v in 1:nv(g)

"""
    ne(g::HeapKNNGraph)

Return the number of edges in `g`.
"""
@inline LightGraphs.ne(g::Union{HeapKNNGraph{V, U}, LockHeapKNNGraph{V, U}}) where {V, U} = g.n_neighbors*nv(g)

"""
    nv(g::HeapKNNGraph)

Return the number of vertices in `g`.
"""
@inline LightGraphs.nv(g::Union{HeapKNNGraph, LockHeapKNNGraph}) = length(g.heaps)

"""
    inneighbors(g::HeapKNNGraph, v)

Return a list of the neighbors connected to `v` by an incoming edge.

**Implementation Notes**

HeapKNNGraph doesn't store inneighbors directly; it must find them by iterating
over the outgoing edges for each vertex and saving those where `v == dst(e)`.
Thus, this has time complexity `𝒪(nv(g)*K)`.
"""
@inline function LightGraphs.inneighbors(g::HeapKNNGraph{V}, v::V) where V
    return collect(src(e) for e in edges(g) if dst(e) == v)
end

"""
    outneighbors(g::HeapKNNGraph, v)

Return a list of the neighbors of `v` connected by outgoing edges.

**Implementation Notes**

HeapKNNGraph stores each vertex's outgoing edges in a heap, so this has a time
complexity of `𝒪(K)`.
"""
@inline LightGraphs.outneighbors(g::HeapKNNGraph{V}, v::V) where V = dst.(g.heaps[v].valtree)

"""
    add_edge!(g::HeapKNNGraph, e::HeapKNNGraphEdge)

Try to add an edge `e`, indicating a new candidate nearest neighbor `e.dest` for
vertex `e.src`, by pushing onto `e.src`'s heap. This will fail (return `false`) if
`e` is already in the heap or if `e.weight > first(heap).weight`. Otherwise return `true`.
"""
function LightGraphs.add_edge!(g::HeapKNNGraph, e::HeapKNNGraphEdge)
    # NOTE we can assume the invariants for heap knn graphs hold
    if e < first(g.heaps[src(e)]) && !has_edge(g, src(e), dst(e))
        # we know this edge is smaller than the first, so we can start by removing that
        pop!(g.heaps[src(e)])
        push!(g.heaps[src(e)], e)
        return true
    end
    return false
end

# KNNGraphs public methods
"""
    knn_diameter(graph, v) -> diameter

Return the diameter of the set of KNNs of vertex `v`.
"""
function knn_diameter(g::Union{HeapKNNGraph{V}, LockHeapKNNGraph{V}}, v::V) where V
    return 2 * weight(first(g.heaps[v]))
end

"""
    knn_matrices(graph) -> indices, distances

Return the indices and distances of the approximate KNNs as dense
matrices where indices[j, i] and distances[j, i] are the index of
and distance to node i's jth nearest neighbor, respectively.
"""
function knn_matrices(g::Union{HeapKNNGraph{V, U}, LockHeapKNNGraph{V, U}}) where {V, U}
    indices = Matrix{V}(undef, g.n_neighbors, nv(g))
    distances = Matrix{U}(undef, g.n_neighbors, nv(g))
    for i in vertices(g)
        i_edges = sort!(node_edges(g, i))
        for (j, e) in enumerate(i_edges)
            indices[j, i] = dst(e)
            distances[j, i] = weight(e)
        end
    end
    return indices, distances
end

"""
    edge_indices(graph) -> CartesianIndices

Return the indices of the KNNs for each node `v` as tuples (v, i). To be used with
`node_edge(graph, v, i)`.
"""
function edge_indices(g::ApproximateKNNGraph)
    return CartesianIndices((nv(g), g.n_neighbors))
end

"""
    node_edge(graph, v, i) -> edge

Return the ith outgoing edge from node `v`. No ordering of the
edges is guaranteed; in particular, node_edge(graph, v, 1) is not
guaranteed to be the edge to v's nearest neighbor.
"""
@inline function node_edge(g::Union{HeapKNNGraph{V}, LockHeapKNNGraph{V}}, v::V, i::V) where V
    return g.heaps[v].valtree[i]
end

"""
    node_edges(graph, v) -> edges

Return all the outgoing edges from node v in an arbitrary order.
"""
@inline function node_edges(g::HeapKNNGraph{V}, v::V) where V
    return edgetype(g)[node_edge(g, v, i) for i in one(V):g.n_neighbors]
end

"""
    update_flag!(graph, v, i, new_flag)

Update the flag of node `v`s ith outgoing edge. Returns the new edge.
Note that since the flags don't influence the edge ordering, this can't invalidate the heap
invariant.
"""
function update_flag!(g::HeapKNNGraph{V}, v::V, i::V, new_flag::Bool) where V
    edge = node_edge(g, v, i)
    new_edge = edgetype(g)(src(edge), dst(edge), weight(edge), new_flag)
    g.heaps[v].valtree[i] = new_edge
    return new_edge
end


text:
function mean(x::Vector{Real})
    n = size(x, 1)
    if n > 0
        return sum(x) / n
    else
        return 0
    end
end

function _exp_c(x::Float64)
    ccall(:exp, Float64, (Float64,), x)
end

function _log_c(x::Float64)
    ccall(:log, Float64, (Float64,), x)
end

function _log1p_c(x::Float64)
    ccall(:log1p, Float64, (Float64,), x)
end

function _sig_c(x::Float64)
    1 / (1 + _exp_c(-x))
end

function _sig_cplus(x::Float64)
    1 / (1 + _exp_c(x))
end

function _gemmblas(A::Matrix{Float64}, B::Matrix{Float64}, ol_max::Matrix{Float64})
    sa = size(A)
    ol = copy(ol_max)
    ccall(
        ("dgemm_64_", "libopenblas64_"),
        Cvoid,
        (
            Ref{UInt8},
            Ref{UInt8},
            Ref{Int64},
            Ref{Int64},
            Ref{Int64},
            Ref{Float64},
            Ptr{Float64},
            Ref{Int64},
            Ptr{Float64},
            Ref{Int64},
            Ref{Float64},
            Ptr{Float64},
            Ref{Int64},
        ),
        'T',
        'N',
        sa[2],
        sa[2],
        sa[1],
        1.0,
        A,
        max(1, stride(A, 2)),
        B,
        max(1, stride(B, 2)),
        -one(Float64),
        ol,
        max(1, stride(ol, 2)),
    )
    return ol::Matrix{Float64}
end

function _gemvblas(
    A::Matrix{Float64},
    x::Vector{Float64},
    cons::Vector{Float64},
    sizex::Int64,
)
    ccall(
        ("dgemv_64_", "libopenblas64_"),
        Cvoid,
        (
            Ref{UInt8},
            Ref{Int64},
            Ref{Int64},
            Ref{Float64},
            Ptr{Float64},
            Ref{Int64},
            Ptr{Float64},
            Ref{Int64},
            Ref{Float64},
            Ptr{Float64},
            Ref{Int64},
        ),
        'N',
        size(A, 1),
        sizex,
        1.0,
        A,
        max(1, stride(A, 2)),
        x,
        1,
        -1.0,
        cons,
        1,
    )
    return cons
end

function _gemvblasT(A::Matrix{Float64}, x::Vector{Float64}, sizeA2::Int64)
    cons = zeros(Float64, sizeA2)
    ccall(
        ("dgemv_64_", "libopenblas64_"),
        Cvoid,
        (
            Ref{UInt8},
            Ref{Int64},
            Ref{Int64},
            Ref{Float64},
            Ptr{Float64},
            Ref{Int64},
            Ptr{Float64},
            Ref{Int64},
            Ref{Float64},
            Ptr{Float64},
            Ref{Int64},
        ),
        'T',
        size(A, 1),
        size(A, 2),
        1.0,
        A,
        max(1, stride(A, 2)),
        x,
        1,
        -1.0,
        cons,
        1,
    )
    return cons
end

function _combinations(T::Int64)
    comb = vcat([[[t, i] for i = (t+1):T] for t = 1:T]...)
end

function _myqle(x::Vector{Float64}, lx::Int64, k::Int64)
    y = ones(k) * Inf
    n = 1
    y[1] = copy(x[1])
    for i = 2:lx
        current = copy(x[i])
        k2 = 1
        while k2 <= k
            if y[k2] > current
                y[(k2+1):n] .= y[k2:(n-1)]
                y[k2] = copy(current)
                k2 = k + 1
                if n < k
                    n += 1
                end
            end
            k2 += 1
        end
    end
    return y[k]
end

function _mycopy(NH::Neighbourhood, NH_new::Neighbourhood)
    NH_new.x = copy(NH.x)
    NH_new.f = copy(NH.f)
    NH_new.obj = copy(NH.obj)
    NH_new.infeas = copy(NH.infeas)
    NH_new.ol = copy(NH.ol)
    NH_new.iu = copy(NH.iu)
    return NH_new::Neighbourhood
end

function _comp_f(NH::Neighbourhood, OptFeas::Float64)
    return (1 - OptFeas) * (sum(NH.infeas + NH.ol) + NH.iu) - OptFeas * (minimum(NH.obj))
end

function _myqle_simp(x::Vector{Float64}, ind::Vector{Float64})
    R = size(x, 1)
    alphaR = Int.(ceil.(R .* (ind)))
    alphaR[findall(alphaR .< 1)] .= 1
    alphaR[findall(alphaR .> R)] .= R
    αQle = Inf
    xsorted = sort(x)
    return sort(x)[alphaR]
end

function _cutR(
    x;
    start = "minimum",
    stop = "maximum",
    nBins = 2,
    returnBreaks = true,
    returnMidPts = false,
)
    if (start == "minimum")
        start = minimum(x)
    end
    if (stop == "maximum")
        stop = maximum(x)
    end
    bw = (stop - start) / (nBins - 1)
    midPts = zeros(nBins)
    for i = 1:nBins
        midPts[i] = start + (i - 1) * bw
    end
    breaks = collect(range(start - (bw / 2); length = nBins + 1, stop = stop + (bw / 2)))
    y = zeros(size(x, 1))
    for j = 1:size(x, 1)
        for i = 1:nBins
            if (x[j] >= breaks[i]) && (x[j] < breaks[i+1])
                y[j] = i
            end
            if i == nBins && x[j] == breaks[i+1]
                y[j] = i
            end
        end
    end
    if (returnBreaks == true || returnMidPts == true)
        if returnMidPts == false
            return (Int.(y), breaks)
        elseif returnBreaks == false
            return (Int.(y), midPts)
        else
            return (Int.(y), breaks, midPts)
        end
    else
        return Int.(y)
    end
end

function _desume_pars(
    pars::DataFrames.DataFrame;
    model = "2PL", #1PL, 2PL, 3PL
    parametrization = "at-b",
) #"at-b, at-ab, at+b, at+ab")
    n_items = size(pars, 1)
    if (model == "1PL")
        ("b" in names(pars)) ||
            ("d" in names(pars)) ||
            error("difficulty parameter b or d not defined")
        c = zeros(Float64, n_items)
        a = ones(Float64, n_items)
    elseif (model == "2PL")
        ("a" in names(pars) || "a1" in names(pars)) ||
            error("discrimination parameter a1 or a not defined")
        ("b" in names(pars)) ||
            ("d" in names(pars)) ||
            error("difficulty parameter b or d not defined")
        if "a" in names(pars)
            a = pars.a
        else
            a = pars.a1
        end
        c = zeros(Float64, n_items)
    elseif (model == "3PL")
        ("a1" in names(pars) || "a" in names(pars)) ||
            error("discrimination parameter a1 not defined")
        ("b" in names(pars)) ||
            ("d" in names(pars)) ||
            error("difficulty parameter b or d not defined")
        ("c" in names(pars)) ||
            ("g" in names(pars)) ||
            error("guessing parameter c or g not defined")
        if "a" in names(pars)
            a = pars.a
        else
            a = pars.a1
        end
        if ("c" in names(pars))
            c = pars.c
        else
            c = pars.g
        end
    elseif (model == "grm")
        ("a1" in names(pars) || "a" in names(pars)) ||
            error("discrimination parameter a1 not defined")
    end
    b = zeros(Float64, n_items, 0)
    nb = 0
    # TODO grm
    for n in names(pars)
        if startswith(string(n), "b")
            nb += 1
            b = hcat(b, pars[!, n])
        end
    end
    # nd = 0
    # d = zeros(Float64, size(pars, 1), 1)
    # for n in names(pars)
    #     if startswith(string(n), "d")
    #         nd += 1
    #         d = hcat(d, pars[!, n])
    #     end
    # end
    # if nd > 0
    #     d = mapslices(x -> b .- x, d; dims = 2)
    # end
    a2 = ones(n_items) #default: "at+b" a*theta + b
    if parametrization == "at-b" #a*theta - b
        b = -b
    elseif parametrization == "at-ab" #a*(theta - b)
        a2 = copy(a)
        b = -b
    elseif parametrization == "at+ab" #a*(theta + b)
        a2 = copy(a)
    end
    return a, a2, b, c
end

"""
	fs_to_items(xₜ, n_items, fs_items)

# Description

Ungroup items grouped by friend sets.

# Arguments

- **`xₜ::Vector{Float64}`** : Required. grouped items 0-1 vector.
- **`n_items::Int64`** : Required. Number if items.
- **`fs_items::Vector{Vector{Int64}}`** : Required. vector of items included in the friend sets.

# Output

- **`x_Iₜ::Vector{Float64}`** Returns a 0-1 vector at item level.

"""
function fs_to_items(xₜ::Vector{Float64}, n_items::Int64, fs_items::Vector{Vector{Int64}})
    xₜ_taken = findall(xₜ .== one(Float64))
    x_Iₜ = zeros(Float64, n_items)
    for i in xₜ_taken
        x_Iₜ[fs_items[i]] .= one(Float64)
    end
    return x_Iₜ::Vector{Float64}
end

"""
	item_char(pars::DataFrames.DataFrame, theta::Vector{Float64}; model = "2PL", parametrization = "at-b", D = 1, derivatives = false)

# Description

Compute the item characteristic function (probability of correct answer).

# Arguments

- **`pars::DataFrames.DataFrame`** : Required. DataFrame with item parameters (difficulty: b or d, discrimination: a or a1, guessing: c or g).
- **`theta::Vector{Float64}`** : Required. Vector of ability points.
- **`model`** : Optional. Default: `"2PL`". Values:  `"1PL`", `"2PL`", `"3PL`". IRT model.
- **`parametrization`** : Optional. Default: \"at-b\". Values:  \"at-b\", \"at-ab\", \"at+b\", \"at+ab\". IRT model parametrization. Ex: at-b is ``Da(\\theta-b)``.
- **`D`** : Optional. Default: 1.
- **`derivatives`** : Optional. Default: false. If it's true compute the derivatives. Ow a `zeros(0,0)` matrix is returned.

# Output

- **`p::Matrix{Float64}`**: Matrix of probabilites.
- **`pder::Matrix{Float64}`**: Matrix of derivatives.
"""
function item_char(
    pars::DataFrames.DataFrame,
    theta::Vector{Float64};
    model = "2PL", #1PL, 2PL, 3PL
    parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
    D = 1,#any number (1, 1.7)
    derivatives = false,
)

    n_items = size(pars, 1)
    pder = zeros(0, 0)
    a, a2, b, c = _desume_pars(pars; model = model, parametrization = parametrization)
    nb = size(b, 2)
    N = size(theta, 1)
    p =
        _sig_c.([
            D * (a[i] * theta[n] + a2[i] * b[i, j]) for i = 1:n_items, n = 1:N, j = 1:nb
        ])
    #pder =  eachslice(((1 .- p) .* p), dims = 1) .* a

    if model != "grm"
        p = c .+ ((1 .- c) .* p)
    else
        # TODO grm
        #generalized PCM, does not work
        # p[:, :, 1] = 1 .- p[:, :, 1]
        # pder[:, :, 1] = .- pder[:, :, 1]
        # for k = (nb - 1) :  -1 : 2
        # 	p[:, :, k] = p[:, :, k] - p[:, :, k-1]
        # 	pder[:, :, k] = pder[:, :, k] - pder[:, :, k-1]
        # end
    end
    if derivatives
        pder = mapslices(x -> (1 .- x) .* ((x .- c) ./ (1 .- c)) .* a, p; dims = 1)
    end
    return p, pder
end

"""
	item_char(pars::DataFrames.DataFrame, theta::Float64; model = "2PL", parametrization = "at-b", D = 1, derivatives = false)

# Description

Compute the item characteristic function (probability of correct answer).

# Arguments

- **`pars::DataFrames.DataFrame`** : Required. DataFrame with item parameters (difficulty: b or d, discrimination: a or a1, guessing: c or g).
- **`theta::Float64`** : Required. Ability point.
- **`model`** : Optional. Default: `"2PL`". Values:  `"1PL`", `"2PL`", `"3PL`". IRT model.
- **`parametrization`** : Optional. Default: \"at-b\". Values:  \"at-b\", \"at-ab\", \"at+b\", \"at+ab\". IRT model parametrization. Ex: at-b is ``Da(\\theta-b)``.
- **`D`** : Optional. Default: 1.
- **`derivatives`** : Optional. Default: false. If it's true compute the derivatives. Ow a `zeros(0,0)` matrix is returned.

# Output

- **`p::Matrix{Float64}`**: Matrix of probabilites.
- **`pder::Matrix{Float64}`**: Matrix of derivatives.
"""
function item_char(
    pars::DataFrames.DataFrame,
    theta::Float64;
    model = "2PL", #1PL, 2PL, 3PL
    parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
    D = 1,#any number (1, 1.6)
    derivatives = false,
)
    return item_char(
        pars,
        [theta],
        model = model, #1PL, 2PL, 3PL
        parametrization = parametrization, #"at-b, at-ab, at+b, at+ab"
        D = D,#any number (1, 1.6)
        derivatives = derivatives,
    )
end

"""
    item_info(
        pars::DataFrames.DataFrame,
        theta::Vector{Float64};
        model = "2PL", #1PL, 2PL, 3PL, grm
        parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
        D = 1,
    )

# Description

Compute the item Fisher information function.

# Arguments

- **`pars::DataFrames.DataFrame`** : Required. DataFrame with item parameters (difficulty: b or d, discrimination: a or a1, guessing: c or g).
- **`theta::Vector{Float64}`** : Required. Vector of ability points.
- **`model`** : Optional. Default: `"2PL`". Values:  `"1PL`", `"2PL`", `"3PL`". IRT model.
- **`parametrization`** : Optional. Default: \"at-b\". Values:  \"at-b\", \"at-ab\", \"at+b\", \"at+ab\". IRT model parametrization. Ex: at-b is ``Da(\\theta-b)``.
- **`D`** : Optional. Default: 1.

# Output

- **`i::Matrix{Float64}`**: Matrix of the item information.
"""
function item_info(
    pars::DataFrames.DataFrame,
    theta::Vector{Float64};
    model = "2PL", #1PL, 2PL, 3PL, grm
    parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
    D = 1,
) #true, false

    n_items = size(pars, 1)
    a, a2, b, c = _desume_pars(pars; model = model, parametrization = parametrization)
    p, pder = item_char(
        pars,
        theta;
        model = model,
        parametrization = parametrization,
        derivatives = true,
    )
    nb = size(p, 3)
    if model != "grm"
        #i = (a.^2 ) .* ((1 .- p) ./ p) .* ((p .- c) ./ (1 .- c)).^2
        i = pder .^ 2 ./ (p .* (1 .- p))
    else
        error("GRM still not supported.")
        #i = pder .^ 2 ./ p
    end
    i = sum(i, dims = 3)[:, :, 1]
    return i::Matrix{Float64}
end

"""
    item_info(
        pars::DataFrames.DataFrame,
        theta::Float64;
        model = "2PL", #1PL, 2PL, 3PL, grm
        parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
        D = 1,
    )

# Description

Compute the item Fisher information function.

# Arguments

- **`pars::DataFrames.DataFrame`** : Required. DataFrame with item parameters (difficulty: b or d, discrimination: a or a1, guessing: c or g).
- **`theta::Float64`** : Required. Vector of ability points.
- **`model`** : Optional. Default: `"2PL`". Values:  `"1PL`", `"2PL`", `"3PL`". IRT model.
- **`parametrization`** : Optional. Default: \"at-b\". Values:  \"at-b\", \"at-ab\", \"at+b\", \"at+ab\". IRT model parametrization. Ex: at-b is ``Da(\\theta-b)``.
- **`D`** : Optional. Default: 1.

# Output

- **`i::Matrix{Float64}`**: Matrix of the item information.
"""
function item_info(
    pars::DataFrames.DataFrame,
    theta::Float64;
    model = "2PL", #1PL, 2PL, 3PL, grm
    parametrization = "at-b", #"at-b, at-ab, at+b, at+ab"
    D = 1,
) #true, false
    return item_info(pars, [theta], model = model, parametrization = parametrization, D = D)[
        :,
        1,
        :,
    ]
end

function student_likelihood(
    f::Float64,
    r::Vector{Float64},
    i_index::Vector{Int64},
    a::Vector{Float64},
    b::Vector{Float64},
    theta::Float64;
    logL = false,
)
    n_items = size(a, 1)
    likel = 0
    for i in i_index
        p = b[i] + a[i] * theta
        likel += r[i] * p - _log1p_c(_exp_c(p))
    end
    return likel::Float64
end

include("resp_gen.jl")


text:
function collision_cross_section_pisigma11(T, nch)

    # pisigma
    pisigma11 = zeros(nch, nch)

    # [A, B, C, D]
    A = 0.0
    B = 0.0
    C = 0.0
    D = 0.0

    # A review of reaction rates and thermodynamic and transport properties
    # for an 11-species air model for chemical and
    # thermal nosequilibrium calculations to 30,000 K

    # N2-N2
    A = 0.0
    B = -0.0112
    C = -0.1182
    D = 4.8464

    pisigma11[1,1] = pisigma(T, A, B, C, D)

    # N2-N
    A = 0.0
    B = -0.0194
    C = 0.0119
    D = 4.1055

    pisigma11[1,2] = pisigma(T, A, B, C, D)
    pisigma11[2,1] = 0.0

    # N-N
    A = 0.0
    B = -0.0033
    C = -0.0572
    D = 5.0452

    pisigma11[2,2] = pisigma(T, A, B, C, D)

    return pisigma11
end

function collision_cross_section_pisigma22(T, nch)

    # pisigma
    pisigma22 = zeros(nch, nch)

    # [A, B, C, D]
    A = 0.0
    B = 0.0
    C = 0.0
    D = 0.0

    # A review of reaction rates and thermodynamic and transport properties
    # for an 11-species air model for chemical and
    # thermal nosequilibrium calculations to 30,000 K

    # N2-N2
    A = 0.0
    B = -0.0203
    C = 0.0683
    D = 4.0900

    pisigma22[1,1] = pisigma(T, A, B, C, D)

    # N2-N
    A = 0.0
    B = -0.0190
    C = 0.0239
    D = 4.1782

    pisigma22[1,2] = pisigma(T, A, B, C, D)
    pisigma22[2,1] = 0.0

    # N-N
    A = 0.0
    B = -0.0118
    C = -0.0960
    D = 4.3252

    pisigma22[2,2] = pisigma(T, A, B, C, D)

    return pisigma22
end

function pisigma(T, A, B, C, D)
    ans = exp(D) * T^(A^2*log(T) + B*log(T) + C)
    return ans
end

function freaction_rate(Tfr, nreaction)
    kf = zeros(nreaction)

    Cr = set_Cr()
    sr = set_sr()
    tr = set_tr()

    for r in 1:nreaction
        kf[r] = reaction_rate_arrhcnius(Tfr, Cr[r], sr[r], tr[r])
    end

    return kf
end

function set_Cr()
    # N2+N2 -> 2N+N2
    # N2+N  -> 2N+N

    Cr = [7.0*10^(21), 3.0*10^(22)]
    return Cr
end

function set_sr()
    # N2+N2 -> 2N+N2
    # N2+N  -> 2N+N

    sr = [-1.60, -1.60]
    return sr
end

function set_tr()
    # N2+N2 -> 2N+N2
    # N2+N  -> 2N+N

    tr = [113200, 113200]
    return tr
end

function reaction_rate_arrhcnius(Tfr, Cr, sr, tr)
    k = Cr * Tfr^(sr) * exp(-tr/Tfr)
    return k
end

function k_eq(T, nsb, nreaction)
    keq = zeros(nreaction)

    ar1, ar2, ar3, ar4, ar5 = set_ar()

    keq[1] = keq_curve_fitting(T, ar1[nsb], ar2[nsb], ar3[nsb], ar4[nsb], ar5[nsb])

    # N2+N  -> 2N+N
    # 上式と一緒
    keq[2] = keq_curve_fitting(T, ar1[nsb], ar2[nsb], ar3[nsb], ar4[nsb], ar5[nsb])

    return keq
end

function set_ar()
    # Equilibrium constants curve-fits
    # C. Park, “Nonequilibrium Hypersonic Aerothermodynamics”, Wiley, New York, 1990.

    # N2+N2 -> 2N+N2
    ar1 = [3.490700,
            2.072300,
            1.606000,
            1.535100,
            1.476600,
            1.476600]

    ar2 = [8.313300,
            1.389700,
            1.573200,
            1.606100,
            1.629100,
            1.629100]

    ar3 = [4.097800,
            2.061700,
            1.392300,
            1.299300,
            1.215300,
            1.215300]

    ar4 = [-1.272800,
            -1.182800,
            -1.153300,
            -1.149400,
            -1.145700,
            -1.145700]

    ar5 = [7.487000,
            1.510500,
            -4.543000,
            -6.980000,
            -9.444000,
            -9.444000]

    return ar1, ar2, ar3, ar4, ar5
end

function keq_curve_fitting(T, A1, A2, A3, A4, A5)
    Z = 10^4/T
    keq = exp(A1/Z + A2 + A3*log(Z) + A4*Z + A5*Z^2)
    keq = min(keq, 700.0)
    return keq
end

function breaction_rate(kf, keq, nreaction)
    kb = zeros(nreaction)
    for i in 1:nreaction
        kb[i] = kf[i] / keq[i]
    end
    return kb
end

function chemical_enthalpy(T)
    # 複数の温度モードを考慮した方がよい
    delta_hs_0 = set_delta_hs_0()
    mw  = set_mw()
    R   = set_gas_const()
    nch = length(delta_hs_0)

    theta_vib = set_thermal_vib()
    num_molecule = length(theta_vib)     # 二原子分子の数

    Rs = zeros(nch)              # 化学種気体定数
    for s in 1:nch
        Rs[s] = R / mw[s]
    end

    hs = zeros(nch)
    for s in 1:num_molecule
        etr   = 1.5 * Rs[s] * T
        erot  = Rs[s] * T
        evib  = Rs[s] * theta_vib[s] / (exp(theta_vib[s]/T)-1)
        hs[s] = etr + erot + evib + delta_hs_0[s] + Rs[s]*T
    end
    for s in num_molecule+1:nch
        etr   = 1.5 * Rs[s] * T
        hs[s] = etr + delta_hs_0[s] + Rs[s]*T
    end

    return hs
end

function set_delta_hs_0()
    # N2, N
    delta_hs_0 = [0.0, 3.364e7]   # [J/kg]
    return delta_hs_0
end

function set_mw()
    # N2, N
    mw  = [28e-3, 14e-3]          # [kg/mol]
    return mw
end

function set_gas_const()
    R   = 8.314                   # [J/K/mol]
    return R
end

function set_thermal_vib()
    # N2
    theta_vib = [3.353]
    return theta_vib
end

function set_fmol()
    # N2 + N2 --> N + N + N2
    # N2 + N  --> N + N + N
    fmol = [[2, 0],
            [1, 1]]
    return fmol
end

function set_bmol()
    # N2 + N2 <-- N + N + N2
    # N2 + N  <-- N + N + N
    bmol = [[1, 2],
            [0, 3]]
    return bmol
end

function set_avgdro_const()
    avgdro  = 6.022*10^23    # アボガドロ定数
    return avgdro
end

function set_boltzmann_const()
    kb = 1.380/10^(23)       # ボルツマン定数
    return kb
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
module aoc_25
    using Test

    export transformer, demo_card_public_key, demo_door_public_key
    export transform_subject, input
    export test_part1
    const transformer = 20201227

    const input = (6270530, 14540258)

#= The card transforms the subject number of 7 according to the card's secret loop size. The result is called the card's public key.
The door transforms the subject number of 7 according to the door's secret loop size. The result is called the door's public key.
The card and door use the wireless RFID signal to transmit the two public keys (your puzzle input) to the other device. Now, the card has the door's public key, and the door has the card's public key. Because you can eavesdrop on the signal, you have both public keys, but neither device's loop size.
The card transforms the subject number of the door's public key according to the card's loop size. The result is the encryption key.
The door transforms the subject number of the card's public key according to the door's loop size. The result is the same encryption key as the card calculated. =#

    const demo_card_public_key = 5764801
    const demo_door_public_key = 17807724

    function transform_subject(subject, max_loop_size, divisor = transformer; value::Int = 1)
        #out = Vector{Int}(undef, max_loop_size )
        for loop_size = 1:max_loop_size
            value *= subject
            value = rem( value, divisor )
            #out[ loop_size ] = value
        end
        value
    end

    export get_loop_numbers
    function get_loop_numbers(public_key::Int; subject = 7)
        loop_size = 0
        value = 1
        n = 1
        while true
            value *= subject
            value = rem( value, transformer )
            if value == public_key
                loop_size = n
                break
            end
            n += 1
        end
        n
    end
    function get_loop_numbers(public_key_1, public_key_2; subject = 7)
        if public_key_1 == public_key_2
            loop_size = get_loop_numbers(public_key_1)
            return loop_size, loop_size
        end
        loop_size_1 = 0
        loop_size_2 = 0
        value = Int(1)
        n = 1
        while loop_size_1 == 0 || loop_size_2 == 0
            value *= subject
            value = rem( value, transformer )
            if value == public_key_1
                loop_size_1 = n
            elseif value == public_key_2
                loop_size_2 = n
            end
            n += 1
        end
        loop_size_1, loop_size_2
    end

    export part1
    function part1()
        loop_sizes = get_loop_numbers.(input)
        encryption_key = transform_subject.(reverse(input), loop_sizes)
        encryption_key[1]
    end

    function test_part1()
        @test 8 == get_loop_numbers( demo_card_public_key )
        @test 11 == get_loop_numbers( demo_door_public_key )
        @test 16311885 == part1()
    end

end