text:
export LaurentDirichlet

##
# LaurentDirichlet represents functions in the basis
#       1, z^(-1)+1,z+1,z^(-2)+z^(-1),z^2+z,…
# which incorporates decay at z = -1
# TODO: restructure as SumSpace made out of HardyDirichlet
##

struct LaurentDirichlet{DD,RR} <: Space{DD,RR}
    domain::DD
    LaurentDirichlet{DD,RR}(d::DD) where {DD,RR} = new(d)
end
LaurentDirichlet(d::Domain) = LaurentDirichlet{typeof(d),complex(eltype(d))}(d)
LaurentDirichlet() = LaurentDirichlet(PeriodicSegment())

spacescompatible(a::LaurentDirichlet,b::LaurentDirichlet) = domainscompatible(a,b)

canonicalspace(S::LaurentDirichlet) = Laurent(domain(S))

Conversion(a::LaurentDirichlet{DD,RR},b::Laurent{DD,RR}) where {DD,RR} = ConcreteConversion(a,b)
bandwidths(::ConcreteConversion{LaurentDirichlet{DD,RR},Laurent{DD,RR}}) where {DD,RR} = 0,2

function getindex(C::ConcreteConversion{LaurentDirichlet{DD,RR},Laurent{DD,RR}},k::Integer,j::Integer) where {DD,RR}
    if k==1 && j==2
        one(eltype(C))
    elseif k==j || j==k+2
        one(eltype(C))
    else
        zero(eltype(C))
    end
end

conversion_rule(b::LaurentDirichlet,a::Laurent{DD,RR}) where {DD,RR}=b

differentiate(f::Fun{LaurentDirichlet{DD,RR}}) where {DD,RR} = differentiate(Fun(f,Laurent))

for op in (:+,:-,:*)
    @eval begin
        $op(f::Fun{Laurent{DD,RR}},g::Fun{LaurentDirichlet{DD,RR}}) where {DD,RR} = $op(f,Fun(g,Laurent))
        $op(f::Fun{LaurentDirichlet{DD,RR}},g::Fun{Laurent{DD,RR}}) where {DD,RR} = $op(Fun(f,Laurent),g)
        $op(f::Fun{LaurentDirichlet{DD,RR}},g::Fun{LaurentDirichlet{DD,RR}}) where {DD,RR} = $op(Fun(f,Laurent),Fun(g,Laurent))
    end
end

Base.real(f::Fun{LaurentDirichlet{DD,RR}}) where {DD,RR} = real(Fun(f,Laurent))
Base.imag(f::Fun{LaurentDirichlet{DD,RR}}) where {DD,RR} = imag(Fun(f,Laurent))

coefficients(v::AbstractVector,::Laurent,::LaurentDirichlet) = laurentdirichlettransform!(copy(v))
coefficients(v::AbstractVector,::LaurentDirichlet,::Laurent) = laurentidirichlettransform!(copy(v))

function laurentdirichlettransform!(w::AbstractVector)
    if length(w) > 1
        for k=length(w)-2:-1:1
            @inbounds w[k] -= w[k+2]
        end
        @inbounds w[1] -= w[2]
    end

    w
end

function laurentidirichlettransform!(w::AbstractVector)
    if length(w) == 1
    elseif length(w) == 2
        @inbounds w[1] += w[2]
    else
        @inbounds w[1] += w[2]+w[3]
        for k=4:length(w)
            @inbounds w[k-2] += w[k]
        end
    end

    w
end

##
# CosDirichlet represents
# 1,cos(θ)+1,cos(2θ)+cos(θ),…
#

struct CosDirichlet{DD,RR} <: Space{DD,RR}
    domain::DD
end

CosDirichlet(d::Domain) = CosDirichlet{typeof(d),real(eltype(d))}(d)
CosDirichlet() = CosDirichlet(PeriodicSegment())

spacescompatible(a::CosDirichlet,b::CosDirichlet) = domainscompatible(a,b)

canonicalspace(S::CosDirichlet)=CosSpace(domain(S))

Conversion(a::CosSpace,b::CosDirichlet) = ConcreteConversion(a,b)
bandwidths(::ConcreteConversion{CD,CS}) where {CS<:CosSpace,CD<:CosDirichlet}=0,1
getindex(C::ConcreteConversion{CD,CS},k::Integer,j::Integer) where {CS<:CosSpace,CD<:CosDirichlet} =
    (k==j||j==k+1) ? one(eltype(C)) : zero(eltype(C))

conversion_rule(b::CosDirichlet,a::CosSpace)=b


text:
using Pkg; Pkg.activate(".")
using DataFrames
using Printf
using DelimitedFiles
using JSON
using GraphFusedElasticNet
using NPZ

##
splits = readtable("../gfen-reproduce/processed_data/splits_opt_pt.csv")

##
split_to_fit = 3
split = splits[split_to_fit, :]

##
fname = @sprintf("../gfen-reproduce/productivity_splits/%02d.csv", split_to_fit + 1)
splitdata = readdlm(fname, ',')

## find best lambda
fname = @sprintf("../gfen-reproduce/modelfit_metrics/cvloss_%02d.csv", split_to_fit + 1)
best_lams = readtable(fname)
best_row = argmax(best_lams.final_pred)
λsl1 = best_lams.λsl1[best_row]
λsl2 = best_lams.λsl2[best_row]
λtl1 = best_lams.λtl1[best_row]
λtl2 = best_lams.λtl2[best_row]

##
edges_df = readtable("../gfen-reproduce/processed_data/spatiotemporal_graph.csv")

##
num_nodes = size(splitdata, 1)
num_edges = size(edges_df, 1)

## nbr map
edges = [(r.vertex1 + 1, r.vertex2 + 1) for r in eachrow(edges_df)]
tv1 = [(r.temporal == 1) ? λtl1 : λsl1 for r in eachrow(edges_df)]
tv2 = [(r.temporal == 1) ? λtl2 : λsl2 for r in eachrow(edges_df)]

##
mod = BayesianBinomialGFEN(edges, tv1=tv1, tv2=tv2)

## load data and map estimate for fast mixing
fname = @sprintf("../gfen-reproduce/best_betas/betas_%02d.csv", split_to_fit)
init = vec(readdlm(fname, ','))
s = splitdata[:, 1]
a = splitdata[:, 2]

## fit model
n = 1_000
thinning = 5
burnin = 0.5

##
chain = sample_chain(mod, s, a, n, init=init, verbose=true, async=true)

##
nstart = ceil(Int, size(chain, 2) * burnin)
chain = chain[:, nstart:end]

##
writedlm()

ext: jl


text:
using LeastSquaresEstimation
using Plots

function main()
    # Define basis function vector
    h = [t->sin(t), t->cos(t), t->t, t->cos(t)*sin(t), t->t^2]

    # Instantiate Linear Least Squares Estimator
    llse = LinearLeastSquares(h)

    # Generate Measurements
    m    = 101
    y(t) = 0.3*sin(t) + 0.5*cos(t) + 0.1*t + sqrt(0.001)*randn()

    ts = zeros(m)
    ys = zeros(m)
    @inbounds for i in 1:m
        if i != 1
            ts[i] = ts[i-1] + 0.1
        end
        ys[i] = y(ts[i])
    end

    # Feed estimator batch of measurements
    FeedMeasurementBatch!(llse, ys, ts)

    # Compute estimate
    ComputeEstimate!(llse)

    # Compute estimated measurements and plot
    yhat = zeros(m)
    @inbounds for i in 1:m
        yhat[i] = GetMeasurementEstimate(llse, ts[i])
    end

    # Display
    println("Estimated parameters:")
    display(GetEstimate(llse))

    # Plot generated and estimated measurements
    plot(ts, ys; label = "Simulated measurements")
    plot!(ts, yhat; label = "Estimated measurements")
end

main()


text:
# # Motzkin

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/motzkin.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/motzkin.ipynb)
# **Adapted from**: (3.6) and (3.19) of [BPT12]
#
# [BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
# *Semidefinite Optimization and Convex Algebraic Geometry*.
# Society for Industrial and Applied Mathematics, **2012**.

# The first explicit example of nonnegative polynomial that is not a sum of squares was found by Motzkin in 1967. By the [Arithmetic-geometric mean](https://en.wikipedia.org/wiki/Arithmetic%E2%80%93geometric_mean),
# $$ \frac{x^4y^2 + x^2y^4 + 1}{3} \ge \sqrt[3]{x^4y^2 \cdot x^2y^4 \cdot 1} = x^2y^2 $$
# hence
# $$ x^4y^2 + x^2y^4 + 1 - 3x^2y^2 \ge 0. $$
# The code belows construct the Motzkin polynomial using [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl).

using Test #src
using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2

# The Motzkin polynomial is nonnegative but is not a sum of squares as we can verify numerically as follows.
# We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
@constraint(model, motzkin >= 0) # We constraint `motzkin` to be a sum of squares

optimize!(model)

# We see that the problem is detected as infeasible...

@test termination_status(model) == MOI.INFEASIBLE #src
termination_status(model)

# ... and that the dual solution is a certificate of the infeasibility of the problem.

@test dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE #src
dual_status(model)

# Even if the Motzkin polynomial is not a sum of squares, it can still be certified to be nonnegative using sums of squares.
# Indeed a polynomial is certified to be nonnegative if it is equal to a fraction of sums of squares.
# The Motzkin polynomial is equal to a fraction of sums of squares whose denominator is $x^2 + y^2$.
# This can be verified numerically as follows:

model = SOSModel(solver)
@constraint(model, (x^2 + y^2) * motzkin >= 0) # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares

optimize!(model)

# Now the problem is declared feasible by the solver...

@test termination_status(model) == MOI.OPTIMAL #src
termination_status(model)

# ... and the primal solution is a feasible point, hence it is a certificate of nonnegativity of the Motzkin polynomial.

@test primal_status(model) == MOI.FEASIBLE_POINT #src
primal_status(model)

# One may consider ourself lucky to have had the intuition that $x^2 + y^2$ would work as denominator.
# In fact, the search for the denominator can be carried out in parallel to the search of the numerator.
# In the example below, we search for a denominator with monomials of degrees from 0 to 2.
# If none is found, we can increase the maximum degree 2 to 4, 6, 8, ...
# This gives a hierarchy of programs to try in order to certify the nonnegativity of a polynomial by identifying it with a fraction of sum of squares polynomials.
# In the case of the Motzkin polynomial we now that degree 2 is enough since $x^2 + y^2$ works.

model = SOSModel(solver)
X = monomials([x, y], 0:2)

# We create a quadratic polynomial that is not necessarily a sum of squares since this is implied by the next constraint: `deno >= 1`.

@variable(model, deno, Poly(X))

# We want the denominator polynomial to be strictly positive, this prevents the trivial solution deno = 0 for instance.

@constraint(model, deno >= 1)
@constraint(model, deno * motzkin >= 0)
optimize!(model)

@test termination_status(model) == MOI.OPTIMAL #src
termination_status(model)

@test primal_status(model) == MOI.FEASIBLE_POINT #src
primal_status(model)

# We can check the denominator found by the program using `JuMP.value`

value(deno)

# Because a picture is worth a thousand words let's plot the beast.
# We can easily extend `Plots` by adding a recipe to plot bivariate polynomials.

using RecipesBase
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
    x, y, (x, y) -> p(variables(p) => [x, y])
end
import Plots
Plots.plot(
    range(-2, stop=2, length=100),
    range(-2, stop=2, length=100),
    motzkin,
    st = [:surface],
    seriescolor=:heat,
    colorbar=:none,
    clims = (-10, 80)
)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Joseki, HTTP

#=
    This is a simple Battlesnake written in Julia
    For instructions see https://github.com/Nettogrof/starter-snake-julia/blob/main/README.md

=#

#=
    Receive a Get request at  /
    Return info about Battlesnake Personalize
=#
function index(req::HTTP.Request)
    apiversion = "1"
    author = ""  # TODO: Your Battlesnake Username
    color = "#888888"  # TODO: Personalize
    headType = "default"  # TODO: Personalize
    tailType = "default"  # TODO: Personalize

    dict = Dict("color" => color, "head" => headType, "tail" => tailType, "apiversion" => apiversion, "author" => author)
    simple_json_responder(req,dict)
end

 #=
    At start of the game, your snake will receive a /start request, and this function will execute.
    Response will be ignored
 =#
function start(req::HTTP.Request)
    println("START")
    render(Text, "OK")
end

#=
    Your snake will receive at /move request at every turn. Use the information provided to determine
    how your Battlesnake will move on that turn, either up, down, left, or right.
=#
function move(req::HTTP.Request)
    #=
        Here some head-up to learn how to parse the request.  (See https://docs.battlesnake.com/references/api#move  for detail about the request)

        Your code could look like:

        json = body_as_dict(req)
        boardHeight = json["board"]["height"]
        boardWidth = json["board"]["width"]
        yourSnakeName = json["you"]["name"]

    =#

    directions = ["up", "down", "left", "right"]    #all directions possibles
    direction = rand(directions)                    #Choose randomly one direction

    shout = ""                                      #your Snake's shout

    response = Dict("move"=>direction , "shout"=>shout)
    simple_json_responder(req,response)             #Sending the response
end

 #=
    At end of the game, your snake will receive a /end request, and this function will execute.
    Response will be ignored
 =#
function endGame(req::HTTP.Request)
    println("END")
    render(Text, "OK")
end

 endpoints = [
     (index, "GET", "/"),
     (start, "POST", "/start"),
     (move, "POST", "/move"),
     (endGame, "POST", "/end")

 ]
 r = Joseki.router(endpoints)

 # Fire up the server
 println("Server up")
 HTTP.serve(r, "0.0.0.0", 8080; verbose=false, reuseaddr=true)

max_stars_repo_name: proteanblank/starter-snake-julia

max_issues_repo_name: proteanblank/starter-snake-julia

max_forks_repo_name: proteanblank/starter-snake-julia


text:
using Test

@testset "all" begin
    include("Prerequisites/Prerequisites.jl")
    include("Mathematics/Mathematics.jl")
    include("Essentials/Essentials.jl")
end


text:
module TestConnectorComp

using Mimi
using Test

import Mimi: compdef, compdefs, components, getspan

@defcomp Long begin
    x = Parameter(index=[time])
end

late_start = 2005

@defcomp Short begin
    a = Parameter()
    b = Variable(index=[time])

    function run_timestep(p, v, d, t)
        if t >= TimestepValue(late_start)
            v.b[t] = p.a * t.t
        end
    end
end

years = 2000:2010
year_dim = Mimi.Dimension(years)

#------------------------------------------------------------------------------
#  1. Use the connect_param! method with backup data (ConnectorComp gets added
#       under the hood during build)
#------------------------------------------------------------------------------

model1 = Model()
set_dimension!(model1, :time, years)
add_comp!(model1, Short, first=late_start)
add_comp!(model1, Long)
update_param!(model1, :Short, :a, 2.)
connect_param!(model1, :Long, :x, :Short, :b, zeros(length(years)))

run(model1)

@test length(components(model1.mi)) == 3    # ConnectorComp is added to the list of components in the model isntance
@test length(compdefs(model1.md)) == 2      # The ConnectorComp shows up in the model instance but not the model definition

b = model1[:Short, :b]
x = model1[:Long, :x]

# Test that all allocated datum arrays are the full length of the time dimension
@test length(b) == length(years)
@test length(x) == length(years)

@test all(ismissing, b[1:year_dim[late_start]-1])
@test all(iszero, x[1:year_dim[late_start]-1])

# Test the values are right after the late start
offset = late_start - years[1]
@test b[year_dim[late_start]:end] == x[year_dim[late_start]:end]
@test b[year_dim[late_start]:end] == (collect(year_dim[late_start]:year_dim[years[end]]) .- 5) * 2.0

@test Mimi.datum_size(model1.md, Mimi.compdef(model1.md, :Long), :x) == (length(years),)

# Test the dataframe size
b = getdataframe(model1, :Short, :b)
@test size(b) == (length(years), 2)

# Try connecting using backup_offset
connect_param!(model1, :Long, :x, :Short, :b, zeros(length(years)), backup_offset = 1)
run(model1)

b = model1[:Short, :b]
x = model1[:Long, :x]

@test all(ismissing, b[1:year_dim[late_start]-1])
@test all(iszero, x[1:year_dim[late_start]]) # one more zero because used backup an extra year
@test b[year_dim[late_start]+1:end] == x[year_dim[late_start]+1:end] # the rest of the data shoudl match

#------------------------------------------------------------------------------
#  2. Test with a short component that ends early (and test Variable timesteps)
#------------------------------------------------------------------------------

years_variable = [2000:2004..., 2005:5:2030...]
dim_variable = Mimi.Dimension(years_variable)

early_last = 2010

model2 = Model()
set_dimension!(model2, :time, years_variable)
add_comp!(model2, Short; last=early_last)
add_comp!(model2, Long)
update_param!(model2, :Short, :a, 2.)
connect_param!(model2, :Long, :x, :Short, :b, zeros(length(years_variable)))

run(model2)

@test length(components(model2.mi)) == 3
@test length(compdefs(model2.md)) == 2      # The ConnectorComp shows up in the model instance but not the model definition

b = model2[:Short, :b]
x = model2[:Long, :x]

# Test that all allocated datum arrays are the full length of the time dimension
@test length(b) == length(years_variable)
@test length(x) == length(years_variable)

@test all(ismissing, b[dim_variable[early_last]+1 : end])
@test all(iszero, x[dim_variable[early_last]+1 : end])

# Test the values are right after the late start
@test b[dim_variable[late_start] : dim_variable[early_last]] ==
    x[dim_variable[late_start] : dim_variable[early_last]] ==
    [2 * i for i in dim_variable[late_start]:dim_variable[early_last]]

#------------------------------------------------------------------------------
#  3. A model that requires multiregional ConnectorComps
#------------------------------------------------------------------------------

@defcomp Long_multi begin
    regions = Index()

    x = Parameter(index = [time, regions])
end

@defcomp Short_multi begin
    regions = Index()

    a = Parameter(index=[regions])
    b = Variable(index=[time, regions])

    function run_timestep(p, v, d, t)
        for r in d.regions
            v.b[t, r] = t.t + p.a[r]
        end
    end
end

regions = [:A, :B]

model3 = Model()
set_dimension!(model3, :time, years)
set_dimension!(model3, :regions, regions)
add_comp!(model3, Short_multi, first=late_start)
add_comp!(model3, Long_multi)
update_param!(model3, :Short_multi, :a, [1,2])
connect_param!(model3, :Long_multi, :x, :Short_multi, :b, zeros(length(years), length(regions)))

run(model3)

@test length(components(model3.mi)) == 3
@test length(compdefs(model3.md)) == 2      # The ConnectorComp shows up in the model instance but not the model definition

b = model3[:Short_multi, :b]
x = model3[:Long_multi, :x]

# Test that all allocated datum arrays are the full length of the time dimension
@test size(b) == (length(years), length(regions))
@test size(x) == (length(years), length(regions))

@test all(ismissing, b[1:year_dim[late_start]-1, :])
@test all(iszero, x[1:year_dim[late_start]-1, :])

# Test the values are right after the late start
late_yr_idxs = year_dim[late_start]:year_dim[end]

@test b[late_yr_idxs, :] == x[year_dim[late_start]:end, :]

offset = late_start - years[1]
@test b[late_yr_idxs, :] == [[i + 1 - offset for i in late_yr_idxs] [i + 2 - offset for i in late_yr_idxs]]

#------------------------------------------------------------------------------
#  4. Test where the short component starts late and ends early
#------------------------------------------------------------------------------

first, last = 2002, 2007

model4 = Model()
set_dimension!(model4, :time, years)
set_dimension!(model4, :regions, regions)
add_comp!(model4, Short_multi; first=first, last=last)
add_comp!(model4, Long_multi)

update_param!(model4, :Short_multi, :a, [1,2])
connect_param!(model4, :Long_multi=>:x, :Short_multi=>:b, zeros(length(years), length(regions)))

run(model4)

@test length(components(model4.mi)) == 3
@test length(compdefs(model4.md)) == 2      # The ConnectorComp shows up in the model instance but not the model definition

b = model4[:Short_multi, :b]
x = model4[:Long_multi, :x]

# Test that all allocated datum arrays are the full length of the time dimension
@test size(b) == (length(years), length(regions))
@test size(x) == (length(years), length(regions))

@test all(ismissing, b[1:year_dim[first]-1, :])
@test all(ismissing, b[year_dim[last]+1:end, :])
@test all(iszero, x[1:year_dim[first]-1, :])
@test all(iszero, x[year_dim[last]+1:end, :])

# Test the values are right after the late start

offset = first - years[1]
yr_idxs = year_dim[first]:year_dim[last]
@test b[yr_idxs, :] == x[yr_idxs, :]
@test b[yr_idxs, :] == [[i + 1 - offset for i in yr_idxs] [i + 2 - offset for i in yr_idxs]]

#------------------------------------------------------------------------------
#  5. Test errors with backup data
#------------------------------------------------------------------------------

late_start_long = 2002

model5 = Model()
set_dimension!(model5, :time, years)
add_comp!(model5, Short; first = late_start)
add_comp!(model5, Long; first = late_start_long)
update_param!(model5, :Short, :a, 2)

# A. test wrong size (needs to be length of model, not length of component)
@test_throws ErrorException connect_param!(model5, :Long=>:x, :Short=>:b, zeros(getspan(model5, :Long)))
@test_throws ErrorException connect_param!(model4, :Long_multi=>:x, :Short_multi=>:b, zeros(getspan(model4, :Long_multi)[1], length(regions)+1)) # test case with >1 dimension

# B. test no backup data provided
@test_throws ErrorException connect_param!(model5, :Long=>:x, :Short=>:b)   # Error because no backup data provided

#------------------------------------------------------------------------------
#  6. Test connecting Short component to Long component (does not add a
#       connector component)
#------------------------------------------------------------------------------

@defcomp foo begin
    par = Parameter(index=[time])
    var = Variable(index=[time])
    function run_timestep(p, v, d, t)
        v.var[t] = p.par[t]
    end
end

model6 = Model()
set_dimension!(model6, :time, years)
add_comp!(model6, foo, :Long; rename=[:var => :long_foo])
add_comp!(model6, foo, :Short; rename=[:var => :short_foo],first=late_start)
@test_throws ErrorException connect_param!(model6, :Short => :par, :Long => :var, backup_offset = 1) # can't use backup_offset without backup
connect_param!(model6, :Short => :par, :Long => :var)
update_param!(model6, :Long, :par, years)

run(model6)

@test length(components(model6.mi)) == 2

short_par = model6[:Short, :par]
short_var = model6[:Short, :var]

@test short_par == years    # The parameter has values instead of `missing` for years when this component doesn't run,
                            # because they are coming from the longer component that did run

@test all(ismissing, short_var[1:year_dim[late_start]-1])
@test short_var[year_dim[late_start]:end] == years[year_dim[late_start]:end]

end #module


text:
include("AbstractTypes/AbstractType.jl")

mutable struct CapacityExpansionNode <: ScenarioNode
    year::Int64
    incrementalQuantity::Float64
end

mutable struct CapacityDivisionNode <: ScenarioNode
    year::Int64
    probability::Float64
    capacityMix::Array{PhysicalNode,1}
end

mutable struct CapacityExpansionEdge <: ScenarioEdge
    #from::ScenarioNode
    #to::ScenarioNode
    probability::Float64
    gapCoverPercentage::Float64
end

mutable struct CapacityDivisionEdge <: ScenarioEdge
    #from::ScenarioNode
    #to::ScenarioNode
    probability::Float64
    typeShare::Dict{Symbol, Float64}
end

global gapCoverProbability = [(0.7,0.3),(0.7,0.3),(0.7,0.3),(0.7,0.3),(0.7,0.3)]
#global gapCoverPercentage = [(0.9,0.9),(0.9,0.9),(0.9,0.9),(0.9,0.9),(0.9,0.9)]

#capacityDivision expectation as seen by each GenCo
global divisionProbability = [(0.3,0.3,0.4),(0.3,0.3,0.4),(0.3,0.3,0.4),(0.3,0.3,0.4),(0.3,0.3,0.4)]

function NewCase(newExpansion::Float64,newDivision::Tuple{Float64,Float64,Float64})
    gapCoverPercentage = Array{Tuple,1}(undef,0)
    divisionTypeShare = Array{Tuple,1}(undef,0)
    for genCoCount in 1:5
        thisGapCover = (newExpansion,newExpansion)
        thisDivisionType = (Dict([(:Base, newDivision[1]),(:Mid, newDivision[2]),(:Peak, newDivision[3])]),
                            Dict([(:Base, newDivision[1]),(:Mid, newDivision[2]),(:Peak, newDivision[3])]),
                            Dict([(:Base, newDivision[1]),(:Mid, newDivision[2]),(:Peak, newDivision[3])]))
        push!(gapCoverPercentage,thisGapCover)
        push!(divisionTypeShare,thisDivisionType)
    end
    return gapCoverPercentage,divisionTypeShare
end


text:
function test_renaming_queries(section, queries, tab)
    @testset "$section" begin
        for (sel, res) in queries
            @test s.select(tab, sel...) == _select_rename(tab, res)
        end
    end
end

function test_queries(section, queries, tab)
    @testset "$section" begin
        for (sel, res) in queries
            @test s.select(tab, sel...) == _select_subset(tab, res)
        end
    end
end

function test_queries_and_errors(section, queries, errors, tab)
    @testset "$section" begin
        @testset "queries" begin
            for (sel, res) in queries
                @test s.select(tab, sel...) == _select_subset(tab, res)
            end
        end
        @testset "errors" begin
            for (sel, res) in errors
                @test_throws res s.select(tab, eval(sel)...)
            end
        end
    end
end

# https://github.com/JuliaLang/julia/issues/25612
function test_renaming_queries_warn(section, queries, tab)
    @testset "$section" begin
        for (query, msg) in queries
            (sel, res) = query
            @test_logs (:warn, msg) (@test s.select(tab, sel...) == _select_rename(tab, res))
        end
    end
end

function test_transforms(section, tab)
    base = Tables.columntable(tab).a
    @testset "$section" begin
        @test all((_pull(select(tab, :a => bycol.(x->x + 1)), :a) .- base) .== 1)
        @test all((_pull(select(tab, :a => byrow.((r,c)->r[c] + 1)), :a) .- base) .== 1)
        @test all((_pull(select(tab, :a => bycol(x->x .+ 1)), :a) .- base) .== 1)
        @test all((_pull(select(tab, :a => byrow((r,c)->getindex.(r, c) .+ 1)), :a) .- base) .== 1)
        @test all((_pull(select(tab, :a => bytab((t,c)->t[c] .+ 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), :a => bycol!.(x->x + 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), :a => byrow!.((r,_)->r.a + 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), :a => bycol!(x->x .+ 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), :a => byrow!((r,c)->getindex.(r, c) .+ 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), :a => bytab!((t,_)->t.a .+ 1)), :a) .- base) .== 1)
    end
end

function test_transforms(section, first_trans, tab)
    base = Tables.columntable(tab).a
    @testset "$section" begin
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => bycol.(x->x + 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => byrow.((r,c)->r[c] + 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => bycol(x->x .+ 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => byrow((r,c)->getindex.(r, c) .+ 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => bytab((t,c)->t[c] .+ 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => bycol!.(x->x + 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => byrow!.((r,_)->r.a + 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => bycol!(x->x .+ 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => byrow!((r,c)->getindex.(r, c) .+ 1)), :a) .- base) .== 2)
        @test all((_pull(select(deepcopy(tab), (:a, :b) => first_trans, :a => bytab!((t,_)->t.a .+ 1)), :a) .- base) .== 2)
    end
end

function test_newcols(section, tab)
    base = Tables.columntable(tab).a
    @testset "$section" begin
        @test all((_pull(select(tab, t = :a => bycol.(x->x + 1)), :t) .- base) .== 1)
        @test all((_pull(select(tab, t = :a => byrow.((r,c)->r.a + 1)), :t) .- base) .== 1)
        @test all((_pull(select(tab, t = :a => bycol(x->x .+ 1)), :t) .- base) .== 1)
        @test all((_pull(select(tab, t = :a => byrow((r,c)->getindex.(r, :a) .+ 1)), :t) .- base) .== 1)
        @test all((_pull(select(tab, t = :a => bytab((t,c)->t.a .+ 1)), :t) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), a = :a => bycol!.(x->x + 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), a = :a => byrow!.((r,_)->r.a + 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), a = :a => bycol!(x->x .+ 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), a = :a => byrow!((r,c)->getindex.(r, :a) .+ 1)), :a) .- base) .== 1)
        @test all((_pull(select(deepcopy(tab), a = :a => bytab!((t,_)->t.a .+ 1)), :a) .- base) .== 1)
        @test_throws KeyError select(deepcopy(tab), t = :a => bycol!.(x->x + 1))
        @test_throws KeyError select(deepcopy(tab), t = :a => byrow!.((r,_)->r.a + 1))
        @test_throws KeyError select(deepcopy(tab), t = :a => bycol!(x->x .+ 1))
        @test_throws KeyError select(deepcopy(tab), t = :a => byrow!((r,c)->getindex.(r, :a) .+ 1))
        @test_throws KeyError select(deepcopy(tab), t = :a => bytab!((t,_)->t.a .+ 1))
    end
end

symbol_queries = [
    (:a,) => [:a],
    (-cols(:b),) => [:a, :c1],
    (not(:c1),) => [:a, :b],
    (not(:a, :b, :c1),) => [],
    (:b, :a) => [:b, :a],
    ((:b, :a)) => [:b, :a],
    (-cols(:b, :c1),) => [:a],
    (not(:b, :c1),) => [:a],
    (-(cols(:b) | cols(:c1)),) => [:a],
    (not(cols(:b) & cols(:c1)),) => [:a, :b, :c1],
    (not(not(:a)),) => [:a]
]

symbol_errors = [
    :(:A,) => KeyError,
    :(cols(:A),) => KeyError,
    :((:a, :A)) => KeyError,
    :(-cols(:A),) => KeyError
]

int_queries = [
    (1,) => [:a],
    (-cols(:b),) => [:a, :c1],
    (not(3),) => [:a, :b],
    ((-3,)) => [:a, :b],
    (2, 1) => [:b, :a],
    ((2, 1)) => [:b, :a],
    (-cols(2, 3),) => [:a],
    (not(2, 3),) => [:a],
    (-(cols(2) | cols(3)),) => [:a],
    (not(cols(2) & cols(3)),) => [:a, :b, :c1]
]

int_errors = [
    :((4,)) => BoundsError,
    :((0,)) => ArgumentError,
    :((4, 1)) => BoundsError,
    :((0, -2)) => ArgumentError,
    :(-cols(4),) => BoundsError
]

mixed_queries = [
    (1, :b) => [:a , :b],
    ((2, :a)) => [:b , :a],
    ([2, :a],) => [:b, :a],
    (-cols(2, :c1),) => [:a],
    (!cols(2, :c1),) => [:a],
    (~cols(2, :c1),) => [:a],
    (cols(:b, 1),) => [:b, :a],
    ([2, 1],) => [:b, :a]
]

mixed_errors = [
    :(:c1, 4) => BoundsError,
    :(0, :a) => ArgumentError,
    :(1, :D) => KeyError
]

int_range_queries = [
    (1:1,) => [:a],
    (!cols(2:2),) => [:a, :c1],
    (-3:-3,) => [:a, :b],
    (1:2,) => [:a, :b],
    (1:2:3,) => [:a, :c1],
    (-3:2:-1,) => [:b]
]

int_range_errors = [
    :(0:3,) => ArgumentError,
    :(4:4,) => BoundsError,
    :(-4:-4) => BoundsError
]

symbol_range_queries = [
    (colrange(:a, :a),) => [:a], (!colrange(:b, :b),) => [:a, :c1],
    (colrange(:c1, :b),) => [:c1, :b], (colrange(:a, :b),) => [:a, :b],
    (colrange(:c1, :a, by=2),) => [:c1, :a],
    (colrange(1, 1),) => [:a], (!colrange(2, 2),) => [:a, :c1],
    (colrange(3, 2),) => [:c1, :b], (colrange(1, 2),) => [:a, :b],
    (colrange(3, 1, by=2),) => [:c1, :a]
]

symbol_range_errors = [
    :(colrange(:a, :d),) => KeyError,
    :(colrange(:d, :d),) => KeyError,
    :(-colrange(:d, :d),) => KeyError,
    :(colrange(1, 4),) => BoundsError,
    :(colrange(0, 0),) => ArgumentError,
    :(-colrange(4, 0),) => ArgumentError
]

bool_queries = [
    ([true, true, true],) => [:a, :b, :c1], (:b, [false, false, false], :a) => [:b, :a],
    ([false, true, false],) => [:b],  (!cols([true, false, true]),) => [:b],
    (.![false, true, false],) => [:a, :c1]
]

bool_errors = [
    :([true],) => BoundsError,
    :([false, false, false, false],) => BoundsError
]

regex_queries = [
    (if_matches(r"."),) => [:a, :b, :c1], (if_matches(r"a"),) => [:a],
    (-if_matches(r"b"),) => [:a, :c1], (!if_matches(r"c"),) => [:a, :b],
    (if_matches("c"),) => [:c1], (-if_matches('b'),) => [:a, :c1],
    (if_keys(x->length(x) == 2),) => [:c1], (~if_keys(x->length(x) == 1),) => [:c1]
]

predicate_queries = [
    (if_values(x -> !any(ismissing.(x))),) => [:a, :b, :c1],
    (-if_values(x -> !any(ismissing.(x))),) => [],
    (!if_values(x -> !any(ismissing.(x))),) => [],
    (if_eltype(Char),) => [:b],
    (~if_eltype(Int),) => [:b, :c1],
    (if_pairs((k,v) -> k == "a" || eltype(v) == Char),) => [:a, :b],
    (!if_pairs((k,v) -> k == "a" || eltype(v) == Char),) => [:c1]
]

combined_queries = [
    (if_eltype(Int) | if_matches(r"b"),) => [:a, :b],
    (-if_eltype(Int) | if_matches(r"b"),) => [:b, :c1],
    (-if_eltype(Int) & if_matches(r"b"),) => [:b],
    (cols(:c1, (not(:a) & if_matches(r"b"))) | cols(:a),) => [:c1, :b, :a],
    (not(:c1) & not(:b),) => [:a],
    (cols(:c1) | cols(:a),) => [:c1, :a]
]

context_queries = [
    (:b, all_cols()) => [:b, :a, :c1],
    (:c1, :a, all_cols(), :b) => [:c1, :a, :b],
    (all_cols(),) => [:a, :b, :c1],
    (:b, other_cols()) => [:b, :a, :c1],
    (:c1, :a, other_cols(), :b) => [:c1, :a, :b],
    (:b, else_cols()) => [:b, :a, :c1],
    (:c1, :a, else_cols(), :b) => [:c1, :a, :b]
]

context_errors = [:(cols(:a) & all_cols()) => MethodError]

renaming_queries = [
    (all_cols() => key_map(uppercase),) => [:a => :A, :b => :B, :c1 => :C1],
    (:a => :X,) => [:a => :X],
    ((cols(:a) | cols(:b)) => key_map(uppercase),) => [:a => :A, :b => :B],
    (1:2 => [:X, :Y],) => [:a => :X, :b => :Y],
    (1:2 => key_map(uppercase),) => [:a => :A, :b => :B],
    (1 => key_prefix("A"), :b => key_suffix("B"),) => [:a => :Aa, :b => :bB],
    (:c1 => key_suffix("C"), not(:b) => key_suffix("D"),) => [:c1 => :c1CD, :a => :aD]
]

renaming_queries_warn = [
    (query = (:a => [:X, :Y],) => [:a => :a],
     msg = "Renaming array had different length (2) than target selections (1), renaming skipped."),
    (query = (1:3 => [:X, :Y],) => [:a => :a, :b => :b, :c1 => :c1],
     msg = "Renaming array had different length (2) than target selections (3), renaming skipped."),
    (query = (1:2 => [:X],) => [:a => :a, :b => :b],
     msg = "Renaming array had different length (1) than target selections (2), renaming skipped."),
    (query = (1:2 => :X,) => [:a => :X, :b => :X1],
     msg = "Following columns' renaming was modified to preserve uniqueness"),
    (query = (:a => :b, :b, :c1 => :b) => [:a => :b1, :b => :b, :c1 => :b2],
     msg = "Following columns' renaming was modified to preserve uniqueness")
]

@testset "select" begin
    for (tab, pkg) in [(df, :DataFrames), (jdb, :JuliaDB), (tt, :TypedTables)]
        @testset "$(pkg)" begin
            @time test_queries_and_errors("Symbol", symbol_queries, symbol_errors, tab)
            @time test_queries_and_errors("Int", int_queries, int_errors, tab)
            @time test_queries_and_errors("Int+Symbol", mixed_queries, mixed_errors, tab)
            @time test_queries_and_errors("IntRange", int_range_queries, int_range_errors, tab)
            @time test_queries_and_errors("SymbolRange", symbol_range_queries, symbol_range_errors, tab)
            @time test_queries_and_errors("Bool", bool_queries, bool_errors, tab)
            @time test_queries("Matches", regex_queries, tab)
            @time test_queries("Predicate", predicate_queries, tab)
            @time test_queries("Chaining", combined_queries, tab)
            @time test_queries_and_errors("Complement", context_queries, context_errors, tab)
            @time test_renaming_queries("Renaming", renaming_queries, tab)
            @time test_renaming_queries_warn("Renaming Warning", renaming_queries_warn, tab)
            @time test_transforms("transform", tab)
            @time test_transforms("transform bycol", bycol(x -> x .+ 1), tab)
            @time test_transforms("transform bycol!", bycol!(x -> x .+ 1), tab)
            @time test_transforms("transform bycol.", bycol.(x -> x + 1), tab)
            @time test_transforms("transform bycol!.", bycol!.(x -> x + 1), tab)
            @time test_transforms("transform byrow", byrow((x, c) -> getindex.(x, c) .+ 1), tab)
            @time test_transforms("transform byrow!", byrow!((x, c) -> getindex.(x, c) .+ 1), tab)
            @time test_transforms("transform byrow.", byrow.((x, c) -> x[c] + 1), tab)
            @time test_transforms("transform byrow!.", byrow!.((x, c) -> x[c] + 1), tab)
            @time test_transforms("transform bytab", bytab((x, c) -> x[c] .+ 1), tab)
            @time test_transforms("transform bytab!", bytab!((x, c) -> x[c] .+ 1), tab)
            @time test_newcols("add new column", tab)
        end
    end
end


text:
## Assumptions
## http://docs.sympy.org/0.7.2/modules/assumptions/index.html

"""
     refine

Simplify an expression using assumptions; [refine](https://docs.sympy.org/dev/modules/assumptions/refine.html).
"""
refine(ex, assumpts...) = sympy.refine(ex, assumpts...)

"""
    ask(query)

Returns `true`, `false`, or `nothing`; [ask](https://docs.sympy.org/dev/modules/assumptions/ask.html)

Example:

```jldoctest ask
julia> using SymPyCall

julia> @syms x::integer y::integer
(x, y)

julia> ask(𝑄.integer(x*y), And(𝑄.integer(x), 𝑄.integer(y)))
ERROR: UndefVarError: And not defined
Stacktrace:
 [1] top-level scope
   @ none:1

julia> ## really slow isprime:
       filter(x -> ask(𝑄.prime(x)), 1:10)
4-element Vector{Int64}:
 2
 3
 5
 7
```

"""
function ask(x::Sym, args...)
    val = sympy.ask(x, args...)
    Bool(getpy(val) == pybuiltins.None) && return nothing
    Bool(getpy(val) == pybuiltins.True) && return true
    Bool(getpy(val) == pybuiltins.False) && return false
    throw(DomainError("Argument not none, true, or false"))
end

ask(x::Bool, args...) = x
ask(x::Nothing, args...) = x

# ## should we support & and | for (sym,sym) pairs? Not sure
# ## dependso on what happens to | (x==0) and ex == x usage for ==
# ## for now, we can combine terms logically with And, Or, Not...
# ## these are in logic module

# ## We make a module Q to hold the assumptions
# ## this follows this page http://docs.sympy.org/0.7.5/_modules/sympy/assumptions/ask.html
# """
#     𝑄 [\itQ]
#     SymPy.Q

# Documentation for the `SymPy.Q` module, exported as `𝑄`.

# SymPy allows for
# [assumptions](https://docs.sympy.org/latest/modules/assumptions/index.html)
# on variables. These may be placed on free sympols at construction.

# For example, the following creates a real value variable `x` and a postive, real variable `y`:

# ```jldoctest 𝑄
# julia> using SymPy

# julia> @syms x::real
# (x,)

# julia> @syms y::(real,positive)
# (y,)
# ```

# The `Q` module exposes a means to *q*uery the assumptions on a
# variable. For example,

# ```jldoctest 𝑄
# julia> ask(𝑄.positive(y))  # true
# true

# julia> ask(𝑄.negative(y))  # false
# false

# julia> ask(SymPyCall.Q.positive(x))  # `nothing`

# julia> ask(SymPyCall..Q.positive(x^2)) # `nothing` -- might be 0

# julia> ask(SymPyCall.Q.positive(1 + x^2)) # true  -- must be postive now.
# true
# ```

# The ask function uses tri-state logic, returning one of 3 values:
# `true`; `false`; or `nothing`, when the query is indeterminate.

# The construction of predicates is done through `Q` methods. These can
# be combined logically. For example, this will be `true`:

# ```jldoctest 𝑄
# julia> ask(𝑄.positive(y) & 𝑄.negative(-x^2 - 1))

# ```

# The above use `&` as an infix operation for the binary operator
# `And`. Values can also be combined with `Or`, `Not`, `Xor`, `Nand`,
# `Nor`, `Implies`, `Equivalent`, and `satisfiable`.

# !!! note "Matrix predicates"
#     As `SymPy.jl` converts symbolic matrices into Julia's `Array`
# type and not as matrices within Python, the predicate functions from SymPy for
# matrices are not used, though a replacement is given.
# """
# module Q
# import SymPyCall: sympy
# import PythonCall
# import LinearAlgebra: det, norm

# ##http://docs.sympy.org/dev/_modules/sympy/assumptions/ask.html#ask
# Q_predicates = (:antihermitian,
#                 :bounded, :finite, # bounded deprecated
#                 :commutative,
#                 :complex,
#                 :composite,
#                 :even,
#                 :extended_real,
#                 :hermitian,
#                 :imaginary,
#                 :infinitesimal,
#                 :infinity, :infinite, # :infinity deprecated
#                 :integer,
#                 :irrational,
#                 :rational,
#                 :algebraic,
#                 :transcendental,
#                 :negative,
#                 :nonzero, :zero,
#                 :positive,
#                 :prime,
#                 :real,
#                 :odd,
#                 :is_true,
#                 :nonpositive,
#                 :nonnegative
# #                :symmetric,
# #                :invertible,
# #                :singular,
# #                :orthogonal,
# #                :unitary,
# #                :normal,
# #                :positive_definite,
# #                :upper_triangular,
# #                :lower_triangular,
# #                :diagonal,
# #                :triangular,
# #                :unit_triangular,
# #                :fullrank,
# #                :square,
# #                :real_elements,
# #                :complex_elements,
# #                :integer_elements
# )

# for meth in Q_predicates
#    nm = string(meth)
#       @eval begin
#           ($meth)(x) = sympy.Q.$(nm)(x)
#               #PyCall.pycall(SymPy.sympy.Q.$nm, SymPy.Sym, x)::SymPy.Sym
#    end
# end

# symmetric(M::Array{T,2}) where {T <: SymPyCall.Sym} = SymPyCall.issymmetric(M)
# function invertible(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     d = det(M)
#     pos = SymPyCall.ask(positive(d))
#     if pos == nothing
#         return nothing
#     elseif pos == true
#         return true
#     end
#     neg = SymPyCall.ask(negative(d))
#     if neg == nothing
#         return nothing
#     end
#     z = SymPyCall.ask(zero(d))
#     if z == nothing
#         return nothing
#     elseif z == true
#         return false
#     end

#     return true
# end

# function singular(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     !invertible(M)
# end

# function orthogonal(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     vals = SymPyCall.simplify.(SymPyCall.simplify.(M*transpose(M)) .== one(T))
#     no_nothing = 0
#     for val in vals
#         a = SymPyCall.ask(zero(val))
#         if a == nothing
#             no_nothing += 1
#         elseif a == false
#             return false
#         end
#     end

#     no_nothing > 0 && return nothing
#     return true
# end

# function unitary(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     vals = SymPyCall.simplify.(SymPyCall.simplify.(M*ctranspose(M)) .== one(T))
#     no_nothing = 0
#     for val in vals
#         a = SymPyCall.ask(zero(val))
#         if a == nothing
#             no_nothing += 1
#         elseif a == false
#             return false
#         end
#     end

#     no_nothing > 0 && return nothing
#     return true
# end

# function normal(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     lhs = ctranspose(M) * M
#     rhs = M * ctranspose(M)
#     vals = zero.(SymPyCall.simplify.(lhs - rhs))
#     no_nothing = 0
#     for val in vals
#         a = SymPyCall.ask(val)
#         if a == nothing
#             no_nothing += 1
#         elseif a == false
#             return false
#         end
#     end

#     no_nothing > 0 && return nothing
#     return true
# end

# # Use [Sylvester's](https://en.wikipedia.org/wiki/Sylvester%27s_criterion) Criteria
# function positive_definite(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     !SymPyCall.ask(square(M)) && return false
#     !SymPyCall.ask(symmetric(M))  && return false
#     m, n = size(M)
#     no_false = 0
#     no_nothing = 0
#     for i in 1:m
#         a = SymPyCall.ask(Q.positive(det(M[1:i, 1:i])))
#         if a == nothing no_nothing += 1 end
#         if a == false no_false += 1 end
#     end
#     if no_false > 0
#         return false
#     elseif no_nothing > 0
#         return nothing
#     else
#         return true
#     end
# end

# upper_triangular(M::Array{T,2}) where {T <: SymPyCall.Sym} = SymPyCall.istriu(M)
# lower_triangular(M::Array{T,2}) where {T <: SymPyCall.Sym} = SymPyCall.istril(M)
# diagonal(M::Array{T,2}) where {T <: SymPyCall.Sym} = upper_triangular(M) && lower_triangular(M)
# triangular(M::Array{T,2}) where {T <: SymPyCall.Sym} = upper_triangular(M) || lower_triangular(M)

# ## This is likely not the best way as it is a bit fidgety due
# ## to the call to rref.
# function full_rank(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     m,n = size(M)
#     m <= n || return full_rank(transpose(M))

#     rr, p = SymPyCall.rref(M)
#     lr = rr[end, :] # is this zero?
#     no_nothing = 0
#     no_nonzero = 0
#     for val in lr
#         a = SymPyCall.ask(nonzero(val))
#         if a == nothing
#             no_nothing += 1
#         end
#         if a == true
#             no_nonzero += 1
#         end
#     end
#     if no_nothing > 0
#         return nothing
#     elseif no_nonzero == 0
#         return false
#     else
#         return true
#     end

# end

# function square(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     m,n = SymPyCall.size(M)
#     m == n
# end

# function real_elements(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     vals = real.(M)
#     for val in vals
#         a = SymPyCall.ask(real(val))
#         (a == nothing || a == false) && return false
#     end
#     return true
# end

# function complex_elements(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     vals = real.(M)
#     for val in vals
#         a = SymPyCall.ask(SymPyCall.sympy."Q".complex(val))
#         (a == nothing || a == false) && return false
#     end
#     return true

# end

# function integer_elements(M::Array{T,2}) where {T <: SymPyCall.Sym}
#     vals = real.(M)
#     for val in vals
#         a = SymPyCall.ask(integer(val))
#         (a == nothing || a == false) && return false
#     end
#     return true
# end

# end

# ## Issue  #354; request to *not*  export  Q
# ## export
# #export Q

# const 𝑄 = Q

# """
#     𝑄

# Exported  symbol for  [`SymPy.Q`](@ref), a  Julia  module implementing `sympy.Q`. "Questions" can be asked through the patterns
# `𝑄.query(value)` (𝑄 is entered as  [slash]itQ[tab]) or `SymPy.Q.query(value)` *but  not* as `sympy.Q.query(value)`

# !!! note
#     At one time, the symbol `Q` was exported for this. To avoid namespace clutter, the unicode alternative is now used. Legacy code would need a definition like `const Q = SymPy.Q`  to work.

# """
# 𝑄
# export  𝑄

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
const _type_pairs = [
    ("f64" => Float64),
    ("f32" => Float32),
    ("f16" => Float16),
    ("i64" => Int64),
    ("i32" => Int32),
    ("i16" => Int16),
    ("i8"  => Int8),
    ("u64" => UInt64),
    ("u32" => UInt32),
    ("u16" => UInt16),
    ("u8"  => UInt8)
]

const _type_map_jl2sigmf = Dict(reverse.(_type_pairs))
const _type_map_sigmf2jl = Dict(_type_pairs)

function datatype_julia2sigmf(::Type{Complex{T}}) where T
    return string("c", _type_map_jl2sigmf[T], "_le")
end

"""
    datatype_sigmf2julia(::String) -> Type, Bool

Return a tuple of: the Julia representation of a SigMF datatype, and whether format is little endian (true)
"""
function datatype_sigmf2julia(s::String)
    #dataset-format = (real / complex) ((type endianness) / byte)
    #real = "r"
    #complex = "c"
    #type = floating-point / signed-integer / unsigned-integer
    #floating-point = "f32"
    #signed-integer = "i32" / "i16"
    #unsigned-integer = "u32" / "u16"
    #endianness = little-endian / big-endian
    #little-endian = "_le"
    #big-endian = "_be"
    #byte = "i8" / "u8"

    # get the numeric repr of the type
    e = findnext('_', s, 1)
    isnothing(e) && throw(SigMFDataTypeError(s))

    # 'real' seems to be default if unspecified
    repr_start = s[1] in ('r', 'c') ? 2 :  s[1] in ('f', 'i', 'u') ? 1 : throw(SigMFDataTypeError(s))

    numeric_repr = s[repr_start:e-1]
    haskey(_type_map_sigmf2jl, numeric_repr) || throw(SigMFDataTypeError(s))

    # endianness
    little_endian = s[e+1:end] == "le" ? true : false

    # return either complex/real depending on first character
    if s[1] == 'r' || repr_start == 1
        return _type_map_sigmf2jl[numeric_repr], little_endian
    elseif s[1] == 'c'
        return Complex{_type_map_sigmf2jl[numeric_repr]}, little_endian
    else
        throw(SigMFDataTypeError(s))
    end
end


text:
@inline function ExaTronKernels.dgpnorm(n::Int, x::ROCDeviceArray{Float64,1}, xl::ROCDeviceArray{Float64,1},
                         xu::ROCDeviceArray{Float64,1}, g::ROCDeviceArray{Float64,1})
    tx = workitemIdx().x

    inf_norm = @amdlocalmem(Float64, 1)
    AMDGPU.sync_workgroup()
    res = 0.0

    v = 0.0
    if tx == 1
        inf_norm[1] = 0.0
        for i in 1:n
            @inbounds begin
                if xl[i] != xu[i]
                    if x[i] == xl[i]
                        v = min(g[i], 0.0)
                        v = v*v
                    elseif x[i] == xu[i]
                        v = max(g[i], 0.0)
                        v = v*v
                    else
                        v = g[i]*g[i]
                    end

                    v = sqrt(v)
                    if inf_norm[1] > v
                        inf_norm[1] = inf_norm[1]
                    else
                        inf_norm[1] = v
                    end
                end
            end
        end
    end

    AMDGPU.sync_workgroup()
    res = inf_norm[1]
    AMDGPU.sync_workgroup()
    return res
end


text:
"""
    replica_fidelity(df::DataFrame; p_field = :hproj, skip = 0)
Compute the fidelity of the average coefficient vector and the projector defined in
`p_field` from the result of replica [`lomc!()`](@ref) passed as  argument `df`,
using replicas `_1` and `_2`.
Calls [`ratio_of_means()`](@ref) to perform a blocking analysis
on a ratio of the means of separate time series and returns a
[`RatioBlockingResult`](@ref).
The first `skip` steps in the time series are skipped.

The fidelity of states `|ψ⟩` and `|ϕ⟩` is defined as
```math
F(ψ,ϕ) = \\frac{|⟨ψ|ϕ⟩|^2}{⟨ψ|ψ⟩⟨ϕ|ϕ⟩} .
```
Specifically, `replica_fidelity` computes
```math
F(\\mathbf{v},⟨\\mathbf{c}⟩) =
    \\frac{⟨(\\mathbf{c}_1⋅\\mathbf{v})(\\mathbf{v}⋅\\mathbf{c}_1)⟩}
    {⟨\\mathbf{c}_1⋅\\mathbf{c}_1⟩} ,
```
where `v` is the projector specified by `p_field`, which is assumed to be normalised to
unity with the two-norm (i.e. `v⋅v == 1`), and ``\\mathbf{c}_1`` and ``\\mathbf{c}_2``
are two replica coefficient vectors.
"""
function replica_fidelity(df::DataFrame; p_field = :hproj, skip = 0,  args...)
    p_field_1 = Symbol(p_field, :_1)
    p_field_2 = Symbol(p_field, :_2)
    fid_num = conj(getproperty(df, p_field_1)) .* getproperty(df, p_field_2)
    fid_num = fid_num[skip+1:end]
    # denominator
    fid_den = df.c1_dot_c2[skip+1:end]

    return ratio_of_means(fid_num, fid_den; args...)
end


text:
using PastaQ
using PastaQ.ITensors
using Test
using LinearAlgebra
using Random

@testset "noise models: pauli channel" begin
  E = gate("pauli_channel")
  K = [E[:,:,k] for k in 1:size(E,3)]
  @test size(E,3) == 4
  @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,2,2)
  for N in 2:5
    probs = rand(4^N)
    probs = probs ./ sum(probs)
    E = gate("pauli_channel", Tuple(repeat([2],N)); error_probabilities = probs)
    @test size(E,3) == 4^N
    @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,1<<N,1<<N)
  end
end

@testset "noise models: bitflip" begin
  E = gate("bit_flip"; p = 0.1)
  K = [E[:,:,k] for k in 1:size(E,3)]
  @test size(E,3) == 2
  @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,2,2)
  for N in 2:5
    E = gate("bit_flip", Tuple(repeat([2],N)); p = 0.1)
    @test size(E,3) == 2^N
    @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,1<<N,1<<N)
  end
end

@testset "noise models: phaseflip" begin
  E = gate("phase_flip"; p = 0.1)
  K = [E[:,:,k] for k in 1:size(E,3)]
  @test size(E,3) == 2
  @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,2,2)
  for N in 2:5
    E = gate("phase_flip", Tuple(repeat([2],N)); p = 0.1)
    @test size(E,3) == 2^N
    @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,1<<N,1<<N)
  end
end

@testset "noise models: amplitude damping" begin
  E = gate("AD"; γ = 0.1)
  K = [E[:,:,k] for k in 1:size(E,3)]
  @test size(E,3) == 2
  @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,2,2)
  for N in 2:5
    E = gate("AD", Tuple(repeat([2],N)); γ = 0.1)
    @test size(E,3) == 2^N
    @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,1<<N,1<<N)
  end

  N = 2
  circuit1 = randomcircuit(N; depth = 3, layered = false)
  circuit2 = copy(circuit1)
  push!(circuit1, ("AD",1,(γ=0.1,)))
  push!(circuit1, ("AD",2,(γ=0.1,)))
  push!(circuit2, ("AD",(1,2),(γ=0.1,)))
  ρ₀ = MPO(productstate(N))
  ρ1 = runcircuit(ρ₀, circuit1)
  ρ2 = runcircuit(ρ₀, circuit2)
  @test PastaQ.array(ρ1) ≈ PastaQ.array(ρ2)
end

@testset "noise models: phase damping" begin
  E = gate("PD"; γ = 0.1)
  K = [E[:,:,k] for k in 1:size(E,3)]
  @test size(E,3) == 2
  @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,2,2)
  for N in 2:5
    E = gate("PD", Tuple(repeat([2],N)); γ = 0.1)
    @test size(E,3) == 2^N
    @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,1<<N,1<<N)
  end

  N = 2
  circuit1 = randomcircuit(N; depth = 3, layered = false)
  circuit2 = copy(circuit1)
  push!(circuit1, ("PD",1,(γ=0.1,)))
  push!(circuit1, ("PD",2,(γ=0.1,)))
  push!(circuit2, ("PD",(1,2),(γ=0.1,)))
  ρ₀ = MPO(productstate(N))
  ρ1 = runcircuit(ρ₀, circuit1)
  ρ2 = runcircuit(ρ₀, circuit2)
  @test PastaQ.array(ρ1) ≈ PastaQ.array(ρ2)
end

@testset "noise models: depolarizing" begin
  E = gate("DEP"; p = 0.1)
  K = [E[:,:,k] for k in 1:size(E,3)]
  @test size(E,3) == 4
  @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,2,2)
  for N in 2:5
    E = gate("DEP", Tuple(repeat([2],N)); p = 0.1)
    @test size(E,3) == 4^N
    @test sum([E[:,:,k]' * E[:,:,k] for k in 1:size(E,3)]) ≈ Matrix{Float64}(I,1<<N,1<<N)
  end
end

@testset "insertnoise: iid noise" begin
  N = 6
  depth = 4

  circuit = randomcircuit(N; depth = depth, twoqubitgates = ["CX","CZ"], onequbitgates = ["Rn","X"], layered = false)
  ngates = length(circuit)
  nCZ,nCX,nR,nX = 0,0,0,0
  for g in circuit
    g[1] == "CX" && (nCX +=1)
    g[1] == "CZ" && (nCZ +=1)
    g[1] == "Rn" && (nR +=1)
    g[1] == "X"  && (nX +=1)
  end

  noisemodel = ("DEP", (p=0.01,))
  noisycircuit = insertnoise(circuit, noisemodel)
  @test length(noisycircuit) == 2*ngates
  for k in 2:2:length(circuit)
    g = noisycircuit[k]
    @test g[1] == noisemodel[1]
    @test g[3] == noisemodel[2]
  end

  noisycircuit = insertnoise(circuit, noisemodel; gate = "CX")
  @test length(noisycircuit) == ngates + nCX
  for k in 1:length(circuit)
    g = noisycircuit[k]
    if g[1] == "CX"
      @test noisycircuit[k+1][1] == noisemodel[1]
      @test noisycircuit[k+1][3] == noisemodel[2]
    else
      @test noisycircuit[k+1][1] != noisemodel[1]
    end
  end
  noisycircuit = insertnoise(circuit, noisemodel; gate = ["CX","CZ"])
  @test length(noisycircuit) == ngates + nCX + nCZ

end

@testset "insertnoise: one and two qubit noise" begin
  N = 6
  depth = 4

  circuit = randomcircuit(N; depth = depth, twoqubitgates = ["CX","CZ"], onequbitgates = ["Rn","X"], layered = false)
  ngates = length(circuit)
  nCZ,nCX,nR,nX = 0,0,0,0
  for g in circuit
    g[1] == "CX" && (nCX +=1)
    g[1] == "CZ" && (nCZ +=1)
    g[1] == "Rn" && (nR +=1)
    g[1] == "X"  && (nX +=1)
  end

  noise1Q = ("AD", (γ=0.01,))
  noise2Q = ("DEP", (p=0.05,))

  noisemodel = (1 => noise1Q, 2 => noise2Q)
  noisycircuit = insertnoise(circuit, noisemodel)
  @test length(noisycircuit) == 2*ngates
  for k in 1:2:length(circuit)
    g = noisycircuit[k]
    if g[2] isa Int
      @test noisycircuit[k+1][1] == noise1Q[1]
      @test noisycircuit[k+1][3] == noise1Q[2]
    else
      @test noisycircuit[k+1][1] == noise2Q[1]
      @test noisycircuit[k+1][3] == noise2Q[2]
    end
  end
end


text:
using Ansillary

@enum Direction Up Down Left Right

mutable struct Game
    fruit::Cursor.Coordinate
    snake::Vector{Cursor.Coordinate}
    clear::Union{Nothing, Cursor.Coordinate}
    direction::Direction
    finished::Bool
end

function Game(rows, columns)
    crow = div(rows - 3, 2)
    ccol = div(columns - 2, 2)
    snake = [
        Cursor.Coordinate(crow, ccol - 2),
        Cursor.Coordinate(crow, ccol - 1),
        Cursor.Coordinate(crow, ccol),
        Cursor.Coordinate(crow, ccol + 1),
    ]

    fruit = new_fruit(snake, rows, columns)

    Game(fruit, snake, nothing, Right, false)
end

function draw(rows, columns)
    Cursor.move!(Cursor.Coordinate(1, 1))

    print("┏", "━"^(columns - 2), "┓")

    for row in 2:(rows - 2)
        Cursor.move!(Cursor.Coordinate(row, 1))
        print("┃")

        Cursor.move!(Cursor.Column(columns))
        print("┃")
    end

    Cursor.move!(Cursor.Coordinate(rows - 1, 1))
    print("┗", "━"^(columns - 2), "┛")

    Cursor.move!(Cursor.Coordinate(rows, 1))
    print("Arrow Keys: Change direction. - Esc: Quit.")
end

function draw(game::Game)
    if game.clear !== nothing
        Cursor.move!(game.clear)
        print(' ')
    end

    draw(game.snake)
    draw(game.fruit)
end

function draw(snake::Vector{Cursor.Coordinate})
    char = 'O'
    for coord in Iterators.reverse(snake)
        Cursor.move!(coord)
        print(char)
        char = '*'
    end
end

function draw(fruit::Cursor.Coordinate)
    Cursor.move!(fruit)
    print('$')
end

function update!(game, rows, columns)
    game.finished && return

    old = game.snake[end]

    new = if game.direction == Up
        Cursor.Coordinate(old.row - 1, old.column)
    elseif game.direction == Down
        Cursor.Coordinate(old.row + 1, old.column)
    elseif game.direction == Left
        Cursor.Coordinate(old.row, old.column - 1)
    elseif game.direction == Right
        Cursor.Coordinate(old.row, old.column + 1)
    end

    if new.row == 1
        new = Cursor.Coordinate(rows - 2, new.column)
    end

    if new.row == rows - 1
        new = Cursor.Coordinate(2, new.column)
    end

    if new.column == 1
        new = Cursor.Coordinate(new.row, columns - 1)
    end

    if new.column == columns
        new = Cursor.Coordinate(new.row, 2)
    end

    if new in game.snake
        game.finished = true
        return
    end

    push!(game.snake, new)

    if new == game.fruit
        game.fruit = new_fruit(game.snake, rows, columns)
        game.clear = nothing
    else
        game.clear = popfirst!(game.snake)
    end
end

function new_fruit(snake, rows, columns)
    while true
        fruit = Cursor.Coordinate(rand(2:(rows - 2)), rand(2:(columns - 1)))
        if !(fruit in snake)
            return fruit
        end
    end
end

function main()
    size = Screen.raw(Screen.size)
    main(size.row, size.column)
end

function main(rows, columns)
    Screen.alternative() do
        Cursor.hide() do
            game = Game(rows, columns)

            draw(rows, columns)
            draw(game)

            for event in Inputs.EventLoop(Inputs.Millisecond(200))
                if event == Inputs.Esc()
                    break
                elseif event == Inputs.Up()
                    game.direction = Up
                elseif event == Inputs.Down()
                    game.direction = Down
                elseif event == Inputs.Left()
                    game.direction = Left
                elseif event == Inputs.Right()
                    game.direction = Right
                elseif event == Inputs.Tick()
                    update!(game, rows, columns)
                    draw(game)
                end
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    function usage(code = 0)
        println(PROGRAM_FILE, " [-h|--help] [ROWS COLUMNS]")
        exit(code)
    end

    if "-h" in ARGS || "--help" in ARGS
        usage()
    elseif isempty(ARGS)
        main()
    elseif length(ARGS) == 2
        main(parse(UInt16, ARGS[1]), parse(UInt16, ARGS[2]))
    else
        usage(1)
    end
end


text:
include("header.jl")

@testset "sparse" begin

    # Issue #114 (a): plus for Sparse.
    using AutoGrad: Sparse, full, addto!
    a = Sparse(zeros(3,4), [ [1.,1.], [1.,1.], 1., 1. ], [ ([1,2],), (3:4,), (2,2), (1,) ])
    b = a + a
    @test b isa Sparse
    @test full(b) == full(a) + full(a)
    addto!(a, a)
    @test a isa Sparse
    @test full(a) == full(b)
    b = a + a
    a .+= a
    @test full(a) == full(b)

    # Issue #114 (b): lmul! ambiguous for Sparse, breaks gclip.
    using LinearAlgebra
    foo(w) = (s = 0.0; for i=1:length(w); s+=w[i]; end; s)
    w = Param(randn(2,2))
    J = @diff foo(w)
    @test lmul!(0.5, grad(J,w)) isa Sparse

end


text:
module ZirconAges

# Julia translation of a R-script provided by Gregor Weber & Luca Caricchi to compute zircon ages from Tt paths
# Used in the Nat. Comm. publication "Estimating the current size and state of subvolcanic magma reservoirs"
# 15/04/2022, Nico Riel & Boris Kaus

import Base.Threads
using Parameters
using Loess, Statistics, StatsBase, KernelDensity, Loess

export  ZirconAgeData,
        compute_zircon_age_PDF,  compute_zircons_Ttpath,  # calculation routines
        zircon_age_PDF,
        compute_zircons_convert_vecs2mat

"""
    Declare function to calculate zircon fraction as function of a temperature profile
"""
function zircon_fraction(T::AbstractArray{_T}, max_x_zr::_T) 	where _T
    A = (1.62.-1.8*(10^4)*exp.((-10^4)./(T .+ 273.15))).*max_x_zr
    A[A .<= 0.0] .= 0.0

    return A
end

"""
    ZirconAgeData

Struct that holds default parameters for the calculations
"""
@with_kw_noshow struct ZirconAgeData
    Tsat::Float64 			= 825.0		# Maximum zircon saturation temperature [C]
    Tmin::Float64 			= 690.0 	# [C] Minimum zircon saturation Temperature [C]
    Tsol::Float64 			= 690.0		# [C] Solidus temperature [C]
    Tcal_max::Float64		= 800.0		# max temperature to calculate zircon fraction
    Tcal_step::Float64 		= 1.0		# temperature step to caclulate zircon fraction (resolution of Zircon saturation curve discretization)
    max_x_zr::Float64 		= 0.001		# max fraction zircons at solidus
    zircon_number::Int64    = 100.0		# number of required zircons
    time_zr_growth::Float64 = 0.7e6		# Minimum time within T saturation range (This is what the method used in the R script, boils down too)
                                        # -> remain in the Zr saturation zone more than 1/3 of the time the Tt path with the longest time in the saturation zone
end

"""
    Loess fit of zircon number over saturation range
"""
function loess_fit_zircon_sat(ZirconData::ZirconAgeData)
    @unpack Tsol, Tcal_step, Tcal_max, Tsat, max_x_zr, zircon_number, Tmin = ZirconData

    # Get cumulative Zircon fraction (after Tierney et al., 2016; Geology)
    T 				= range(Tsol, stop = Tcal_max, step = Tcal_step  )
    Tfit 			= range(Tsol, stop = Tsat, 	   length = length(T))

    x_zircon 		= zircon_fraction(T,max_x_zr)			# compute zircon fraction
    n_zircon 		= zero(x_zircon)						# number of zircons

    n_zircon[2:end] = -diff(x_zircon)						# differentiate zircon fraction
    n_zircon[1] 	= n_zircon[2]							# set first value

    n_zircon		= ceil.(( (n_zircon*zircon_number)/maximum(n_zircon)).-minimum(floor.((n_zircon*zircon_number)/maximum(n_zircon))))

    # fit n_zircon_N with loess regression
    n_zircon_N_fit 	= loess(Tfit, n_zircon, span=1.0)

    return n_zircon_N_fit
end

"""

This employs a loess fit function to compute the number of zircons for each of the Tt paths
"""
function  compute_number_zircons!(n_zr::AbstractArray{_T,N}, Tt_paths_Temp::AbstractArray{_T,N}, ZirconData::ZirconAgeData) where {_T,N}
    @unpack Tmin, Tsat, Tsol = ZirconData

    n_zircon_N_fit  = 	loess_fit_zircon_sat(ZirconData)			# loess fit through zircon saturation
    Threads.@threads for i in 1:size(Tt_paths_Temp,2)
        for j in 1:size(Tt_paths_Temp,1)
            T = Tt_paths_Temp[j,i]
            if (T > Tsol) & (T < Tsat)
                # The following line takes a lot of time, because Loess is not type-stable.
                # There is an open PR in the package that may fix it
                dat::_T      = Loess.predict(n_zircon_N_fit, T)
                n_zr[j,i]    = floor(dat)
            else
                n_zr[j,i]    = 0.0
            end
        end
    end

    nothing
end

"""
    prob, ages_eruptible, number_zircons, T_av_time, T_sd_time  compute_zircons_Ttpath(time_years::AbstractArray{Float64,1}, Tt_paths_Temp::AbstractArray{Float64,2}; ZirconData::ZirconAgeData)

This computes the number of zircons produced from a series of temperature-time path's.
The Tt-paths are stored in a 2D matrix `Tt_paths_Temp` with rows being the temperature at time `time_years`.

Input:
====
- `time_years` : vector of length `nt` with the time in years (since the beginning of the simulation) of the points provided
- `Tt_paths_Temp` : array of size `(nt,npaths)`` with the temperature of every path.

Output:
- `prob` : a vector that gives the relative probability that a zircon with a given age exists
- `ages_eruptible` : age of eruptble magma
- `number_zircons` : 2D array of size `(nt,)`
- `T_av_time`: vector of size `nt` that contains the average T of the paths
- `T_sd_time`: vector of size `nt` that contains the standard deviation of the T of the paths

This routine was developed based on an R-routine provided as electronic supplement in the paper:
- Weber, G., Caricchi, L., Arce, J.L., Schmitt, A.K., 2020. Determining the current size and state of subvolcanic magma reservoirs. Nat Commun 11, 5477. https://doi.org/10.1038/s41467-020-19084-2

"""
function compute_zircons_Ttpath(time_years::AbstractArray{_T,1}, Tt_paths_Temp::AbstractArray{_T,2}; ZirconData::ZirconAgeData = ZirconAgeData()) where _T

    @unpack Tmin, Tsat, Tsol, time_zr_growth = ZirconData

    Tt_paths_Temp1 = copy(Tt_paths_Temp)

    Δt 				=	diff(time_years)[1]							# timestep [yrs]
    time_er_min 	= 	maximum(time_years)							# backward count

    # find all the Tt paths that go through the zircon saturation range & are at the end of the path still below Tsat (otherwise Zr are not yet crystallized)
    ID_col_er 		= findall( (maximum(Tt_paths_Temp1,dims=1).>Tmin) .& (Tt_paths_Temp1[end,:]' .< Tsat))
    n_zr			= zero(Tt_paths_Temp1)
    compute_number_zircons!(n_zr, Tt_paths_Temp1, ZirconData)		# computes the number of zircons for every path

    # find the number of timesteps for every path, during which the temperature is > Tmin and < Tsat
    length_trace 	= Vector{Float64}(undef,length(ID_col_er))
    for i in 1:length(ID_col_er)
        length_trace[i]    = length( findall( (Tt_paths_Temp1[:,ID_col_er[i][2]] .> Tmin) .& (Tt_paths_Temp1[:,ID_col_er[i][2]] .< Tsat)) )
    end

    # the next several lines can likely be achieved in a more elegant way...
    id				= findall( length_trace .== maximum(length_trace))
    ID_col_lgst_tr 	= ID_col_er[ id[1] ][2]

    id 				= findall( Tt_paths_Temp1[:,ID_col_lgst_tr] .< Tsat)
    VALID_min_time 	= findmin( Tt_paths_Temp1[id,ID_col_lgst_tr])
    ID_min_time		= VALID_min_time[2]

    max_age_spread	= maximum(length_trace)*Δt					# This is defined among all selected paths

    T_av_time_1 	= replace!(Tt_paths_Temp1, 0.0 => NaN)
    T_av_time 		= zeros(size(T_av_time_1,1))
    T_sd_time 		= zeros(size(T_av_time_1,1))

    # get the average temperature of the Tt paths and the standard deviation
    for i in 1:size(Tt_paths_Temp1,1)
        T_av_time[i]	= mean(filter(!isnan, T_av_time_1[i,:]))
        T_sd_time[i]	= std(filter( !isnan, T_av_time_1[i,:]))
    end

    # I clarified the R function because the minimum step length to grow a zircon is simply a ratio of the maximum trace between Tmin and Tsol
    # This makes sense as we only deal with fractions here. Because no mass is provided the real zircon size cannot possibly be determined
    min_step_n		= floor( (time_zr_growth/max_age_spread)*(max_age_spread/Δt) )

    # find the Tt paths that have a number of timesteps in the saturation range greater than the defined min_step_n ()
    # this is to mimic that it takes some time to grow zircons
    id				= findall( length_trace .> min_step_n)
    if isempty(id)
        max_Ptpath = maximum(length_trace)*Δt
        error("I don't have a single Pt-path that is sufficiently long within time_zr_growth (=$(time_zr_growth) yrs).
            The Longest Pt-path I have is $(max_Ptpath) years.
            Decrease this value within the ZirconDataAge struct with ZirconData=ZirconAgeData(time_zr_growth=0.1e6) & rerun.")
    end
    ID_col_er_1		= getindex.(ID_col_er[id], [2])

    int_zr_sat		= collect(Float64,  ID_min_time:1.0:(time_er_min/Δt)-min_step_n)
    int_zr_sat		= floor.(Int64,int_zr_sat)

    T_av_time_slct 	= Vector{Float64}(undef,length(int_zr_sat)-1) .= 0.0

    for i in 1:(size(int_zr_sat,1)-1)
        id2				= ID_col_er_1[ findall( (Tt_paths_Temp1[int_zr_sat[i],ID_col_er_1[:]] .> Tmin) .& (Tt_paths_Temp1[int_zr_sat[i],ID_col_er_1[:]] .< Tsat)) ]
        if isempty(id2) == true
            T_av_time_slct[i] = NaN
        else
            T_av_time_slct[i] = median(filter(!isnan, Tt_paths_Temp1[int_zr_sat[i],id2]))
        end
    end

    replace!(Tt_paths_Temp1, NaN => 0.0)
    ID_col_er			= getindex.(ID_col_er, [2])
    for i in 1:length(ID_col_er)
        ind = findall( (Tt_paths_Temp1[:,ID_col_er[i]]) .== 0.0 );
        if ~isempty(ind)
            k 				= maximum(ind)
            Tt_paths_Temp1[1:k,ID_col_er[i]] .= 0.0
        end
    end

    zr_select			= zero(Tt_paths_Temp1)
    zr_select[Tt_paths_Temp1 .> 0.0] .= 1.0
    n_zrc2_0			= zr_select.*n_zr						# filters out those Tt path that are still >Tsat @ the end
    number_zircons      = n_zrc2_0[:,ID_col_er_1];
    n_measurable_ages 	= sum(number_zircons, dims=2)
    sz 					= size(number_zircons,1)
    ages_eruptible		= collect(Float64,  1.0:Δt:sz*Δt)

    # probability that a certain zircon is sampled, dependens on how many of a given age ara available:
    prob 				= n_measurable_ages/sum(n_measurable_ages)
    prob 				= prob[:,1]

    return prob, ages_eruptible, number_zircons, T_av_time, T_sd_time
end

"""
    time_years, prob, ages_eruptible, number_zircons, T_av_time, T_sd_time  = compute_zircons_Ttpath(time_years::Vector{Vector{Float64}}, Tt_paths_Temp::Vector{Vector{Float64}}; ZirconData::ZirconAgeData = ZirconAgeData())

This accepts Vector{Vector} as input for time and temperature of each Tt-path. Here, the length of the vector can be variable between different points.

Internally, we interpolate this into a 2D matrix and a longer vector that includes all paths and a single vector with times

"""
function compute_zircons_Ttpath(time_years_vecs::Vector{Vector{_T}}, Tt_paths_Temp_vecs::Vector{Vector{_T}}; ZirconData::ZirconAgeData = ZirconAgeData()) where _T

    # convert to a vector with time & matrix with T values at every timestep
    time_years, Tt_paths_Temp = compute_zircons_convert_vecs2mat(time_years_vecs, Tt_paths_Temp_vecs)

    # call main routine
    prob, ages_eruptible, number_zircons, T_av_time, T_sd_time = compute_zircons_Ttpath(time_years, Tt_paths_Temp, ZirconData=ZirconData)

    # return, including time_years
    return time_years, prob, ages_eruptible, number_zircons, T_av_time, T_sd_time

end

"""
    time_years, Ttpaths_mat = compute_zircons_convert_vecs2mat(time_years_vecs::Vector{Vector{Float64}}, Tt_paths_Temp_vecs::Vector{Vector{Float64}})

This converts a vector with Vectors contain time and temperature path's to a single time vector and a matrix that combines all vectors
"""
function compute_zircons_convert_vecs2mat(time_years_vecs::Vector{Vector{_T}}, Tt_paths_Temp_vecs::Vector{Vector{_T}}) where _T

    # Create a single vector with the time values:
    time_years = unique(reduce(vcat,unique(time_years_vecs)));

    # Add the vectors to an array with T values
    Tt_paths_mat = zeros(_T,length(time_years), length(Tt_paths_Temp_vecs))
    for i in eachindex(Tt_paths_Temp_vecs)
        istart = findall(time_years .== time_years_vecs[i][1])
        iend   = findall(time_years .== time_years_vecs[i][end])
        Tt_paths_mat[istart[1]:iend[1], i] .= Tt_paths_Temp_vecs[i]
    end

    return time_years, Tt_paths_mat

end

"""
    zircon_age_PDF(ages_eruptible::AbstractArray{Float64,1}, number_zircons::AbstractArray{Float64,2}, bandwidth=1e5, n_analyses=300)

Compute probability density functions for zircon age path's describes in `number_zircons` with age `ages_eruptible` (both computed ).
`bandwidth` is the smoothening window of the resulting curves (in years), whereas `n_analyses` are the number of analyses done.

"""
function zircon_age_PDF(ages_eruptible::AbstractArray{_T,1}, number_zircons::AbstractArray{_T,2}; bandwidth=1e5, n_analyses=300) where _T

    # compute PDF for each of the zircon Tt-paths:
    PDF_zircons = []
    time_Ma = [];
    for i in 1:size(number_zircons,2)
        n_meas 			 = number_zircons[:,i]
        px	  			 = n_meas/sum(n_meas)		# probability to have a certain age

        # random numbers selected according to the probability
        smp				 = sample( (maximum(ages_eruptible) .- ages_eruptible)/1e6, Weights(px), n_analyses, replace=true)
        y 				 = kde(smp, bandwidth=bandwidth/1e6)

        # store data
        push!(PDF_zircons, 	y.density)
        push!(time_Ma, 		y.x)
    end

    n_measurable_ages   = sum(number_zircons, dims=2)
    pxAv	  			= n_measurable_ages[:,1]./sum(n_measurable_ages[:,1])
    smpAv				= sample( (maximum(ages_eruptible) .- ages_eruptible)/1e6, Weights(pxAv), n_analyses, replace=true)
    yAv 				= kde(smpAv, bandwidth=bandwidth/1e6)
    time_Ma_average     = Vector(yAv.x);
    PDF_zircon_average  = Vector(yAv.density);

    return time_Ma, PDF_zircons, time_Ma_average, PDF_zircon_average
end

"""
    time_Ma, PDF_zircons, time_Ma_average, PDF_zircon_average, time_years, prob, ages_eruptible, number_zircons, T_av_time, T_sd_time = compute_zircon_age_PDF(time_years_vecs::Vector{Vector}, Tt_paths_Temp_vecs::Vector{Vector}; ZirconData::ZirconAgeData = ZirconAgeData(), bandwidth=bandwidth, n_analyses=300)

This computes the PDF (probability density function) with zircon age data from Vectors with Tt-paths

"""
function compute_zircon_age_PDF(time_years_vecs::Vector{Vector{_T}}, Tt_paths_Temp_vecs::Vector{Vector{_T}}; ZirconData::ZirconAgeData = ZirconAgeData(), bandwidth=1e5, n_analyses=300) where _T

    # Compute the probability that a zircon of certain age is sampled:
    time_years, prob, ages_eruptible, number_zircons, T_av_time, T_sd_time = compute_zircons_Ttpath(time_years_vecs, Tt_paths_Temp_vecs, ZirconData=ZirconData);

    # Use this to compute PDF curves:
    time_Ma, PDF_zircons, time_Ma_average, PDF_zircon_average  = zircon_age_PDF(ages_eruptible, number_zircons, bandwidth=bandwidth, n_analyses=n_analyses)

    return time_Ma, PDF_zircons, time_Ma_average, PDF_zircon_average, time_years, prob, ages_eruptible, number_zircons, T_av_time, T_sd_time

end

end


text:
################################################################################
#
#   factorbase.jl : functions for computing factor base logarithms
#
################################################################################

export factorBaseDeg2

"""
    bracket{T <: PolyElem}(D::Int , A::T, B::T, h0::T, h1::T)

Compute the `D`-bracket of `A` and `B`.

This is equal to ``h1^D(B×A(h0/h1) - A×B(h0/h1))``.
"""
function bracket{T <: PolyElem}(D::Int , A::T, B::T, h0::T, h1::T)

    # We create zero polynomials
    R = parent(A)
    Ahomo = R()
    Bhomo = R()

    # We compute h1^D × A(h0/h1)
    for j in 0:degree(A)
        Ahomo += coeff(A, j)*(h0^j)*(h1^(D-j))
    end

    # We compute h1^D × B(h0/h1)
    for j in 0:degree(B)
        Bhomo += coeff(B, j)*(h0^j)*(h1^(D-j))
    end

    # We return the result h1^D(B×A(h0/h1) - A×B(h0/h1))
    return B*Ahomo - A*Bhomo
end

"""
    fillMatrixGKZ!{T <: PolyElem}(M::MatElem, j::Int, A::T, B::T, f1::T,
                                  D::Int, fact::Nemo.Fac{T},
                                  tmpFact::Nemo.Fac{T}, unknowns::Dict{T,
                                                                       Int})

Fill the `j`-th column of the matrix M, according to the relation generated by
`Α` and `B`.
"""
function fillMatrixGKZ!{T <: PolyElem}(M::MatElem, j::Int, A::T, B::T, f1::T,
                                       D::Int, fact::Nemo.Fac{T},
                                       tmpFact::Nemo.Fac{T}, unknowns::Dict{T,
                                                                            Int})

    # We define the finite field
    F = base_ring(A)

    # We fill the matrix M to represent the factors involved in the equation
    # generated by `A` and `B`:
    #
    # B×Π_{a\in F} (A-aB) = 1/h1^D×([A, B]_D)
    # where [A, B]_D = h1^D(A(h0/h1)×B - B(h0/h1)×A)

    # We begin with the factors in the product Π (A-aB)
    for a in F
        tmp = factor(A-a*B)
        for f in tmp
            M[unknowns[f[1]],j] += f[2]
        end
    end

    # We add the factors of B
    tmp = factor(B)
    for f in tmp
        M[unknowns[f[1]],j] += f[2]
    end

    # Then the factors of h1 = X(X+t), where f1 := X+t
    M[unknowns[f1], j] += D
    M[unknowns[gen(parent(A))], j] += D

    # We add the factors on the other side of the equation
    # First the factors of h1×X-h0, indeed we always have h1×X-h0 | [A, B]_D
    for f in fact
        M[unknowns[f[1]],j] -= f[2]
    end

    # Finally the factors of [A, B]_D / (h1×X-h0)
    for f in tmpFact
        M[unknowns[f[1]],j] -= f[2]
    end
end

"""
    fillMatrixGKZ2!{T <: PolyElem}(M::MatElem, j::Int, A::T, B::T, f1::T,
                                   D::Int, tmpFact::Nemo.Fac{T},
                                   unknowns::Dict{T, Int})

Fill the matrix `M` according to the relation generated by `A` and `B`, in
the case where h1×X-h0 is irreducible.
"""
function fillMatrixGKZ2!{T <: PolyElem}(M::MatElem, j::Int, A::T, B::T, f1::T,
                                       D::Int, tmpFact::Nemo.Fac{T},
                                       unknowns::Dict{T, Int})

    # We define the finite field
    F = base_ring(A)

    # We fill the matrix M to represent the factors involved in the equation
    # generated by `A` and `B`:
    #
    # B×Π_{a\in F} (A-aB) = 1/h1^D×([A, B]_D)
    # where [A, B]_D = h1^D(A(h0/h1)×B - B(h0/h1)×A)

    # We begin with the factors in the product Π (A-aB)
    for a in F
        tmp = factor(A-a*B)
        for f in tmp
            M[unknowns[f[1]],j] += f[2]
        end
    end

    # We add the factors of B
    tmp = factor(B)
    for f in tmp
        M[unknowns[f[1]],j] += f[2]
    end

    # Then the factors of h1 = X(X+t), where f1 := X+t
    M[unknowns[f1], j] += D
    M[unknowns[gen(parent(A))], j] += D

    # We add the factors on the other side of the equation
    # First the irreducible polynomial h1×X-h0, indeed we always have h1×X-h0 | [A, B]_D
    M[length(unknowns)+1,j] -= 1

    # Finally the factors of [A, B]_D / (h1×X-h0)
    for f in tmpFact
        M[unknowns[f[1]],j] -= f[2]
    end
end

"""
    factorBaseDeg2{T <: PolyElem}(K::GkzContext{T})

Compute the logarithms of the elements of degree ≤ 2 in the field `K`.
"""
function factorBaseDeg2{T <: PolyElem}(K::GkzContext{T})

    # We setup some variables
    h0, h1 = K.h0, K.h1
    g = K.gen
    defPol = K.definingPolynomial
    D = max(degree(h0), degree(h1))
    F = base_ring(h0)
    card = ZZ(length(F))^degree(defPol)
    q::Int = length(F)
    X = gen(parent(h0))
    f1 = divexact(h1, X)
    systFactor = h1*X-h0

    # We factor the systematic factor of [A, B]_D
    fact = factor(systFactor)

    # We create a bijection between irreducible polynomials and an iteration
    # of them
    unknownsToCoord, coordToUnknowns = irreduciblesDeg2(parent(h0))
    nbukn = length(unknownsToCoord)

    # We create the dictionary containing the logarithms
    dlogs = Dict{T, Nemo.fmpz}()
    j = 0

    # We first deal with the case where the factors of the systematic factor of
    # [A, B]_D are of degree ≤ 2, hence are already in the unknowns
    if !isIrreducible(systFactor)

        # We create a zero matrix
        Sp = MatrixSpace(ZZ, nbukn, q^2)
        M = zero(Sp)

        # And we generate equations Π (A-aB) = 1/h1^D [A, B]_D
        # If [A, B]_D is smooth, we fill the matrix accordingly
        for a in F
            for b in F
                A = X^D + a
                B = X^(D-1) + b
                N = bracket(D, A, B, h0, h1)
                O = divexact(N, systFactor)

                tmpFact = factor(O)
                if isSmooth(tmpFact, 2)

                    j += 1
                    fillMatrixGKZ!(M, j, A, B, f1, D, fact, tmpFact,
                                   unknownsToCoord)

                end
            end
        end

        # Then we erase the blank columns and we transpose the matrix to later perform
        # a Gaussian elimination
        M = subMatrix(M, nbukn, j)
        M = transpose(M)

        # We start by finding the logarithms of the elements in the factor
        # base modulo the small factors, using Pohlig Hellman algorithm
        n = ZZ()

        for k in keys(unknownsToCoord)
            dlogs[k], n = pohligHellman(card, g, k, defPol)
        end

        # One we have dealt with the small factors, we compute the big ones
        v = div(card-1, n)
        bigFactors = factor(v)

        # And we compute the logarithms modulo each factor, concluding
        # using the Chinese Remaindering Theorem. With each factor, we obtain
        # a one-dimensional kernel so we have to scale by the generator.
        for p in bigFactors
            Mred = rrefMod(M, p[1])[1]
            ker = Nemo.fmpz[Mred[i, nbukn] for i in 1:(nbukn-1)]
            s = gcdinv(ker[unknownsToCoord[g]], p[1])[2]
            for j in 1:(nbukn-1)
                dlogs[coordToUnknowns[j]] =
                crt(dlogs[coordToUnknowns[j]], n, ker[j]*s, p[1])
            end
            dlogs[coordToUnknowns[nbukn]] =
            crt(dlogs[coordToUnknowns[nbukn]], n, -s, p[1])
            n *= p[1]
        end

        return dlogs

    # If the systematic factor is irreducible, we consider as an other
    # unknown
    else

        # We create a zero matrix
        Sp = MatrixSpace(ZZ, nbukn+1, q^2)
        M = zero(Sp)

        # And we generate equations Π (A-aB) = 1/h1^D [A, B]_D
        # If [A, B]_D is smooth, we fill the matrix accordingly
        for a in F
            for b in F
                A = X^D + a
                B = X^(D-1) + b
                N = bracket(D, A, B, h0, h1)
                O = divexact(N, systFactor)

                tmpFact = factor(O)
                if isSmooth(tmpFact, 2)

                    j += 1
                    fillMatrixGKZ2!(M, j, A, B, f1, D, tmpFact,
                                   unknownsToCoord)

                end
            end
        end

        # Then we erase the blank columns and we transpose the matrix to later perform
        # a Gaussian elimination
        M = subMatrix(M, nbukn+1, j)
        M = transpose(M)

        # We start by finding the logarithms of the elements in the factor
        # base modulo the small factors, using Pohlig Hellman algorithm
        n = ZZ()

        for k in keys(unknownsToCoord)
            dlogs[k] = pohligHellman(card, g, k, defPol)[1]
        end

        dlogs[systFactor], n = pohligHellman(card, g, systFactor, defPol)

        # One we have dealt with the small factors, we compute the big ones
        v = div(card-1, n)
        bigFactors = factor(v)

        # And we compute the logarithms modulo each factor, concluding
        # using the Chinese Remaindering Theorem. With each factor, we obtain
        # a one-dimensional kernel so we have to scale by the generator.
        for p in bigFactors
            Mred = rrefMod(M, p[1])[1]
            ker = Nemo.fmpz[Mred[i, nbukn+1] for i in 1:nbukn]
            s = gcdinv(ker[unknownsToCoord[g]], p[1])[2]
            for j in 1:nbukn
                dlogs[coordToUnknowns[j]] =
                crt(dlogs[coordToUnknowns[j]], n, ker[j]*s, p[1])
            end
            dlogs[systFact] = crt(dlogs[systFact], n, -s, p[1])
            n *= p[1]
        end

        return dlogs
    end
end


text:
using Soss
using Test

# write your own tests here
@testset "Soss" begin
    2 == 2
end


text:
register(DataDep(
    "pudl-raw-eia861",
    """
    Dataset: PUDL Raw EIA Form 861
    Website: https://doi.org/10.5281/zenodo.4127029
    Author: Gosnell, Christina, Selvans, Zane
    Date of Publication: October 24, 2020
    License:

    Raw US Energy Information Administration (EIA) Form 861 Data, archived from https://www.eia.gov/electricity/data/eia861/

    EIA Form 861 Annual Electric Power Industry Report, detailed data files.

    $PUDL_DESCRIPTION
    """,
    [
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/datapackage.json",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2001.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2002.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2003.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2004.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2005.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2006.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2007.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2008.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2009.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2010.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2011.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2012.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2013.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2014.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2015.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2016.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2017.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2018.zip",
        "https://zenodo.org/api/files/7670bb0f-ef46-4b0c-96ce-5b728a985e09/eia861-2019.zip",
    ],
    "2007a19aae76c22f6875c39ea3d8b8113602203556a1b0e2de18cb9f96fc17be",
))