text:
using Poly
using Modeler

## Build problem objects
V1 = Variable("VOLT_1", Complex)
V2 = Variable("VOLT_2", Complex)

obj= 192 * V1*conj(V1) - (96.2+481im) * conj(V2)*V1 + (-96.2+481im) * conj(V1)*V2

Load2 = (-96.2 -481im) * conj(V1)*V2 + ( 96.2 +481im) * conj(V2)*V2
Load2_lb = -350 + 350im; Load2_ub = -350 + 350im
Unit1 = ( 96.2 +481im) * conj(V1)*V1 + (-96.2 -481im) * conj(V2)*V1
Unit1_lb = -400im; Unit1_ub = 600 + 400im
Voltm1 = abs2(V1)
Voltm1_lb = 0.90; Voltm1_ub = 1.10
Voltm2 = abs2(V2)
Voltm2_lb = 0.90; Voltm2_ub = 1.05

## Assemble problem
pb = Problem()
add_variable!(pb, V1)
add_variable!(pb, V2)

set_objective!(pb, obj)

add_constraint!(pb, "UNIT_1", Unit1_lb << Unit1 << Unit1_ub)
add_constraint!(pb, "LOAD_2", Load2_lb << Load2 << Load2_ub)
add_constraint!(pb, "VOLTM_1", Voltm1_lb << Voltm1 << Voltm1_ub)
add_constraint!(pb, "VOLTM_2", Voltm2_lb << Voltm2 << Voltm2_ub)

println(pb)

## Applying at a local solution
pt = Point([V1, V2], [-0.635304 + 0.706321im, 0.380961 + 0.898656im])

println("⟶ Constraints at Knitro local OPF solution:")
for (cstrName, cstr) in get_constraints(pb)
  println(cstrName, ": ", cstr.lb, " ≤ ", evaluate(cstr.p, pt), " ≤ ", cstr.ub)
end

println("⟶ Objective value: ", evaluate(pb.objective, pt))
println("⟶ Corresponding slacks: ")
print(get_slacks(pb, pt))
println("... and min slack: ", get_minslack(pb, pt))

## In real variables:
pb_R = pb_cplx2real(pb)
pt_R = Point([Variable("VOLT_1_r", Real),  Variable("VOLT_1_i", Real),  Variable("VOLT_2_r", Real),  Variable("VOLT_2_i", Real)], [-0.635304, 0.706321, 0.380961, 0.898656])

println("\n⟶ Constraints at Knitro local OPF solution, real problem:")
for (cstrName, cstr) in get_constraints(pb_R)
  println(cstrName, ": ", cstr.lb, " ≤ ", evaluate(cstr.p, pt_R), " ≤ ", cstr.ub)
end

println("⟶ Objective value, real problem: ", evaluate(pb_R.objective, pt_R))

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Yao
using YaoExtensions, YaoTensorNetwork, BitBasis
using YaoBlocks.Optimise
using Fire

# the switch for seperating control gates
sepctrl = false
seperatecontrol(sepctrl)
rk = 4-sepctrl

@main function qft(nbit::Int)
    config = bit_literal(ones(Int,nbit)...)   # final output config
    c = qft_circuit(nbit)
    gb = circuit2tn(c, final_config=config)
    dump_graph(joinpath(@__DIR__, "graphs", "qft$nbit-rank$rk"), gb)
end

@main function variational(nbit::Int, depth::Int)
    c = variational_circuit(nbit, depth, pair_square(isqrt(nbit), isqrt(nbit); periodic=false))
    gb = circuit2tn(c)  # collapse to 00000 state
    dump_graph(joinpath(@__DIR__, "graphs", "variational$nbit-d$depth-rank$rk"), gb)
end

@main function variational2d(nx::Int, ny::Int, depth::Int)
    nbit = nx*ny
    c = variational_circuit(nbit, depth, pair_square(nx, ny; periodic=false))
    gb = circuit2tn(c)  # collapse to 00000 state
    dump_graph(joinpath(@__DIR__, "graphs", "variational2d$nbit-d$depth-rank$rk"), gb)
end

# dump supremacy circuit
@main function supremacy2d(nx::Int, ny::Int, depth::Int, seed::Int=1)
    nbit = nx*ny
    c = rand_supremacy2d(nx, ny, depth; seed=seed)
    gb = circuit2tn(simplify_blocktypes(c))  # collapse to 00000 state
    dump_graph(joinpath(@__DIR__, "graphs", "supremacy2d$nbit-d$depth-rank$rk"), gb)
end

# dump supremacy circuit
@main function google53(depth::Int, nbits::Int=53, seed::Int=1)
    c = rand_google53(depth; seed=seed, nbits=nbits)
    fs = FSimGate(0.5π, π/6)
    if sepctrl
        c = replace_block(fs => openbox(fs), c)
    end
    gb = circuit2tn(simplify_blocktypes(c))  # collapse to 00000 state
    dump_graph(joinpath(@__DIR__, "graphs", "google53$nbits-d$depth-rank$rk"), gb)
end

ext: jl


text:
using PyPlot, PyCall
include("amplt.jl")
PyDict(matplotlib["rcParams"])["axes.formatter.use_locale"] = true
sexframes = Dict(1=>readtable(joinpath("data", "mse14.csv")),
	2=>readtable(joinpath("data", "fse14.csv")))
sexlabels = ["män"; "kvinnor"]
cdicts = vcat(map(sex -> map(ncd -> amplt(ratechange(sexframes[sex], 4/5, 1, ncd), sex, "mort") |>
	f->changes(f, amplt(sexframes[sex], sex, "mort")) |>
	c -> Dict(:sex => sex, :ncd => ncd, :cframe => c),
	[1;4/5;3/5]), [2;1])...)
map(cd -> plot(cd[:cframe][:age], cd[:cframe][:relchange],
	label = "$(sexlabels[cd[:sex]]), ncd=$(cd[:ncd])", "-o"), cdicts)
xlabel("Ålder")
ylabel("RR")
grid(1)
title("Relativ risk upp till ålder för cd=0,8, Sverige 2014")
legend()

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
#=
The number, 197, is called a circular prime because all rotations of the digits: 197, 971, and 719, are themselves prime.

There are thirteen such primes below 100: 2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, and 97.

How many circular primes are there below one million?
=#
include("projecteulerutils.jl")
using Primes

function iscircularprime(x)
  i = reverse(digits(x))
  a = copy(i)
  circularprime = true
  for j in i
    if !isprime(nd(a))
      return false
    end
    a = prepend!(a, [pop!(a)])
  end
  true
end

function calc()
  count = 0
  for x in primes(2_000_000)
    if iscircularprime(x)
      count += 1
    end
  end
  count
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using JuMP
using Ipopt
using Random
# for plots
using Gadfly
using DataFrames

Random.seed!(1234);

x = rand(10)
a = rand(10)
b = rand()

projection = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(projection, x0[1:10])
@objective(projection, Min, sum((x - x0) .* (x - x0))) # We minimize the square of the distance here
@constraint(projection, x0' * a == b)                  # Point must lie on the hyperplane

optimize!(projection)
@show objective_value(projection);
@show value.(x0);

A_1 = rand(10, 10)
A_2 = rand(10, 10)
b_1 = rand(10)
b_2 = rand(10)

polyhedra_distance = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(polyhedra_distance, x[1:10])                       # Point closest on the first polyhedron
@variable(polyhedra_distance, y[1:10])                       # Point closest on the second polyhedron
@objective(polyhedra_distance, Min, sum((x - y) .* (x - y))) # We minimize the square of the distance here as above
@constraint(polyhedra_distance, A_1 * x .<= b_1)             # Point x must lie on the first polyhedron
@constraint(polyhedra_distance, A_2 * y .<= b_2)             # Point y must lie on the second polyhedron

optimize!(polyhedra_distance)
@show objective_value(polyhedra_distance);

fixed = [ 1   1  -1 -1    1   -1  -0.2  0.1;         # coordinates of fixed points
          1  -1  -1  1 -0.5 -0.2    -1    1]

M = size(fixed,2)                                    # number of fixed points
N = 6                                                # number of free points

A = [ 1  0  0 -1  0  0    0  0  0  0  0  0  0  0;    # Matrix on links
      1  0 -1  0  0  0    0  0  0  0  0  0  0  0;
      1  0  0  0 -1  0    0  0  0  0  0  0  0  0;
      1  0  0  0  0  0   -1  0  0  0  0  0  0  0;
      1  0  0  0  0  0    0 -1  0  0  0  0  0  0;
      1  0  0  0  0  0    0  0  0  0 -1  0  0  0;
      1  0  0  0  0  0    0  0  0  0  0  0  0 -1;
      0  1 -1  0  0  0    0  0  0  0  0  0  0  0;
      0  1  0 -1  0  0    0  0  0  0  0  0  0  0;
      0  1  0  0  0 -1    0  0  0  0  0  0  0  0;
      0  1  0  0  0  0    0 -1  0  0  0  0  0  0;
      0  1  0  0  0  0    0  0 -1  0  0  0  0  0;
      0  1  0  0  0  0    0  0  0  0  0  0 -1  0;
      0  0  1 -1  0  0    0  0  0  0  0  0  0  0;
      0  0  1  0  0  0    0 -1  0  0  0  0  0  0;
      0  0  1  0  0  0    0  0  0  0 -1  0  0  0;
      0  0  0  1 -1  0    0  0  0  0  0  0  0  0;
      0  0  0  1  0  0    0  0 -1  0  0  0  0  0;
      0  0  0  1  0  0    0  0  0 -1  0  0  0  0;
      0  0  0  1  0  0    0  0  0  0  0 -1  0  0;
      0  0  0  1  0  0    0  0  0  0  0 -1  0  0;
      0  0  0  0  1 -1    0  0  0  0  0  0  0  0;
      0  0  0  0  1  0   -1  0  0  0  0  0  0  0;
      0  0  0  0  1  0    0  0  0 -1  0  0  0  0;
      0  0  0  0  1  0    0  0  0  0  0  0  0 -1;
      0  0  0  0  0  1    0  0 -1  0  0  0  0  0;
      0  0  0  0  0  1    0  0  0  0 -1  0  0  0;]

placement = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(placement, p[1:M + N, 1:2])                     # A variable array for the coordinates of each point
@constraint(placement, p[N + 1:N + M, :] .== fixed')      # We had a constraint for the fixed points
dist = A * p                                         # Matrix of differences between coordinates of 2 points with a link
@objective(placement, Min, sum(dist .* dist))            # We minimize the sum of the square of the distances

optimize!(placement)
@show value.(p);
@show objective_value(placement);

# Plotting the points
df = DataFrame()
df.x = value.(p)[:,1]
df.y = value.(p)[:,2]
df.type = vcat(fill("Free points", N), fill("Fixed points", M))
plt = plot(df, x = "x", y = "y", color = "type", Geom.point)
draw(SVG(6inch, 6inch), plt)

n = 5;

Amin = [                                        # We'll try this problem with 4 times with different minimum area constraints
100 100 100 100 100;
 20  50  80 150 200;
180  80  80  80  80;
 20 150  20 200 110]

r = 1

figs=[]

for i = 1:4
    local A = Amin[i, :]

    floor_planning = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    local x, y, w, h, W, H
    @variables(floor_planning, begin
        x[1:n] >= r
        y[1:n] >= r
        w[1:n] >= 0
        h[1:n] >= 0
        W
        H
    end)

    @constraints(floor_planning, begin
        x[5] + w[5] + r <= W    # No rectangles at the right of Rectangle 5
        x[1] + w[1] + r <= x[3] # Rectangle 1 is at the left of Rectangle 3
        x[2] + w[2] + r <= x[3] # Rectangle 2 is at the left of Rectangle 3
        x[3] + w[3] + r <= x[5] # Rectangle 3 is at the left of Rectangle 5
        x[4] + w[4] + r <= x[5] # Rectangle 4 is at the left of Rectangle 5
        y[4] + h[4] + r <= H    # No rectangles on top of Rectangle 4
        y[5] + h[5] + r <= H    # No rectangles on top of Rectangle 5
        y[2] + h[2] + r <= y[1] # Rectangle 2 is below Rectangle 1
        y[1] + h[1] + r <= y[4] # Rectangle 1 is below Rectangle 4
        y[3] + h[3] + r <= y[4] # Rectangle 3 is below Rectangle 4
        w .<= 5*h               # Aspect ratio constraint
        h .<= 5*w               # Aspect ratio constraint
        A .<= h .*  w           # Area constraint
    end)

    @objective(floor_planning, Min, W + H)

    optimize!(floor_planning)

    @show objective_value(floor_planning);

    D = DataFrame(x = value.(x), y = value.(y), x2 = value.(x) .+ value.(w), y2 = value.(y) .+ value.(h))
    local plt = plot(D, xmin = :x, ymin = :y, xmax = :x2, ymax = :y2, Geom.rect)
    push!(figs, plt)
end

draw(SVG(6inch, 6inch), vstack(hstack(figs[1], figs[2]), hstack(figs[3], figs[4])))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VaultSecretGroup <: SwaggerModel
    sourceVault::Any # spec type: Union{ Nothing, SubResource } # spec name: sourceVault
    vaultCertificates::Any # spec type: Union{ Nothing, Vector{VaultCertificate} } # spec name: vaultCertificates

    function VaultSecretGroup(;sourceVault=nothing, vaultCertificates=nothing)
        o = new()
        validate_property(VaultSecretGroup, Symbol("sourceVault"), sourceVault)
        setfield!(o, Symbol("sourceVault"), sourceVault)
        validate_property(VaultSecretGroup, Symbol("vaultCertificates"), vaultCertificates)
        setfield!(o, Symbol("vaultCertificates"), vaultCertificates)
        o
    end
end # type VaultSecretGroup

const _property_map_VaultSecretGroup = Dict{Symbol,Symbol}(Symbol("sourceVault")=>Symbol("sourceVault"), Symbol("vaultCertificates")=>Symbol("vaultCertificates"))
const _property_types_VaultSecretGroup = Dict{Symbol,String}(Symbol("sourceVault")=>"SubResource", Symbol("vaultCertificates")=>"Vector{VaultCertificate}")
Base.propertynames(::Type{ VaultSecretGroup }) = collect(keys(_property_map_VaultSecretGroup))
Swagger.property_type(::Type{ VaultSecretGroup }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VaultSecretGroup[name]))}
Swagger.field_name(::Type{ VaultSecretGroup }, property_name::Symbol) =  _property_map_VaultSecretGroup[property_name]

function check_required(o::VaultSecretGroup)
    true
end

function validate_property(::Type{ VaultSecretGroup }, name::Symbol, val)
end


text:
"""
This module provides concrete implementations of `Number` that represent
1st, 2nd and general order tensors.

## Why

The main feature of this module is that the provided types do not extend from `AbstractArray`, but from `Number`!

This allows one to work with them as if they were scalar values in broadcasted operations on arrays of `VectorValue` objects (also for `TensorValue` or `MultiValue` objects). For instance, one can perform the following manipulations:
```julia
# Assing a VectorValue to all the entries of an Array of VectorValues
A = zeros(VectorValue{2,Int}, (4,5))
v = VectorValue(12,31)
A .= v # This is posible since  VectorValue <: Number

# Broadcasting of tensor operations in arrays of TensorValues
t = TensorValue(13,41,53,17) # creates a 2x2 TensorValue
g = TensorValue(32,41,3,14) # creates another 2x2 TensorValue
B = fill(t,(1,5))
C = inner.(g,B) # inner product of g against all TensorValues in the array B
@show C
# C = [2494 2494 2494 2494 2494]
```

The exported names are:

$(EXPORTS)
"""
module TensorValues

using DocStringExtensions

using StaticArrays
using StaticArrays: SVector, MVector, SMatrix, MMatrix, SArray, MArray
using Base: @propagate_inbounds, @pure
using Gridap.Helpers
using Gridap.Arrays
using Gridap: CellField
using LinearAlgebra
using Random
using TensorOperations

export MultiValue
export VectorValue
export TensorValue
export SymTensorValue
export SymFourthOrderTensorValue
export ThirdOrderTensorValue

export inner, outer, meas
export mutable
export Mutable
export symmetric_part
export n_components
export num_components
export change_eltype
export diagonal_tensor
export ⊙
export ⊗
export ⋅¹
export ⋅²
export double_contraction
export contract
export @GTensor

import Base: show
import Base: zero, one
import Base: +, -, *, /, \, ==, ≈, isless
import Base: conj
import Base: sum, maximum, minimum
import Base: getindex, iterate, eachindex
import Base: size, length, eltype
import Base: reinterpret
import Base: convert
import Base: CartesianIndices
import Base: LinearIndices
import Base: adjoint
import Base: transpose
import Base: rand

import LinearAlgebra: det, inv, tr, cross, dot, norm
# Reexport from LinearAlgebra (just for convenience)
export det, inv, tr, cross, dot, norm, ×, ⋅

import Gridap.Arrays: get_array

include("MultiValueTypes.jl")

include("VectorValueTypes.jl")

include("TensorValueTypes.jl")

include("SymTensorValueTypes.jl")

include("SymFourthOrderTensorValueTypes.jl")

include("ThirdOrderTensorValueTypes.jl")

include("Indexing.jl")

include("Operations.jl")

include("Reinterpret.jl")

end # module


text:
# Plot the result in space (across all x) for a specific angular frequency
@recipe function plot(simres::FrequencySimulationResult, ω::AbstractFloat;
        phase_time = 0.0, # does not except "t" as name of variable..
        x_indices = axes(simres.x,1),
        ω_index = findmin(abs.(getfield(simres, 3) .- ω))[2],
        field_apply = real, seriestype = :surface)

    x = [x[1] for x in simres.x[x_indices]]

    # y will actually be z for 3D...
    y = [x[end] for x in simres.x[x_indices]]
    ω = getfield(simres, 3)[ω_index]

    phase = exp(-im*ω*phase_time)

    seriescolor --> :balance
    title --> "Field for ω = $ω"
    seriestype --> seriestype
    aspect_ratio --> 1.0

    if seriestype == :contour
        # We could check here to see if x and y have the right structure
        x = unique(x)
        y = unique(y)

        n_x = length(x)
        n_y = length(y)

        fill --> true
        x, y, field_apply.(phase.*transpose(reshape(field(simres)[x_indices,ω_index],n_x,n_y)))

    else
        (x, y, field_apply.(phase.*field(simres)[x_indices,ω_index]))
    end

end

# Plot the result in space (across all x) for a specific angular frequency
@recipe function plot(timres::TimeSimulationResult, t::AbstractFloat;
        x_indices = axes(timres.x,1),
        t_index = findmin(abs.(getfield(timres, 3) .- t))[2],
        field_apply = real, seriestype = :surface)

    x = [x[1] for x in timres.x[x_indices]]
    y = [x[2] for x in timres.x[x_indices]]
    t = getfield(timres, 3)[t_index]

    seriescolor --> :balance
    title --> "Field for time = $t"
    seriestype --> seriestype
    aspect_ratio --> 1.0

    if seriestype == :contour

        # We should really check here to see if x and y have the right structure
        x = unique(x)
        y = unique(y)

        n_x = length(x)
        n_y = length(y)

        fill --> true
        x, y, field_apply.(transpose(reshape(field(timres)[x_indices,t_index],n_y,n_x)))

    else

        (x, y, field_apply.(field(timres)[x_indices,t_index]))

    end

end


text:
function WindUpFactory()
    return WindUpFactory(())
end

function get_wind_up(obj::WindUpFactory, arg0::SatelliteSystem, arg1::jint, arg2::JString)
    return jcall(obj, "getWindUp", WindUp, (SatelliteSystem, jint, JString), arg0, arg1, arg2)
end


text:
################################################################################
#                          EVALUATION METHOD DATATYPES
################################################################################
"""
    AbstractIntegralMethod

An abstract type for integral evaluation methods use in combination with
`integral` and `generate_integral_data`.
"""
abstract type AbstractIntegralMethod end

"""
    Automatic <: AbstractIntegralMethod

An integral evaluation type for automically selecting an appropriate integral
evaluation method. Contains no fields.
"""
struct Automatic <: AbstractIntegralMethod end

"""
    AbstractUnivariateMethod <: AbstractIntegralMethod

An abstract type for integral evaluation methods for 1-dimensional integrals.
"""
abstract type AbstractUnivariateMethod <: AbstractIntegralMethod end

"""
    UniTrapezoid <: AbstractUnivariateMethod

An integral evalution method that uses the trapezoid rule to in combination
with all parameter supports available when the integral is expanded and/or when
the infinite model is optimized, whichever comes first. Note this method will
ignore the `num_supports` keyword argument. Note this is valid only for finite
integral domains. Contains no fields.
"""
struct UniTrapezoid <: AbstractUnivariateMethod end

"""
    UniMCSampling <: AbstractUnivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral. This variant will add more supports to the model as
needed to satisfy `num_supports` and it will include all supports with the
`MCSample` label up till the integral is expanded and/or when
the infinite model is optimized, whichever comes first. Note this is valid only
for finite integral domains. Contains no fields.
"""
struct UniMCSampling <: AbstractUnivariateMethod end

"""
    UniIndepMCSampling <: AbstractUnivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral similar to [`UniMCSampling`](@ref MeasureToolbox.UniMCSampling).
However, this variant will generate its own set of supports and ignore all other
supports with the `MCSample` label. Note this is valid only for finite integral
domains. This is not compatible with individual dependent parameters. Contains no fields.
"""
struct UniIndepMCSampling <: AbstractUnivariateMethod end

"""
    Quadrature <: AbstractUnivariateMethod

A general integral evaluation method that will automatically select the
appropriate quadrature method to approximate the integral. Please note that this
will generate a unique set of parameter supports and will ignore existing supports
when the integral is evaluated and thus should be used with caution. However,
this method is able to handle infinite and semi-infinite integral domains.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct Quadrature <: AbstractUnivariateMethod end

"""
    GaussHermite <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Hermite quadrature to
evaluate integrals. This is valid for infinite integral domains. Note this will
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussHermite <: AbstractUnivariateMethod end

"""
    GaussLegendre <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Legendre quadrature to
evaluate integrals. This is valid for finite integral domains. Note this will
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussLegendre <: AbstractUnivariateMethod end

"""
    GaussRadau <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Radau quadrature to
evaluate integrals. This is valid for finite integral domains. Note this will
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussRadau <: AbstractUnivariateMethod end

"""
    GaussLobatto <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Lobatto quadrature to
evaluate integrals. This is valid for finite integral domains. Note this will
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussLobatto <: AbstractUnivariateMethod end

"""
    GaussJacobi <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Jacobi quadrature to
evaluate integrals. This is valid for finite integral domains. This requires the user
to input the alpha and beta shape parameters for their function. This will then
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters.

**Fields**
- `α::Float64`: The α Jacobi shape parameter.
- `β::Float64`: The β Jacobi shape parameter.
"""
struct GaussJacobi <: AbstractUnivariateMethod
    α::Float64
    β::Float64
    function GaussJacobi(α::Real, β::Real)
        return new(Float64(α), Float64(β))
    end
end

"""
    GaussLaguerre <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Laguerre quadrature to
evaluate integrals. This is valid for semi-infinite integral domains. Note this
will generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussLaguerre <: AbstractUnivariateMethod end

"""
    AbstractMultivariateMethod <: AbstractIntegralMethod

An abstract type for integral evaluation methods for multi-dimensional integrals.
"""
abstract type AbstractMultivariateMethod <: AbstractIntegralMethod end

"""
    MultiMCSampling <: AbstractMultivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral. This variant will add more supports to the model as
needed to satisfy `num_supports` and it will include all supports with the
`MCSample` label up till the integral is expanded and/or when
the infinite model is optimized, whichever comes first. Note this is valid only
for finite integral domains. If an array of independent infinite parameters
is specified, they must use the same amount of supports. Contains no fields.
"""
struct MultiMCSampling <: AbstractMultivariateMethod end

"""
    MultiIndepMCSampling <: AbstractMultivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral similar to [`MultiMCSampling`](@ref MeasureToolbox.MultiMCSampling).
However, this variant will generate its own set of supports and ignore all other
supports with the `MCSample` label. Note this is valid only for finite integral
domains. Contains no fields.
"""
struct MultiIndepMCSampling <: AbstractMultivariateMethod end

################################################################################
#                             INTERNAL SUPPORT LABELS
################################################################################
"""
    InternalGaussLobatto <: InfiniteOpt.InternalLabel

A support label Gauss Lobatto points that are used as generative supports.
"""
struct InternalGaussLobatto <: InfiniteOpt.InternalLabel end

################################################################################
#                       DATA GENERATION METHODS (UNIVARIATE)
################################################################################
"""
    generate_integral_data(
        prefs::Union{InfiniteOpt.GeneralVariableRef, Vector{InfiniteOpt.GeneralVariableRef}},
        lower_bounds::Union{Real, Vector{<:Real}},
        upper_bounds::Union{Real, Vector{<:Real}},
        method::V; [num_supports::Int = InfiniteOpt.DefaultNumSupports,
        weight_func::Function = InfiniteOpt.default_weight,
        extra_kwargs...]
        )::InfiniteOpt.AbstractMeasureData where {V <: AbstractIntegralMethod}

Generate the appropriate concrete realization of `AbstractMeasureData` using
`method`. Here `prefs`, `lower_bounds`, and `upper_bounds` will always have a
1 to 1 correspondence when this is called from `integral`. Please refer
to the method docstrings for an explanation of each one.

User-defined method extensions should first
define a concrete `method` type inheriting from `AbstractUnivariateMethod` or
`AbstractMultivariateMethod` as appropriate and then implement extend this
method using that type for `method`.
"""
function generate_integral_data end

# General fallback
function generate_integral_data(prefs, lb, ub, method; kwargs...)
    error("`generate_integral_data` is not defined for method type `$(method)` " *
          "in accordance with the arguments given.")
end

# Single pref with Automatic
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
    lower_bound::Real,
    upper_bound::Real,
    method::Automatic;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    is_depend = InfiniteOpt._index_type(pref) == DependentParameterIndex
    inf_bound_num = (lower_bound == -Inf) + (upper_bound == Inf)
    if inf_bound_num == 0 # finite interval
        method = UniTrapezoid()
    elseif inf_bound_num == 1 && !is_depend # semi-infinite interval
        method = GaussLaguerre()
    elseif inf_bound_num == 2 && !is_depend # infinite interval
        method = GaussHermite()
    else
        error("Cannot generate measure data for individual dependent parameters " *
              "with infinite or semi-infinite domains.")
    end
    return generate_integral_data(pref, lower_bound, upper_bound, method,
                                  num_supports = num_supports,
                                  weight_func = weight_func)
end

# Univariate trapezoid coefficient function
function _trapezoid_coeff(supps::Vector{<:Real})::Vector{Float64}
    len = length(supps)
    len >= 2 || error("Cannot invoke trapezoid rule on integral if there are " *
                      "less than 2 supports added to its infinite parameter. Ensure " *
                      "all infinite parameters have supports.")
    coeffs = Vector{Float64}(undef, len)
    for i in eachindex(supps)
        if i == 1
            coeffs[i] = 1 / 2 * (supps[2] - supps[1])
        elseif i == len
            coeffs[i] = 1 / 2 * (supps[len] - supps[len-1])
        else
            coeffs[i] = 1 / 2 * (supps[i+1] - supps[i-1])
        end
    end
    return coeffs
end

# Single pref trapezoid rule
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::UniTrapezoid;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    return InfiniteOpt.FunctionalDiscreteMeasureData(pref, _trapezoid_coeff, 0,
                                                     InfiniteOpt.All,
                                                     InfiniteOpt.NoGenerativeSupports(),
                                                     weight_func, lower_bound,
                                                     upper_bound, false)
end

# Useful error function for dependent parameters
function _ensure_independent_param(pref::InfiniteOpt.GeneralVariableRef,
                                   method)::Nothing
    if InfiniteOpt._index_type(pref) == InfiniteOpt.DependentParameterIndex
        error("Cannot generate measure data for individual dependent parameters " *
              "via `$(method)`. Try using `UniTrapezoid` for finite domains.")
    end
    return
end

# Single pref general quadrature dispatch
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::Quadrature;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    inf_bound_num = (lower_bound == -Inf) + (upper_bound == Inf)
    if inf_bound_num == 0 # finite interval
        method = GaussLegendre()
    elseif inf_bound_num == 1 # semi-infinite interval
        method = GaussLaguerre()
    else # infinite interval
        method = GaussHermite()
    end
    return generate_integral_data(pref, lower_bound, upper_bound, method,
                                  num_supports = num_supports, weight_func = weight_func)
end

# Single pref Gauss-Legendre
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussLegendre;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    if lower_bound == -Inf || upper_bound == Inf
        @warn("Gauss Legendre quadrature can only be applied on finite intervals, " *
              "switching to an appropriate method.")
        return generate_integral_data(pref, lower_bound, upper_bound, Quadrature(),
                                      num_supports = num_supports,
                                      weight_func = weight_func)
    end
    (supports, coeffs) = FastGaussQuadrature.gausslegendre(num_supports)
    supports = (upper_bound - lower_bound) / 2 * supports .+ (upper_bound + lower_bound) / 2
    coeffs = (upper_bound - lower_bound) / 2 * coeffs
    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.generate_unique_label(),
                                           weight_func, lower_bound, upper_bound,
                                           false)
end

# Single pref Gauss-Radau
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussRadau;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    if lower_bound == -Inf || upper_bound == Inf
        @warn("Gauss Radau quadrature can only be applied on finite intervals, " *
            "switching to an appropriate method.")
    return generate_integral_data(pref, lower_bound, upper_bound, Quadrature(),
                                        num_supports = num_supports,
                                        weight_func = weight_func)
    end
    (supports, coeffs) = FastGaussQuadrature.gaussradau(num_supports)
    supports = (upper_bound - lower_bound) / 2 * supports .+ (upper_bound + lower_bound) / 2
    coeffs = (upper_bound - lower_bound) / 2 * coeffs
    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                            InfiniteOpt.generate_unique_label(),
                                            weight_func, lower_bound, upper_bound,
                                            false)
end

# Single pref Gauss-Lobatto
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussLobatto;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    if lower_bound == -Inf || upper_bound == Inf
        @warn("Gauss Lobatto quadrature can only be applied on finite intervals, " *
            "switching to an appropriate method.")
    return generate_integral_data(pref, lower_bound, upper_bound, Quadrature(),
                                        num_supports = num_supports,
                                        weight_func = weight_func)
    end
    (supports, coeffs) = FastGaussQuadrature.gausslobatto(num_supports)
    supports = (upper_bound - lower_bound) / 2 * supports .+ (upper_bound + lower_bound) / 2
    coeffs = (upper_bound - lower_bound) / 2 * coeffs
return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                            InfiniteOpt.generate_unique_label(),
                                            weight_func, lower_bound, upper_bound,
                                            false)
end

# Single pref Gauss-Jacobi
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussJacobi;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    if lower_bound == -Inf || upper_bound == Inf
        @warn("Gauss Jacobi quadrature can only be applied on finite intervals, " *
            "switching to an appropriate method.")
    return generate_integral_data(pref, lower_bound, upper_bound, Quadrature(),
                                        num_supports = num_supports,
                                        weight_func = weight_func)
    end
    (supports, coeffs) = FastGaussQuadrature.gaussjacobi(num_supports, method.α, method.β)
    supports = (upper_bound - lower_bound) / 2 * supports .+ (upper_bound + lower_bound) / 2
    coeffs = (upper_bound - lower_bound) / 2 * coeffs
return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                            InfiniteOpt.generate_unique_label(),
                                            weight_func, lower_bound, upper_bound,
                                            false)
end

# Single pref Gauss-Laguerre
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussLaguerre;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    if upper_bound == Inf
        if lower_bound == -Inf
            @warn("Gauss Laguerre quadrature can only be applied on semi-infinite intervals, " *
                  "switching to an appropriate method.")
            return generate_integral_data(pref, lower_bound, upper_bound, GaussHermite(),
                                          num_supports = num_supports,
                                          weight_func = weight_func)
        end
        (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_supports)
        coeffs = coeffs .* exp.(supports)
        supports = supports .+ lower_bound
    elseif lower_bound == -Inf
        (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_supports)
        coeffs = coeffs .* exp.(supports)
        supports = -supports .+ upper_bound
    else
        @warn("Gauss Laguerre quadrature can only be applied on semi-infinite intervals, " *
              "switching to an appropriate method.")
        return generate_integral_data(pref, lower_bound, upper_bound, GaussLegendre(),
                                      num_supports = num_supports,
                                      weight_func = weight_func)
    end

    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.generate_unique_label(),
                                           weight_func, lower_bound, upper_bound,
                                           false)
end

# Single pref Gauss-Hermite
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussHermite;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    if lower_bound != -Inf || upper_bound != Inf
        @warn("Gauss Hermite quadrature can only be applied on infinite intervals, " *
              "switching to an appropriate method.")
        return generate_integral_data(pref, lower_bound, upper_bound, Quadrature(),
                                      num_supports = num_supports,
                                      weight_func = weight_func)
    end
    (supports, coeffs) = FastGaussQuadrature.gausshermite(num_supports)
    coeffs = coeffs .* exp.(supports.^2)
    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.generate_unique_label(),
                                           weight_func, lower_bound, upper_bound,
                                           false)
end

# Single pref uniform Monte Carlo sampling
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::UniMCSampling;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    # check and process the arguments
    if lower_bound == -Inf || upper_bound == Inf
        error("Univariate MC sampling is not supported for (semi-)infinite intervals.")
    elseif InfiniteOpt._index_type(pref) == InfiniteOpt.DependentParameterIndex
        if !(num_supports in [0, InfiniteOpt.DefaultNumSupports])
            @warn("Cannot specify a nonzero minimum supports for an individual " *
                  "dependent parameter. Setting `num_supports = 0`.")
        end
        num_supports = 0
    end
    # make the coefficient function
    function _coeffs(supps::Vector{<:Real})::Vector{Float64}
        len = length(supps)
        return (upper_bound - lower_bound) * ones(len) / len
    end
    # prepare the data
    return InfiniteOpt.FunctionalDiscreteMeasureData(pref, _coeffs, num_supports,
                                                     InfiniteOpt.MCSample,
                                                     InfiniteOpt.NoGenerativeSupports(),
                                                     weight_func,
                                                     lower_bound, upper_bound,
                                                     false)
end

# Single pref independent Monte Carlo sampling
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::UniIndepMCSampling;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight)
    _ensure_independent_param(pref, method)
    if lower_bound == -Inf || upper_bound == Inf
        error("Univariate MC sampling is not applicable to (semi-)infinite intervals.")
    end
    set = InfiniteOpt.IntervalSet(lower_bound, upper_bound)
    supports, _ = InfiniteOpt.generate_support_values(set, InfiniteOpt.MCSample,
                                                      num_supports = num_supports)
    coeffs = (upper_bound - lower_bound) * ones(num_supports) / num_supports
    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.MCSample, weight_func,
                                           lower_bound, upper_bound, false)
end

################################################################################
#                       DATA GENERATION METHODS (MULTIVARIATE)
################################################################################
# Vector of prefs with Automatic
function generate_integral_data(prefs::Vector{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Vector{<:Real},
    upper_bounds::Vector{<:Real},
    method::Automatic;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    return generate_integral_data(prefs, lower_bounds, upper_bounds,
                                  MultiMCSampling(), num_supports = num_supports,
                                  weight_func = weight_func)
end

# Vector of prefs using Monte Carlo samples
function generate_integral_data(prefs::Vector{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Vector{<:Real},
    upper_bounds::Vector{<:Real},
    method::MultiMCSampling;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    # check bounds
    if any(lb == -Inf for lb in lower_bounds) || any(ub == Inf for ub in upper_bounds)
        error("MC Sampling is not applicable to (semi-)infinite intervals.")
    end
    # make the coefficients function
    vol = prod(upper_bounds .- lower_bounds)
    function _coeffs(supps::Array{<:Real, 2})::Vector{Float64}
        num_samples = size(supps, 2)
        return vol * ones(num_samples) / num_samples
    end
    # make the data
    return InfiniteOpt.FunctionalDiscreteMeasureData(prefs, _coeffs,
                                                     num_supports,
                                                     InfiniteOpt.MCSample,
                                                     weight_func,
                                                     lower_bounds, upper_bounds,
                                                     false)
end

## Generate multivariate Monte Carlo samples in accordance with the parameter type
# DependentParameterRefs
function _make_multi_mc_supports(prefs::Vector{InfiniteOpt.DependentParameterRef},
                                 lbs::Vector{<:Real}, ubs::Vector{<:Real},
                                 num_supps::Int)::Matrix{Float64}
    sets = [InfiniteOpt.IntervalSet(lbs[i], ubs[i]) for i in eachindex(lbs)]
    set = InfiniteOpt.CollectionSet(sets)
    supports, _ = InfiniteOpt.generate_support_values(set, InfiniteOpt.MCSample,
                                                      num_supports = num_supps)
    return supports
end

# Vector of prefs with independent Monte Carlo supports
function generate_integral_data(prefs::Vector{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Vector{<:Real},
    upper_bounds::Vector{<:Real},
    method::MultiIndepMCSampling;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    # check bounds
    if any(lb == -Inf for lb in lower_bounds) || any(ub == Inf for ub in upper_bounds)
        error("MC Sampling is not applicable to (semi-)infinite intervals.")
    end
    # prepare the supports
    sets = [InfiniteOpt.IntervalSet(lower_bounds[i], upper_bounds[i])
            for i in eachindex(lower_bounds)]
    set = InfiniteOpt.CollectionSet(sets)
    supports, _ = InfiniteOpt.generate_support_values(set, InfiniteOpt.MCSample,
                                                      num_supports = num_supports)
    # prepare the coefficients
    coeffs = prod(upper_bounds .- lower_bounds) * ones(num_supports) / num_supports
    # prepare the data
    return InfiniteOpt.DiscreteMeasureData(prefs, coeffs, supports,
                                           InfiniteOpt.MCSample, weight_func,
                                           lower_bounds, upper_bounds, false)
end

################################################################################
#                           UNIVARIATE INTEGRALS
################################################################################
# Define default keyword arguments for 1-D integrals
const UniIntegralDefaults = Dict(:eval_method => Automatic(),
                                 :num_supports => InfiniteOpt.DefaultNumSupports,
                                 :weight_func => InfiniteOpt.default_weight)

"""
    uni_integral_defaults()::Dict{Symbol, Any}

Get the default keyword argument values for defining one-dimensional integrals.

```julia-repl
julia> uni_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :num_supports          => 10
  :eval_method           => Automatic()
  :weight_func           => default_weight
```
"""
uni_integral_defaults()::Dict{Symbol, Any} = UniIntegralDefaults

"""
    set_uni_integral_defaults(; kwargs...)::Nothing

Set the default keyword argument settings for one-dimensional integrals.
The keyword arguments of this function will be recorded in the default keyword
argument dictionary. These will determine the default keyword argument values
when calling [`integral`](@ref MeasureToolbox.integral(::JuMP.AbstractJuMPScalar,::InfiniteOpt.GeneralVariableRef,::Real, ::Real))
with a single infinite parameter.

**Example**
```julia-repl
julia> uni_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :num_supports          => 10
  :eval_method           => Automatic()
  :weight_func           => default_weight

julia> set_uni_integral_defaults(num_supports = 5, eval_method = Quadrature(),
                                 new_kwarg = true)

julia> uni_integral_defaults()
Dict{Symbol,Any} with 4 entries:
  :new_kwarg             => true
  :num_supports          => 5
  :eval_method           => Quadrature()
  :weight_func           => default_weight
```
"""
function set_uni_integral_defaults(; kwargs...)::Nothing
    merge!(UniIntegralDefaults, kwargs)
    return
end

"""
    integral(expr::JuMP.AbstractJuMPScalar,
             pref::GeneralVariableRef,
             [lower_bound::Real = lower_bound(pref),
             upper_bound::Real = upper_bound(pref);
             kwargs...])::GeneralVariableRef

Returns a measure reference that evaluates the integral of `expr` with respect
to infinite parameter `pref` from `lower_bound` to `upper_bound`. This thus
considers integrals of the form: ``\\int_{p \\in P} expr(p) w(p) dp`` where ``p``
is an infinite parameter and ``w`` is the weight function is 1 by default. This
function provides a high-level interface that ultimately constructs an appropriate
concrete form of [`AbstractMeasureData`](@ref) via [`generate_integral_data`](@ref)
in accordance with the keyword arugment `eval_method` that is then used with
[`measure`](@ref). Note that it is preferred to call [`@integral`](@ref) when
`expr` is not just a single variable reference. Errors for bad bound input.

The keyword arguments are as follows:
- `eval_method::AbstractUnivariateMethod`: Used to determine the
    numerical evaluation scheme. Possible choices include:
    - [`Automatic`](@ref MeasureToolbox.Automatic)
    - [`UniTrapezoid`](@ref MeasureToolbox.UniTrapezoid)
    - [`UniMCSampling`](@ref MeasureToolbox.UniMCSampling)
    - [`UniIndepMCSampling`](@ref MeasureToolbox.UniIndepMCSampling)
    - [`Quadrature`](@ref MeasureToolbox.Quadrature)
    - [`GaussHermite`](@ref MeasureToolbox.GaussHermite)
    - [`GaussLegendre`](@ref MeasureToolbox.GaussLegendre)
    - [`GaussLageurre`](@ref MeasureToolbox.GaussLaguerre)
    - [`GaussLobatto`](@ref MeasureToolbox.GaussLobatto)
    - [`GaussRadau`](@ref MeasureToolbox.GaussRadau)
    - [`GaussJacobi`](@ref MeasureToolbox.GaussJacobi)
- `num_supports`: The minimum number of supports to be generated (if used by
    `eval_method`)
- `weight_func`: ``w(p)`` above with parameter value inputs and scalar output

See [`set_uni_integral_defaults`](@ref) to update the default keyword argument
values for all one-dimensional integral calls.

**Example**
```julia-repl
julia> @infinite_parameter(model, x in [0, 1])
x

julia> @infinite_variable(model, f(x))
f(x)

julia> int = integral(f, x)
∫{x ∈ [0, 1]}[f(x)]

julia> expand(int)
0.2 f(0.8236475079774124) + 0.2 f(0.9103565379264364) + 0.2 f(0.16456579813368521) + 0.2 f(0.17732884646626457) + 0.2 f(0.278880109331201)
```
"""
function integral(expr::JuMP.AbstractJuMPScalar,
                  pref::InfiniteOpt.GeneralVariableRef,
                  lower_bound::Real = NaN,
                  upper_bound::Real = NaN;
                  kwargs...)::InfiniteOpt.GeneralVariableRef
    # check parameter formatting
    InfiniteOpt._check_params(pref)
    # fill in bounds if needed
    set = infinite_set(pref)
    if isnan(lower_bound) && JuMP.has_lower_bound(set)
        lower_bound = JuMP.lower_bound(pref)
    end
    if isnan(upper_bound) && JuMP.has_upper_bound(set)
        upper_bound = JuMP.upper_bound(pref)
    end
    # ensure valid bounds
    if lower_bound >= upper_bound
        error("Invalid integral bounds, ensure that lower_bound < upper_bound.")
    elseif !InfiniteOpt.supports_in_set([lower_bound, upper_bound], set)
        error("Integral bounds violate the infinite domain.")
    end
    # prepare the keyword arguments and make the measure data
    processed_kwargs = merge(uni_integral_defaults(), kwargs)
    eval_method = pop!(processed_kwargs, :eval_method)
    data = generate_integral_data(pref, lower_bound, upper_bound,
                                  eval_method; processed_kwargs...)
    # make the measure
    return InfiniteOpt.measure(expr, data, name = "integral")
end

"""
    ∫(expr::JuMP.AbstractJuMPScalar,
      pref::GeneralVariableRef,
      [lower_bound::Real = NaN,
      upper_bound::Real = NaN;
      kwargs...])::GeneralVariableRef

A convenient wrapper for [`integral`](@ref). The `∫` unicode symbol is produced
via `\\int`.
"""
function ∫(expr::JuMP.AbstractJuMPScalar,
           pref::InfiniteOpt.GeneralVariableRef,
           lower_bound::Real = NaN,
           upper_bound::Real = NaN;
           kwargs...)::InfiniteOpt.GeneralVariableRef
    return integral(expr, pref, lower_bound, upper_bound; kwargs...)
end

################################################################################
#                            MULTIVARIATE INTEGRALS
################################################################################
# Define default keyword arguments for multi-D integrals
const MultiIntegralDefaults = Dict(:eval_method => Automatic(),
                                   :num_supports => InfiniteOpt.DefaultNumSupports,
                                   :weight_func => InfiniteOpt.default_weight)

"""
    multi_integral_defaults()::Dict{Symbol, Any}

Get the default keyword argument values for defining multi-dimensional integrals.

```julia-repl
julia> multi_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :num_supports          => 10
  :eval_method           => Automatic()
  :weight_func           => default_weight
```
"""
multi_integral_defaults()::Dict{Symbol, Any} = MultiIntegralDefaults

"""
    set_multi_integral_defaults(; kwargs...)::Nothing

Set the default keyword argument settings for multi-dimesnional integrals.
The keyword arguments of this function will be recorded in the default keyword
argument dictionary. These will determine the default keyword argument values
when calling [`integral`](@ref MeasureToolbox.integral(::JuMP.AbstractJuMPScalar,::AbstractArray{GeneralVariableRef},::Union{Real, AbstractArray{<:Real}}, ::Union{Real, AbstractArray{<:Real}}))
with an array of infinite parameters.

**Example**
```julia-repl
julia> multi_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :num_supports          => 10
  :eval_method           => Automatic()
  :weight_func           => default_weight

julia> set_multi_integral_defaults(num_supports = 5, new_kwarg = true)

julia> multi_integral_defaults()
Dict{Symbol,Any} with 4 entries:
  :new_kwarg             => true
  :num_supports          => 5
  :eval_method           => Automatic()
  :weight_func           => default_weight
```
"""
function set_multi_integral_defaults(; kwargs...)::Nothing
    merge!(MultiIntegralDefaults, kwargs)
    return
end

Base.isnan(arr::AbstractArray{<:Real})::Bool = all(Base.isnan.(arr))
"""
    integral(expr::JuMP.AbstractJuMPScalar,
             prefs::AbstractArray{GeneralVariableRef},
             [lower_bounds::Union{Real, AbstractArray{<:Real}} = [lower_bound(pref)...],
             upper_bounds::Union{Real, AbstractArray{<:Real}} = [upper_bound(pref)...];
             kwargs...])::GeneralVariableRef

Returns a measure reference that evaluates the integral of `expr` with respect
to infinite parameters `prefs` from `lower_bounds` to `upper_bounds`. This thus
considers integrals of the form: ``\\int_{p \\in P} expr(p) w(p) dp`` where ``p``
is an infinite parameter and ``w`` is the weight function is 1 by default. This
function provides a high-level interface that ultimately constructs an appropriate
concrete form of [`AbstractMeasureData`](@ref) via [`generate_integral_data`](@ref)
in accordance with the keyword arugment `eval_method` that is then used with
[`measure`](@ref). Note that it is preferred to call [`@integral`](@ref) when
`expr` is not just a single variable reference. Errors when the container types
and dimensions do not match or the bounds are invalid.

The keyword arguments are as follows:
- `eval_method::AbstractMultivariateMethod`: Used to determine the
    numerical evaluation scheme. Possible choices include:
    - [`Automatic`](@ref MeasureToolbox.Automatic)
    - [`MultiMCSampling`](@ref MeasureToolbox.MultiMCSampling)
    - [`MultiIndepMCSampling`](@ref MeasureToolbox.MultiIndepMCSampling)
- `num_supports`: The minimum number of supports to be generated (if used by
    `eval_method`)
- `weight_func`: ``w(p)`` above with parameter value inputs and scalar output

See [`set_multi_integral_defaults`](@ref) to update the default keyword argument
values for all multi-dimensional integral calls.

**Example**
```julia-repl
julia> @infinite_parameter(model, x[1:2] in [0, 1], independent = true);

julia> @infinite_variable(model, f(x));

julia> int = integral(f, x)
∫{x ∈ [0, 1]^2}[f(x)]
```
"""
function integral(expr::JuMP.AbstractJuMPScalar,
                  prefs::AbstractArray{InfiniteOpt.GeneralVariableRef},
                  lower_bounds::Union{Real, AbstractArray{<:Real}} = NaN,
                  upper_bounds::Union{Real, AbstractArray{<:Real}} = NaN;
                  kwargs...)::InfiniteOpt.GeneralVariableRef
    # check parameter formatting
    InfiniteOpt._check_params(prefs)
    # fill in the lower bounds if needed
    if isnan(lower_bounds) && JuMP.has_lower_bound(first(prefs))
        lbs = map(p -> JuMP.lower_bound(p), prefs)
    elseif !isnan(lower_bounds) && lower_bounds isa Real
        lbs = map(p -> lower_bounds, prefs)
    else
        lbs = lower_bounds
    end
    # fill in the upper bounds if needed
    if isnan(upper_bounds) && JuMP.has_upper_bound(first(prefs))
        ubs = map(p -> JuMP.upper_bound(p), prefs)
    elseif !isnan(upper_bounds) && upper_bounds isa Real
        ubs = map(p -> upper_bounds, prefs)
    else
        ubs = upper_bounds
    end
    # do initial bound checks
    if InfiniteOpt._keys(prefs) != InfiniteOpt._keys(lbs) ||
       InfiniteOpt._keys(prefs) != InfiniteOpt._keys(ubs)
        error("Array keys of infinite parameters and bounds do not match.")
    end
    # process the arrays
    vector_prefs = InfiniteOpt._make_ordered_vector(prefs)
    vector_lbs = InfiniteOpt._make_ordered_vector(lbs)
    vector_ubs = InfiniteOpt._make_ordered_vector(ubs)
    # ensure valid bounds
    if any(vector_lbs[i] >= vector_ubs[i] for i in eachindex(vector_lbs))
        error("Invalid integral bounds, ensure that lower_bounds < upper_bounds.")
    end
    dprefs = map(p -> dispatch_variable_ref(p), vector_prefs)
    InfiniteOpt._check_bounds_in_set(dprefs, vector_lbs, vector_ubs)
    # prepare the keyword arguments and make the measure data
    processed_kwargs = merge(multi_integral_defaults(), kwargs)
    eval_method = pop!(processed_kwargs, :eval_method)
    data = generate_integral_data(vector_prefs, vector_lbs, vector_ubs,
                                  eval_method; processed_kwargs...)
    # make the measure
    return InfiniteOpt.measure(expr, data, name = "integral")
end

"""
    ∫(expr::JuMP.AbstractJuMPScalar,
      prefs::AbstractArray{GeneralVariableRef},
      [lower_bounds::Union{Real, AbstractArray{<:Real}} = NaN,
      upper_bounds::Union{Real, AbstractArray{<:Real}} = NaN;
      kwargs...])::GeneralVariableRef

A convenient wrapper for [`integral`](@ref). The unicode symbol `∫` is produced
via `\\int`.
"""
function ∫(expr::JuMP.AbstractJuMPScalar,
           prefs::AbstractArray{InfiniteOpt.GeneralVariableRef},
           lower_bounds::Union{Real, AbstractArray{<:Real}} = NaN,
           upper_bounds::Union{Real, AbstractArray{<:Real}} = NaN;
           kwargs...)::InfiniteOpt.GeneralVariableRef
    return integral(expr, prefs, lower_bounds, upper_bounds; kwargs...)
end

"""
    @integral(expr::JuMP.AbstractJuMPScalar,
              prefs::Union{GeneralVariableRef, AbstractArray{GeneralVariableRef}},
              [lower_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds,
              upper_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds;
              kwargs...])::GeneralVariableRef

An efficient wrapper for [`integral`](@ref integral(::JuMP.AbstractJuMPScalar, ::InfiniteOpt.GeneralVariableRef, ::Real, ::Real))
and [`integral`](@ref integral(::JuMP.AbstractJuMPScalar, ::AbstractArray{InfiniteOpt.GeneralVariableRef}, ::Union{Real, AbstractArray{<:Real}}, ::Union{Real, AbstractArray{<:Real}})).
Please see the above doc strings for more information.
"""
macro integral(expr, prefs, args...)
    _error(str...) = InfiniteOpt._macro_error(:integral, (expr, prefs, args...),
                                              str...)
    extra, kw_args, requestedcontainer = InfiniteOpt._extract_kw_args(args)
    if length(extra) != 0 && length(extra) != 2
        _error("Incorrect number of positional arguments for @integral. " *
               "Must provide both bounds or no bounds.")
    end
    expression = :( JuMP.@expression(InfiniteOpt._Model, $expr) )
    mref = :( integral($expression, $prefs, $(extra...); ($(kw_args...))) )
    return esc(mref)
end

"""
    @∫(expr::JuMP.AbstractJuMPScalar,
       prefs::Union{GeneralVariableRef, AbstractArray{GeneralVariableRef}},
       [lower_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds,
       upper_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds;
       kwargs...])::GeneralVariableRef

A convenient wrapper for [`@integral`](@ref). The unicode symbol `∫` is produced
via `\\int`.
"""
macro ∫(expr, prefs, args...)
    return esc(:( @integral($expr, $prefs, $(args...)) ))
end

################################################################################
#                               HELPER METHODS
################################################################################
#=
"""
    infinite_transform(lb::Number, ub::Number, num_supports::Int;
                       [sub_method::Function = mc_sampling,
                       transform_x::Function = _default_x,
                       transform_dx::Function = _default_dx,
                       t_lb::Number = -convert(Number, lb == -Inf && ub == Inf),
                       t_ub::Number = 1., kwargs...])::Tuple

Returns a tuple that contains supports and coefficients generated for a
parameter in an infinite or semi-infinite interval. It works by transforming
the original unbounded interval to a finite interval, on which a support
generation method for finite intervals is applied. Then, the generated supports
are transformed back to the original interval. The user is allowed to specify
the support generation method for finite intevals to use, as well as the
transform function. The default transform function is
``t \\in [-\\infty, \\infty] \\rightarrow x \\in [-1, 1]: t(x) = \\frac{t}{1-t^2}``
``t \\in [a, \\infty] \\rightarrow x \\in [0, 1]: t(x) = a + \\frac{t}{1-t}``
``t \\in [-\\infty, a] \\rightarrow x \\in [0, 1]: t(x) = a - \\frac{1-t}{t}``

**Example**
```jldoctest; setup = :(using InfiniteOpt)
julia> (supps, coeffs) = infinite_transform(-Inf, Inf, 5, sub_method = gauss_legendre())
([-5.06704059565454, -0.7583532171678754, 0.0, 0.7583532171678754, 5.06704059565454], [13.490960583398396, 1.2245949721571516, 0.5688888888888889, 1.2245949721571516, 13.490960583398396])
```
"""
function infinite_transform(set::InfiniteOpt.IntervalSet,
                            params::Union{InfiniteOpt.ParameterRef,
                            AbstractArray{<:InfiniteOpt.ParameterRef}},
                            num_supports::Int,
                            lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                            ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                            sub_method::Val = Val(mc_sampling);
                            transform_x::Function = _default_x,
                            transform_dx::Function = _default_dx,
                            t_lb::Number = -convert(Number, lb == -Inf && ub == Inf),
                            t_ub::Number = 1.)::Tuple
    # transform (semi-)infinite domain to finite domain
    if lb != -Inf && ub != Inf
        error("The range is not (semi-)infinite. Use evaluation methods for " *
              "bounded domains.")
    end
    (t_supports, t_coeffs) = generate_integral_data(set, params, num_supports, t_lb, t_ub, sub_method)
    supports = transform_x.(t_supports, lb, ub)
    coeffs = t_coeffs .* transform_dx.(t_supports, lb, ub)
    return (supports, coeffs)
end

function _default_x(t::Number, lb::Number, ub::Number)::Number
    if lb > -Inf
        return lb + t / (1 - t)
    elseif ub < Inf
        return ub - (1 - t) / t
    else
        return t / (1 - t^2)
    end
end

function _default_dx(t::Number, lb::Number, ub::Number)::Number
    if lb > -Inf
        return 1 / (1 - t)^2
    elseif ub < Inf
        return 1 / t^2
    else
        return (1 + t^2) / (1 - t^2)^2
    end
end
=#


text:
using Distributions
import Random
import DrWatson
import TimeSeries
using PowerSystems
import CSV
import Dates
Random.seed!(123)

const PSY = PowerSystems
include(joinpath(DrWatson.srcdir(),"data/data_download.jl"))
include(joinpath(DrWatson.srcdir(),"data/make_scenarios.jl"))
include(joinpath(DrWatson.srcdir(),"data/utility_functions.jl"))

sha = "367916992a15d5fccf6ce23220c4abae1c27de3d"
base_dir = DrWatson.datadir()

RTS_ROOT = obtain_raw_data(base_dir, sha)
TimeSeries_DIR = joinpath(RTS_ROOT,"RTS_Data/timeseries_data_files")

include(joinpath(base_dir, "5_bus_base_data.jl"))

################################# Generate UC System Data ##################################
uc_system = System(100.0)
for b in nodes5
    add_component!(uc_system, b)
end

for b in branches5
    add_component!(uc_system, b)
end

for b in thermal_generators5_uc_testing
    add_component!(uc_system, b)
end

for b in loads5
    add_component!(uc_system, b)
end

add_component!(uc_system, il)

for b in renewable_generators5
    add_component!(uc_system, b)
end

add_component!(uc_system, reserve5)

#Attach Time Series data to Unit Commitment
reserves1 = CSV.read(joinpath(TimeSeries_DIR,"Reserves/DAY_AHEAD_regional_Spin_Up_R1.csv"))
reserves2 = CSV.read(joinpath(TimeSeries_DIR,"Reserves/DAY_AHEAD_regional_Spin_Up_R2.csv"))
reserves3 = CSV.read(joinpath(TimeSeries_DIR,"Reserves/DAY_AHEAD_regional_Spin_Up_R3.csv"))
total_reserves = 0.01*(reserves1[!,:Spin_Up_R1] + reserves2[!,:Spin_Up_R2] + reserves3[!,:Spin_Up_R3]);
Time_stamps = Dates.DateTime.(reserves1[!,:Year],reserves1[!,:Month], reserves1[!,:Day], (reserves1[!,:Period] .- 1))
reserves_time_array = TimeSeries.TimeArray(Time_stamps, total_reserves, [:Spinning]);

PowerSystems.add_forecast!(uc_system, reserves_time_array, reserve5, "Spinning", 1.0)

wind = CSV.read(joinpath(TimeSeries_DIR,"WIND/DAY_AHEAD_wind.csv"))
wind[!,Symbol("309_WIND_1")] = wind[!,Symbol("309_WIND_1")]./maximum(wind[!,Symbol("309_WIND_1")])
wind[!,Symbol("317_WIND_1")] = wind[!,Symbol("317_WIND_1")]./maximum(wind[!,Symbol("317_WIND_1")])
wind[!,Symbol("303_WIND_1")] = wind[!,Symbol("303_WIND_1")]./maximum(wind[!,Symbol("303_WIND_1")])
wind[!,Symbol("122_WIND_1")] = wind[!,Symbol("122_WIND_1")]./maximum(wind[!,Symbol("122_WIND_1")])
Time_stamps = Dates.DateTime.(wind[!,:Year],wind[!,:Month], wind[!,:Day], wind[!,:Period] .- 1)

wind_forecast1 = PowerSystems.Deterministic(renewable_generators5[1],"Deterministic",TimeSeries.TimeArray(Time_stamps, wind[!,Symbol("309_WIND_1")]));
wind_forecast2 = PowerSystems.Deterministic(renewable_generators5[2],"Deterministic",TimeSeries.TimeArray(Time_stamps, wind[!,Symbol("317_WIND_1")]));
wind_forecast3 = PowerSystems.Deterministic(renewable_generators5[3],"Deterministic",TimeSeries.TimeArray(Time_stamps, wind[!,Symbol("122_WIND_1")]));
wind_forecasts = [wind_forecast1, wind_forecast2, wind_forecast3]
PowerSystems.add_forecasts!(uc_system, wind_forecasts)

load = CSV.read(joinpath(TimeSeries_DIR,"Load/DAY_AHEAD_regional_Load.csv"));
load[!,Symbol("1")] = load[!,Symbol("1")]./maximum(load[!,Symbol("1")])
load[!,Symbol("2")] = load[!,Symbol("2")]./maximum(load[!,Symbol("2")])
load[!,Symbol("3")] = load[!,Symbol("3")]./maximum(load[!,Symbol("3")])
Time_stamps = Dates.DateTime.(load[!,:Year],load[!,:Month], load[!,:Day], load[!,:Period].- 1)

load_forecast1 = PowerSystems.Deterministic(loads5[1],"Deterministic",TimeSeries.TimeArray(Time_stamps, load[!,Symbol("1")]));
load_forecast2 = PowerSystems.Deterministic(loads5[2],"Deterministic",TimeSeries.TimeArray(Time_stamps, load[!,Symbol("2")]));
load_forecast3 = PowerSystems.Deterministic(loads5[3],"Deterministic",TimeSeries.TimeArray(Time_stamps, load[!,Symbol("3")]));
load_forecast_interruptible = PowerSystems.Deterministic(il, "Deterministic",TimeSeries.TimeArray(Time_stamps, load[!,Symbol("3")]));
load_forecasts = [load_forecast1, load_forecast2, load_forecast3, load_forecast_interruptible]
PowerSystems.add_forecasts!(uc_system, load_forecasts)

PowerSystems.split_forecasts!(uc_system,
                PowerSystems.get_forecasts(Deterministic, uc_system, Dates.Dates.DateTime("2020-01-01T00:00:00")),
                Dates.Hour(24),
                48)

# Write System Data To File
to_json(uc_system, joinpath(base_dir,"uc_system.json"))

################################# Generate SUC System Data ##################################

scenario_forecasts = make_scenarios(uc_system, 100, 0.2)

suc_system = System(100.0)
for b in nodes5
    add_component!(suc_system, b)
end

for b in branches5
    add_component!(suc_system, b)
end

for b in thermal_generators5_uc_testing
    add_component!(suc_system, b)
end

for b in loads5
    add_component!(suc_system, b)
end

add_component!(suc_system, il)

for b in renewable_generators5
    add_component!(suc_system, b)
end

add_component!(suc_system, reserve5)

PowerSystems.add_forecasts!(suc_system, load_forecasts)
PowerSystems.add_forecast!(suc_system, reserves_time_array, reserve5, "Spinning", 1.0)

PowerSystems.split_forecasts!(suc_system,
                PowerSystems.get_forecasts(Deterministic, suc_system, Dates.Dates.DateTime("2020-01-01T00:00:00")),
                Dates.Hour(24),
                48)

PowerSystems.add_forecasts!(suc_system, scenario_forecasts)

PowerSystems.to_json(suc_system, joinpath(base_dir,"suc_system.json"))

################################# Generate ED System Data ##################################

ed_system = System(100.0)
for b in nodes5
    add_component!(ed_system, b)
end

for b in branches5
    add_component!(ed_system, b)
end

for b in thermal_generators5_uc_testing
    add_component!(ed_system, b)
end

for b in loads5
    add_component!(ed_system, b)
end

add_component!(ed_system, il)

for b in renewable_generators5
    add_component!(ed_system, b)
end

add_component!(ed_system, reserve5)

reserves1 = CSV.read(joinpath(TimeSeries_DIR, "Reserves/REAL_TIME_regional_Spin_Up_R1.csv"))
reserves2 = CSV.read(joinpath(TimeSeries_DIR, "Reserves/REAL_TIME_regional_Spin_Up_R2.csv"))
reserves3 = CSV.read(joinpath(TimeSeries_DIR, "Reserves/REAL_TIME_regional_Spin_Up_R3.csv"))
total_reserves = 0.01*(reserves1[!,:Spin_Up_R1] + reserves2[!,:Spin_Up_R2] + reserves3[!,:Spin_Up_R3]);

Time_stamps = Dates.DateTime.(reserves1[!,:Year],reserves1[!,:Month], reserves1[!,:Day]) .+ Dates.Minute(5).*(reserves1[!,:Period])
Time_stamps = Time_stamps .- Dates.Minute(5)

reserves_time_array = TimeSeries.TimeArray(Time_stamps, total_reserves, [:Spinning]);
add_forecast!(ed_system, reserves_time_array, reserve5, "Spinning", 1.0)

wind = CSV.read(joinpath(TimeSeries_DIR,"WIND/REAL_TIME_wind.csv"))
wind[!,Symbol("309_WIND_1")] = wind[!,Symbol("309_WIND_1")]./maximum(wind[!,Symbol("309_WIND_1")])
wind[!,Symbol("317_WIND_1")] = wind[!,Symbol("317_WIND_1")]./maximum(wind[!,Symbol("317_WIND_1")])
wind[!,Symbol("303_WIND_1")] = wind[!,Symbol("303_WIND_1")]./maximum(wind[!,Symbol("303_WIND_1")])
wind[!,Symbol("122_WIND_1")] = wind[!,Symbol("122_WIND_1")]./maximum(wind[!,Symbol("122_WIND_1")])

wind_forecast1 = Deterministic(renewable_generators5[1],"Deterministic",TimeSeries.TimeArray(Time_stamps, wind[!,Symbol("309_WIND_1")]));
wind_forecast2 = Deterministic(renewable_generators5[2],"Deterministic",TimeSeries.TimeArray(Time_stamps, wind[!,Symbol("317_WIND_1")]));
wind_forecast3 = Deterministic(renewable_generators5[3],"Deterministic",TimeSeries.TimeArray(Time_stamps, wind[!,Symbol("122_WIND_1")]));
wind_forecasts = [wind_forecast1, wind_forecast2, wind_forecast3]
add_forecasts!(ed_system, wind_forecasts)

load = CSV.read(joinpath(TimeSeries_DIR,"Load/REAL_TIME_regional_Load.csv"));
load[!,Symbol("1")] = load[!,Symbol("1")]./maximum(load[!,Symbol("1")])
load[!,Symbol("2")] = load[!,Symbol("2")]./maximum(load[!,Symbol("2")])
load[!,Symbol("3")] = load[!,Symbol("3")]./maximum(load[!,Symbol("3")])

load_forecast1 = Deterministic(loads5[1],"Deterministic", TimeSeries.TimeArray(Time_stamps, load[!,Symbol("1")]));
load_forecast2 = Deterministic(loads5[2],"Deterministic", TimeSeries.TimeArray(Time_stamps, load[!,Symbol("2")]));
load_forecast3 = Deterministic(loads5[3],"Deterministic", TimeSeries.TimeArray(Time_stamps, load[!,Symbol("3")]));
load_forecast_interruptible = Deterministic(il, "Deterministic",TimeSeries.TimeArray(Time_stamps, load[!,Symbol("3")]));
load_forecasts = [load_forecast1, load_forecast2, load_forecast3, load_forecast_interruptible]
add_forecasts!(ed_system, load_forecasts)

split_forecasts!(ed_system,
                get_forecasts(Deterministic, ed_system, Dates.Dates.DateTime("2020-01-01T00:00:00")),
                Dates.Hour(1),
                12)

# Write the ED System to JSON
to_json(ed_system, joinpath(base_dir,"ed_system.json"))


text:
using Random, LinearAlgebra
using Stheno: WoodburyMat, Xt_invA_X, LazyPDMat

@testset "woodbury" begin

    let
        rng, N, D = MersenneTwister(123456), 10, 5
        X, Σ_, σ = randn(rng, D, N), randn(rng, D, D), exp(randn(rng))
        Σ = LazyPDMat(Σ_ * Σ_' + 1e-9I, 0)
        A = WoodburyMat(X, Σ, σ)
        B, b = randn(N, 2), randn(N)

        # Test that conversion and indexing are consistent with one another.
        Am, Al = Matrix(A), LazyPDMat(A)
        @test size(A) == (N, N)
        @test Matrix(Al) == Am
        @test Am[2, 3] ≈ A[2, 3]
        @test Am[3, 2] ≈ A[3, 2]

        # Check unary operations.
        @test logdet(A) ≈ logdet(Matrix(A))

        # Test basics of matrix multiplication.
        @test Matrix(A * B) ≈ Matrix(A) * B
        @test Vector(A * b) ≈ Matrix(A) * b
        Bt = B'
        @test Matrix(Bt * A) ≈ B' * Matrix(A)

        # Test basics of linear systems solving.
        @test Matrix(A \ B) ≈ Matrix(A) \ B
        @test Vector(A \ b) ≈ Matrix(A) \ b

        @test Matrix(B' / A) ≈ B' / Matrix(A)

        # Test quadratic form solving.
        @test Xt_invA_X(A, B) isa Matrix{<:Real}
        @test Xt_invA_X(A, B) ≈ B' * (Matrix(A) \ B)

        @test Xt_invA_X(A, b) isa Real
        @test Xt_invA_X(A, b) ≈ b' * (Matrix(A) \ b)
    end

end


text:
# do the work

function a_readable_number(str::String, rns::ReadableNumStyle)
    local integral_part, fractional_part, integral_readable, fractional_readable, readable
    integral_part, fractional_part = split(str, FRACPOINT)[1:2]
    integral_readable = ifelse( integral_part == "", "0",
                    readable_integer(integral_part, rns.integral_digits_spanned, rns.between_integral_spans) )
    fractional_readable = ifelse( fractional_part == "", "" ,
                    readable_fraction(fractional_part, rns.fractional_digits_spanned, rns.between_fractional_spans) )
    readable = (fractional_part == "") ? integral_readable : string(integral_readable, rns.between_parts, fractional_readable)
    return readable
end

function readable_nonneg_integer(str::String, digits_spanned::I, group_separator::Char) where {I <: Integer}
    n = length(str)
    n==0 && return "0"

    sinteger, sexponent =
        if occursin("e", str)
           split(str,'e')
        else
           str, ""
        end

    n = length(sinteger)

    fullgroups, finalgroup = divrem(n, digits_spanned)

    sv = Vector{Char}(sinteger)
    p = repeat(" ", n+(fullgroups-1)+(finalgroup!=0))
    pretty = Vector{Char}(p)

    sourceidx = n
    targetidx = length(pretty)
    for k in fullgroups:-1:1
        pretty[(targetidx-digits_spanned+1):targetidx] = sv[(sourceidx-digits_spanned+1):sourceidx]
        sourceidx -= digits_spanned
        targetidx -= digits_spanned
        if k > 1
            pretty[targetidx] = group_separator
            targetidx -= 1
        end
    end

    if finalgroup > 0
        if fullgroups > 0
            pretty[targetidx] = group_separator
            targetidx -= 1
        end
        pretty[(targetidx-finalgroup+1):targetidx] = sv[(sourceidx-finalgroup+1):sourceidx]
    end

    prettystring = String(pretty)

    if length(sexponent) != 0
       string(prettystring,"e",sexponent)
    else
       prettystring
    end
end

function readable_integer(str::String, digits_spanned::I, group_separator::Char) where {I <: Integer}
    if str[1] != "-"
       readable_nonneg_integer(str, digits_spanned, group_separator)
    else
       s1 = string(s[2:end])
       integral_readable = readable_nonneg_integer(s1, digits_spanned, group_separator)
       string("-", integral_readable)
    end
end

function readable_fraction(str::String, digits_spanned::I, group_separator::Char) where {I <: Integer}
    sfrac, sexponent =
        if occursin("e", str)
           split(str,'e')
        else
           str, ""
        end

    fractional_readable = reverse(readable_nonneg_integer(reverse(sfrac), digits_spanned, group_separator))

    if length(sexponent) != 0
       string(fractional_readable,"e",sexponent)
    else
       fractional_readable
    end
end


text:
using JLD
using Gen
using Parameters: @with_kw

include("../trace-translator/hmm.jl")

abstract type Discretization end

@with_kw mutable struct CustomDiscretization <: Discretization
    a::Float64
    b::Float64
    c::Float64
    d::Float64
    e::Float64
    f::Float64
    g::Float64
    h::Float64
    i::Float64
    j::Float64
    k::Float64
    l::Float64
    bounds::Bounds
    bin_bounds::Vector{Bounds}
end

function set_bin_bounds!(disc::CustomDiscretization)
    bin_bounds = Vector{Bounds}(undef, get_num_bins(disc))

    # above a
    bin_bounds[1] = Bounds(0.0, disc.d, disc.a, 1.0)
    bin_bounds[2] = Bounds(disc.d, disc.f, disc.a, 1.0)
    bin_bounds[3] = Bounds(disc.f, 1.0, disc.a, 1.0)

    # above b
    bin_bounds[4] = Bounds(0.0, disc.d, disc.b, disc.a)
    bin_bounds[5] = Bounds(disc.d, disc.f, disc.b, disc.a)
    bin_bounds[6] = Bounds(disc.f, disc.g, disc.b, disc.a)
    bin_bounds[7] = Bounds(disc.g, 1.0, disc.b, disc.a)

    # above e, left of g
    bin_bounds[8] = Bounds(0.0, disc.d, disc.e, disc.b)
    bin_bounds[9] = Bounds(disc.d, disc.h, disc.e, disc.b)
    bin_bounds[10] = Bounds(disc.h, disc.g, disc.e, disc.b)

    # above c, left of g
    bin_bounds[11] = Bounds(0.0, disc.d, disc.c, disc.e)
    bin_bounds[12] = Bounds(disc.d, disc.h, disc.c, disc.e)
    bin_bounds[13] = Bounds(disc.h, disc.g, disc.c, disc.e)

    # above i, right of g
    bin_bounds[14] = Bounds(disc.g, 1.0, disc.i, disc.b)

    # below i, right of g
    bin_bounds[15] = Bounds(disc.g, disc.k, disc.c, disc.i)
    bin_bounds[16] = Bounds(disc.k, 1.0, disc.c, disc.i)

    # below c
    bin_bounds[17] = Bounds(0.0, disc.d, 0.0, disc.c)
    bin_bounds[18] = Bounds(disc.d, disc.j, 0.0, disc.c)
    bin_bounds[19] = Bounds(disc.j, disc.l, 0.0, disc.c)
    bin_bounds[20] = Bounds(disc.l, 1.0, 0.0, disc.c)

    disc.bin_bounds = bin_bounds
end

get_num_bins(disc::CustomDiscretization) = 20

function fix(disc::CustomDiscretization, pt::Point)
    return Point(
        max(disc.bounds.xmin + 1e-10, min(disc.bounds.xmax, pt.x) - 1e-10),
        max(disc.bounds.ymin + 1e-10, min(disc.bounds.ymax, pt.y) - 1e-10))
end

function get_bin(disc::CustomDiscretization, pt::Point)
    if pt.y > disc.b
        if pt.y > disc.a
            if pt.x < disc.d
                return 1
            elseif pt.x < disc.f
                return 2
            else
                return 3
            end
        else
            if pt.x < disc.f
                if pt.x < disc.d
                    return 4
                else
                    return 5
                end
            else
                if pt.x < disc.g
                    return 6
                else
                    return 7
                end
            end
        end
    else

        if pt.y < disc.c
            if pt.x < disc.j
                if pt.x < disc.d
                    return 17
                else
                    return 18
                end
            else
                if pt.x < disc.l
                    return 19
                else
                    return 20
                end
            end
        else
            if pt.x < disc.g
#Point(0.20696641618060418, 0.3581357819832387)
                if pt.y < disc.e
                    if pt.x < disc.d
                        return 11
                    else
                        if pt.x < disc.h
                            return 12
                        else
                            return 13
                        end
                    end
                else
                    if pt.x < disc.d
                        return 8
                    else
                        if pt.x < disc.h
                            return 9
                        else
                            return 10
                        end
                    end
                end
            else
                if pt.y < disc.i
                    if pt.x < disc.k
                        return 15
                    else
                        return 16
                    end
                else
                    return 14
                end
            end
        end
    end
    println(pt)
    @assert false
end

function get_bounds(disc::CustomDiscretization, bin_idx::Int)
    bnds = disc.bin_bounds[bin_idx]
    return (bnds.xmin, bnds.xmax, bnds.ymin, bnds.ymax)
end

const custom_disc = CustomDiscretization(
    a = 0.85,
    b = 0.59,
    c = 0.185,
    d = 0.242,
    e = 0.385,
    f = 0.435,
    g = 0.662,
    h = 0.540,
    i = 0.342,
    j = 0.61,
    k = 0.731,
    l = 0.845,
    bounds = Bounds(0.0, 1.0, 0.0, 1.0),
    bin_bounds = Vector{Bounds}())

set_bin_bounds!(custom_disc)

function load_hmms(start::Point)
    d = load("../trace-translator/trained_hmms_new.jld")
    ps = d["prior"]
    ts = d["transition"]
    es = d["emission"]
    num_waypoints = length(ts)
    hmms = Vector{HiddenMarkovModel}()
    start_bin = get_bin(custom_disc, start)
    for i in 1:num_waypoints
        prior_vec = 1e-6 * ones(get_num_bins(custom_disc))
        prior_vec[start_bin] = 1.0
        prior_vec = prior_vec / sum(prior_vec)
        push!(hmms, HiddenMarkovModel(
            prior_vec, ts[i], es[i]))
    end
    return hmms
end

struct CollapsedHMMTrace <: Trace
    obs::Vector{Int}
    score::Float64
end

Gen.get_score(trace::CollapsedHMMTrace) = trace.score
Gen.project(trace::CollapsedHMMTrace, ::EmptySelection) = 0.0
function Gen.get_choices(trace::CollapsedHMMTrace)
    choices = choicemap()
    for t in 1:length(trace.obs)
        choices[t] = trace.obs[t]
    end
    return choices
end
Gen.get_retval(trace::CollapsedHMMTrace) = nothing

struct CollapsedHMM <: GenerativeFunction{Nothing,CollapsedHMMTrace}
end

const lfprobs = Matrix{Float64}(undef, 10, get_num_bins(custom_disc)) # 10 is the number of time steps
const arr = Vector{Float64}(undef, get_num_bins(custom_disc))

function Gen.generate(::CollapsedHMM, args::Tuple, constraints::ChoiceMap)
    hmm, T = args
    obs = Vector{Int}(undef, T)
    for t in 1:T
        obs[t] = constraints[t]
    end
    score = hmm_log_marginal_likelihood(hmm, obs)
    trace = CollapsedHMMTrace(obs, score)
    return (trace, score)
end

const collapsed_hmm = CollapsedHMM()

@gen function discrete_single_waypoint_model(disc::Discretization, start_idx::Int, T::Int)
    waypoint_idx ~ uniform_discrete(1, get_num_bins(disc))
    observations ~ collapsed_hmm(hmms[waypoint_idx], T) # :observations => 1, ... :observations => T
end

function discrete_single_waypoint_posterior(disc::Discretization, start_idx::Int, obs::Vector{Int})
    choices = choicemap()
    for (i, v) in enumerate(obs)
        choices[:observations => i] = v
    end
    probs = Vector{Float64}(undef, get_num_bins(disc))
    for waypoint_idx in 1:get_num_bins(disc)
        choices[:waypoint_idx] = waypoint_idx
        _, weight = generate(discrete_single_waypoint_model, (disc, start_idx, length(obs)), choices)
        probs[waypoint_idx] = weight
    end
    probs = exp.(probs .- logsumexp(probs))
    @assert isapprox(sum(probs), 1.0, atol=1e-2)
    probs = probs / sum(probs)
    return probs
end

function load_hmms(start::Point)
    d = load("../trace-translator/trained_hmms_new.jld")
    ps = d["prior"]
    ts = d["transition"]
    es = d["emission"]
    num_waypoints = length(ts)
    hmms = Vector{HiddenMarkovModel}()
    start_bin = get_bin(custom_disc, start)
    for i in 1:num_waypoints
        prior_vec = 1e-6 * ones(get_num_bins(custom_disc))
        prior_vec[start_bin] = 1.0
        prior_vec = prior_vec / sum(prior_vec)
        push!(hmms, HiddenMarkovModel(
            prior_vec, ts[i], es[i]))
    end
    return hmms
end


text:
# Copied from Query.jl to simplify the creation of anonymous functions

function helper_replace_anon_func_syntax(ex)
	if !(isa(ex, Expr) && ex.head==:->)
		new_symb = gensym()
		new_ex = MacroTools.postwalk(ex) do x
			if isa(x, Symbol) && x==:_
				return new_symb
			else
				return x
			end
		end
		return :($new_symb -> $(new_ex) )
	else
		return ex
	end
end

function replace_selector(s::S, f, sym::Symbol) where {S<:AbstractSelector}
    fields = fieldnames(S)
    new_fields = Tuple(field == sym ? f : getfield(s, field) for field in fields)
	(S <: Selector) ? Selector(new_fields...) : ColumnSelector(new_fields...)
end

macro splitby(s, arg)
    anon_func = helper_replace_anon_func_syntax(arg)
    Expr(:call, :_splitby, esc(s), esc(anon_func))
end

_splitby(s::AbstractSelector, f) = replace_selector(s, f, :splitby)

macro across(s, arg)
    p = (arg == QuoteNode(:all)) ? QuoteNode(:all) :
        helper_replace_anon_func_syntax(arg)
    Expr(:call, :_across, esc(s), esc(p))
end

function _across(s::AbstractSelector, f)
    if f == :all
        s.kw[:acrossall] = true
        s2 = s
    else
        s2 = replace_selector(s, f, :across)
    end
	haskey(s2.kw, :compute_error) && isa(s2.kw[:compute_error], Number) &&
		error("@bootstrap and @across are not supported in the same pipeline")
    s2.kw[:compute_error] = :across
	s2.kw[:xreduce] = get(s2.kw, :xreduce, mean)
    s2.kw[:yreduce] = get(s2.kw, :yreduce, mean)
    s2
end

macro bootstrap(s, arg = 1000)
    Expr(:call, :_bootstrap, esc(s), esc(arg))
end

function _bootstrap(s::AbstractSelector, n)
	haskey(s.kw, :compute_error) && s.kw[:compute_error] == :across &&
		error("@bootstrap and @across are not supported in the same pipeline")
    s.kw[:compute_error] = n
	s
end

functionize(f) = f
functionize(f::Tuple) = x -> map(g->g(x), f)

macro x(s, x, args...)
    anon_func = helper_replace_anon_func_syntax(x)
    Expr(:call, :_x, esc(s), esc(anon_func), (esc(arg) for arg in args)...)
end

function _x(s::AbstractSelector, f, args...)
    s2 = replace_selector(s, f, :x)
    kws = [:axis_type, :nbins]
	if args == () || isa(args[1], Symbol)
	    for (ind, val) in enumerate(args)
	        s2.kw[kws[ind]] = val
	    end
	else
		s2.kw[:xreduce] = functionize(args[1])
	end
    return s2
end

macro xlims(s, arg)
    p = helper_replace_anon_func_syntax(arg)
    Expr(:call, :_xlims, esc(s), esc(p))
end

_xlims(s, p) = _xlims(ProcessedTable(s), p)
function _xlims(s::ProcessedTable, p)
    filter(s) do i
        lower, upper = p(i)
        x = i.x isa Tuple ? i.x[1] : i.x
        lower <= x <= upper
    end
end

macro ylims(s, arg)
    p = helper_replace_anon_func_syntax(arg)
    Expr(:call, :_ylims, esc(s), esc(p))
end

_ylims(s, p) = _ylims(ProcessedTable(s), p)
function _ylims(s::ProcessedTable, p)
    filter(s) do i
        lower, upper = p(i)
        y = i.y isa Tuple ? i.y[1] : i.y
        lower <= y <= upper
    end
end

_kw(ex) = (isa(ex, Expr) && ex.head == :(=)) ? Expr(:kw, ex.args...) : ex

store_kws(; kwargs...) = kwargs

macro y(s, y, args...)
    if @capture(y, :(sym_))
        p = y
    else
        p = helper_replace_anon_func_syntax(y)
    end
    Expr(:call, :_y, esc(s), esc(p), (esc(_kw(arg)) for arg in args)...)
end

function _y(s::Selector, f, args...; kwargs...)
    if isa(f, Symbol)
        s2 = f == :custom && length(args)>1 ? replace_selector(s, args[2], :y) : s
		s2.kw[:f] = f == :custom ? args[1] : f
    else
        s2 = replace_selector(s, f, :y)
        s2.kw[:f] = :locreg
		s2.kw[:axis_type] = get(s2.kw, :axis_type, :pointbypoint)
		(length(args) > 0) && (s2.kw[:yreduce] = functionize(args[1]))
    end
    s2.kw[:fkwargs] = store_kws(; kwargs...)
    return s2
end

function _y(s::ColumnSelector, f, args...; kwargs...)
    if get(s.kw, :axis_type, :pointbypoint) != :pointbypoint
		s2 = (f == :locreg) ? replace_selector(s, args[1], :y) : s
		s2.kw[:f] = f
    else
		s2 = replace_selector(s, f, :y)
		s2.kw[:f] = :locreg
		s2.kw[:axis_type] = get(s2.kw, :axis_type, :pointbypoint)
        s2.kw[:axis_type] == :auto && (s2.kw[:axis_type] = :pointbypoint)
		(length(args) > 0) && (s2.kw[:yreduce] = functionize(args[1]))
    end
    s2.kw[:fkwargs] = store_kws(; kwargs...)
    return s2
end

macro xy(s, arg)
	return esc(:(GroupedErrors.@y(GroupedErrors.@x($s, $arg), $arg)))
end

macro xy(s, arg, f)
	return esc(:(GroupedErrors.@y(GroupedErrors.@x($s, $arg, $f), $arg, $f)))
end

macro compare(s, f)
	anon_func = helper_replace_anon_func_syntax(f)
	Expr(:call, :_compare, esc(s), esc(anon_func))
end

function _compare(s::AbstractSelector, f)
	s2 = replace_selector(s, f, :compare)
	s2.kw[:compare] = true
	s2
end

macro summarize(s, trend, variation)
    Expr(:call, :_summarize, esc(s), esc(trend), esc(variation))
end

function _summarize(s::AbstractSelector, trend, variation)
    s.kw[:summarize] = (trend, variation)
    return s
end

macro set_attr(s, arg, f)
    func = helper_replace_anon_func_syntax(f)
    Expr(:call, :_set_attr, esc(s), esc(arg), esc(func))
end

function _set_attr(s::AbstractSelector, arg, f)
    s.kw[:plot_kwargs] = get(s.kw, :plot_kwargs, [])
    push!(s.kw[:plot_kwargs], (arg, f))
	s
end

for f in (:_splitby, :_across, :_bootstrap, :_x, :_y, :_summarize, :_compare, :_set_attr)
    @eval ($f)(s, args...; kwargs...) = ($f)(convert(Selector, s), args...; kwargs...)
end


text:
struct P2P2BearingReuse
  measvec::Vector{Float64}
  predvec::Vector{Float64}
  resid::Vector{Float64}
  P2P2BearingReuse() = new(zeros(2),zeros(2),zeros(2))
end

"""
    $TYPEDEF

Single dimension bearing constraint from Pose2 to Point2 variable.
"""
struct Pose2Point2Bearing{B <: IIF.SamplableBelief} <: IIF.AbstractManifoldMinimize
    bearing::B
    reuse::Vector{P2P2BearingReuse}
    Pose2Point2Bearing{B}() where B = new{B}()
    Pose2Point2Bearing{B}(x1::B) where {B <: IIF.SamplableBelief} = new{B}(x1, [P2P2BearingReuse() for i in 1:Threads.nthreads()])
end
Pose2Point2Bearing(x1::B) where {B <: IIF.SamplableBelief} = Pose2Point2Bearing{B}(x1)

getManifold(::Pose2Point2Bearing) = SpecialOrthogonal(2)

function getSample(cfo::CalcFactor{<:Pose2Point2Bearing})
  return rand(cfo.factor.bearing)
end

function (cfo::CalcFactor{<:Pose2Point2Bearing})(Xc, p, l)
  #embed l in SE2
  M = SpecialEuclidean(2)
  q = ProductRepr(l, identity_element(SpecialOrthogonal(2), p.parts[2]))
  x,y = Manifolds.compose(M, inv(M, p), q).parts[1]
  # TODO Xc is a coordinate (ie angle), maybe change to X ϵ so2
  # m̂ = exp(so{N}(hat(SpecialOrthogonal(N), SO{N}()[], atan(y, x))))
  # distance(m, m̂)/sqrt(2)

  # FIXME, this should be a vee(log()), not linear - and add to test
  return  Manifolds.sym_rem(Xc[1] - atan(y, x))
end
# define the conditional probability constraint
# function (cfo::CalcFactor{<:Pose2Point2Bearing})(meas, _xi, lm)
#   #
#   M = SpecialEuclidean(2)
#   ϵ = identity_element(M)
#   xi = vee(M, ϵ, log(M, ϵ, _xi))

#   reuse = cfo.factor.reuse[Threads.threadid()]
#   reuse.measvec[1] = cos(meas[1] + xi[3])
#   reuse.measvec[2] = sin(meas[1] + xi[3])

#   # @simd for i in 1:2
#     reuse.predvec[1:2] .= lm[1:2].-xi[1:2]
#   # end
#   normalize!(reuse.predvec)
#   reuse.resid .= reuse.measvec
#   reuse.resid .-= reuse.predvec

#   # must return result of length zDim==1 in this case
#   return sum(abs.(reuse.resid))
#   # reuse.resid .^= 2
#   # res[1] = reuse.resid[1]
#   # res[1] += reuse.resid[2]

# end

# Packing and Unpacking
mutable struct PackedPose2Point2Bearing <: IncrementalInference.PackedInferenceType
    bearstr::String
    PackedPose2Point2Bearing() = new()
    PackedPose2Point2Bearing(s1::AS) where {AS <: AbstractString} = new(string(s1))
end
function convert(::Type{PackedPose2Point2Bearing}, d::Pose2Point2Bearing{B}) where {B <: IIF.SamplableBelief}
  return PackedPose2Point2Bearing(convert(PackedSamplableBelief, d.bearing))
end
function convert(::Type{Pose2Point2Bearing}, d::PackedPose2Point2Bearing)
  Pose2Point2Bearing(convert(SamplableBelief, d.bearstr))
end