text:
#=
  http://cplint.eu/example/inference/viral.swinb
  """
  Viral Marketing

  A firm is interested in marketing a new product to its customers. These are connected in a social network
  that is known to the firm: the network represents the trust relationships between customers. The firm has
  decided to adopt a marketing strategy that involves giving the product for free to a number of its
  customers, in the hope that these influence the other customers and entice them to buy the product.
  The firm wants to choose the customers to which marketing is applied so that its return is maximized.
  This involves computing the probability that the non-marketed customers will acquire the product given the
  action to the marketed customers.

  This viral marketing scenario is inspired by [1].

  We can model this domain with an LPAD where the predicate trust/2 encodes the links between customers in the
  social network and predicate has/1 is true of customers that possess the product, either received as a gift
  or bought. Predicate trust/2 is defined by a number of certain facts, while predicate has/1 is defined by
  two rules, one expressing the prior probability of a customer to buy the product and one expressing the
  fact that if a trusted customer has the product, then there is a certain probability that the trusting
  customer buys the product.
  """

  Reference G. V. den Broeck, I. Thon, M. van Otterlo, L. D. Raedt,
  Dtproblog: A decision-theoretic probabilistic prolog, in: M. Fox, D. Poole (Eds.),
  24th AAAI Conference on Artificial Intelligence, AAAI’10, Atlanta, Georgia, USA,
  July 11-15, 2010, AAAI Press, 2010, pp. 1217–1222.

  """
  We want to compute the probability that customer 2 buys the product if we perform the action of
  giving the product to customer 3. We need to use causal reasoning so we use the action do(has(3))
  as evidence:
  ?- prob(has(2),do(has(3)),P).

  P = 0.136

  If instead we compute the effect of the action using regular probabilistic inference, we get:
  ?- prob(has(2),has(3),P).
  P = 0.4065135474609725
  So not distinguishing seeing from doing leads to an overly optimistic estimate.

  ""

  Cf ~/blog/viral_marketing.blog
     ~/psi/viral_marketing.psi
     ~/webppl/viral_marketing.wppl

  Distribution of values for has(p*)
  The BLOG model (-n 100000): given has(p3)
                  has(p1): 0.25285605483625506
                  has(p2): 0.32750708813609203
                  has(p3): 1.0
                  has(p4): 0.39824124623193263

  The PSI model: given has(p3)
                 has(p1): 0.254777070063694
                 has(p2): 0.331210191082803,
                 has(p3): 1,
                 has(p4): 0.4

  This Turing model,using
     chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)

  Summary Statistics
  parameters      mean       std   naive_se      mcse           ess      rhat
      Symbol   Float64   Float64    Float64   Float64       Float64   Float64

    has_a[1]    0.2512    0.4337     0.0007    0.0033    11421.0186    1.0006
    has_a[2]    0.3407    0.4739     0.0007    0.0035    10540.0175    1.0002
    has_a[3]    1.0000    0.0045     0.0000    0.0000   122493.4137    1.0000
    has_a[4]    0.3973    0.4893     0.0008    0.0034    12741.2279    1.0003

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

    has_a[1]    0.0000    0.0000    0.0000    1.0000    1.0000
    has_a[2]    0.0000    0.0000    0.0000    1.0000    1.0000
    has_a[3]    1.0000    1.0000    1.0000    1.0000    1.0000
    has_a[4]    0.0000    0.0000    0.0000    1.0000    1.0000

Distributions of variable (num:0)
(0.0, 0.0, 1.0, 0.0)    =>  122148 (0.30537)
(0.0, 0.0, 1.0, 1.0)    =>  81756 (0.20439)
(0.0, 1.0, 1.0, 0.0)    =>  56870 (0.142175)
(0.0, 1.0, 1.0, 1.0)    =>  38753 (0.0968825)
(1.0, 0.0, 1.0, 0.0)    =>  36105 (0.0902625)
(1.0, 1.0, 1.0, 0.0)    =>  25975 (0.0649375)
(1.0, 0.0, 1.0, 1.0)    =>  23722 (0.059305)
(1.0, 1.0, 1.0, 1.0)    =>  14663 (0.0366575)
(0.0, 0.0, 0.0, 1.0)    =>  4 (1.0e-5)
(1.0, 0.0, 0.0, 1.0)    =>  3 (7.5e-6)
(1.0, 0.0, 0.0, 0.0)    =>  1 (2.5e-6)
 21.266611 seconds (229.77 M allocations: 11.048 GiB, 23.33% gc time)

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function viral_marketing()
    p1,p2,p3,p4 = 1:4
    people = [p1,p2,p3,p4]

    trustMatrix = [[0,0,0,0], # p1
                   [1,0,0,0], # p2
                   [1,1,0,0], # p3
                   [1,0,1,0]  # p4
                  ]

    function trusts(a, b)
        return trustMatrix[a][b]
    end

    has_a = tzeros(4)
    #
    # Note: It seems that one have to create the extra array has_a
    #       to get the proper probabilities for has(p1..p4)
    #       Using
    #         # ...
    #         if s > 0
    #           t ~ flip(0.4)
    #         else
    #           t ~ flip(0.1)
    #         end
    #         return t
    #       don't differentiate between the p's.
    #       And this means that we can convert this to a for loop...
    #
    #=
    function has(p)
        s = sum(Int64[trusts(p,q)==1 && has(q)==1 for q in people])
        if s > 0
            has_a[p] ~ flip(0.4)
        else
            has_a[p] ~ flip(0.1)
        end
        return has_a[p]
    end
    =#

    # This is a little faster than using has(p)
    for p in people
        # s = sum(Int64[trusts(p,q)==1 && has_a[q]==1 for q in people])
        s = count(q->trusts(p,q)==1 && has_a[q]==1, people)
        if s > 0
            has_a[p] ~ flip(0.4)
        else
            has_a[p] ~ flip(0.1)
        end
    end

    # true ~ Dirac(has(p3)==1)
    # return has(p1),has(p2),has(p3),has(p4)
    true ~ Dirac(has_a[p3]==1)
    return has_a[p1],has_a[p2],has_a[p3],has_a[p4]

end

model = viral_marketing()

num_chains = 4

# chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, MH(), 10_000)

# chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)

# chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 10_000)
# chains = sample(model, IS(), 10_000)
#
display(chains)

genq = generated_quantities(model, chains)
show_var_dist_pct(genq)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text: abstract type AExpressionOperator <: AOperator end

ext: jl


text:
using Plots
using Random
using Colors

# seed
# Random.seed!(3)

random_tenosr_dist = "uniform"
#random_tenosr_dist = "gauss"

# experimental setup
rep_max = 10

# zero cutoff KL divergence
eps_thred = 1.0e-12

# save file name option
n_digits = 5

# Real datasets
datasets = [
            "AttFace",      #(112,92,400)
           #"4DLFD"    #(9,9,512,512,3)
           ]
sort!(datasets)

dataset_label_dict = Dict(
                   "AttFace" => "AttFace",
                   #"4DLFD" => "4DLFD",
                )
dataset_label_dict = sort(dataset_label_dict)

# approximation methods
methods_flags = Dict(
             "LT1R" => false,
             "LTR" => true,
             "NTD_KL" => true,
             "NTD_LS" => true,
             "lraSNTD" => true,
            )
non_negative_methods = ["LT1R", "LTR", "NTD_KL", "NTD_LS", "lraSNTD"]
methods = [k for k = keys(methods_flags) if methods_flags[k] == true]

# move to "TLrRR" to top
sort!(methods)
methods = insert!( setdiff(methods, ["LTR"]), 1, "LTR" )
std_method = "LTR"

# #################################
# Experimental design for real data
# #################################

plans = Dict( dataset for dataset = zip(datasets, [ [] for i = 1:length(datasets) ] ) )
# AttFace (112,92,400)
plans["AttFace"] =
[[1,1,1], [3,3,3], [5,5,5], [10,10,10], [15,15,15], [20,20,20], [30,30,30], [40,40,40], [50,50,50], [60,60,60], [70,70,70], [80,80,80]]

#"4DLFD"    #(9,9,512,512,3)
plans["4DLFD"] =
[[1,1,1,1,1], [2,2,2,2,1],[3,3,4,4,1],[3,3,5,5,1],[3,3,6,6,1],[3,3,7,7,1],[3,3,8,8,1],[3,3,16,16,1],[3,3,20,20,1],
 [3,3,40,40,1],[3,3,60,60,1],[3,3,80,80,1]]

# ##########################
# config for other methods
# ##########################

NTD_init = "SVD"

# ##########################
# config for Plots
# ##########################

# plt options for AI data
# fonts list : Helvetica, palatino, courier, times, newcentryschlbk, advantgrade
img_size = (500, 500)
fnt1 = font(20, "times")
fnt2 = font(12, "times")
fnt3 = font(15, "times")

linewidth = 2.5
markerstrokewidth = 1.0
markersize = 12
markerwidth = 2

legend = :none
#legend = :bottomleft
legendfont = fnt1

lw_y0 = 3

# plot options for real data
img_size_wide = (1050, 300)

# color setting
colors = palette(:tab20)

colors_dict = Dict(
                   "LT1R" => :red,
                   "TLrRR" => :red,
                   "NTD_KL" => :black,
                   "NTD_LS" => :blue,
                   "lraSNTD" => :brown,
                )

label_dict = Dict(
                   "LT1R" => "TL1R",
                   "LTR" => "LTD(proposed)",
                   "NTD_KL" => "NTD_KL",
                   "NTD_LS" => "NTD_LS",
                   "lraSNTD" => "lraSNTD",
                )

markershapes_dict = Dict(
                   "LT1R" => :circle,
                   "TLrRR" => :circle,
                   "NTD_KL" => :rect,
                   "NTD_LS" => :dtriangle,
                   "lraSNTD" => :utriangle,
                )

linetype_dict = Dict(
                   "LT1R" => :dashdot,
                   "LTR" => :dashdot,
                   "NTD_KL" => :dashdot,
                   "NTD_LS" => :dot,
                   "lraSNTD" => :dashdot,
                )

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
######################################################################
# constant_initialier.jl
######################################################################
"""
`ci = ConstantInitializer(action::Vector{Float64})`

Just takes a constant action.
"""
mutable struct ConstantInitializer <: Initializer
	action::Vector{Float64}
end

function initialize(ci::ConstantInitializer, em::ErgodicManager, tm::TrajectoryManager)
	xd =  Array{Vector{Float64}}(undef, tm.N+1)#1.0 .*ones(tm.N+1,2,2)#Array{Vector{Float64}}(tm.N+1)
	ud =  Array{Vector{Float64}}(undef, tm.N)#1.0 .*ones(tm.N,2,2)#Array{Vector{Float64}}(tm.N)
	xd[1] = deepcopy(tm.x0)
	ud[1] = deepcopy(ci.action)
	for i = 1:(tm.N-1)
		xd[i+1] = integrate(tm, xd[i], ud[i])
		ud[i+1] = deepcopy(ci.action)
	end
	xd[tm.N+1] = integrate(tm, xd[tm.N], ud[tm.N])

	return xd, ud
end


text:
"""
    RationalQuadratic(lnℓ::AbstractFloat) <: RadialKernel{SqEuclidean}

The rational quadratic kernel with parameters
``\\exp(2 \\ln\\ell) = \\ell^2 > 0`` and ``\\exp(\\ln\\alpha) = \\alpha > 0``.

The rational quadratic kernel may be used to model functions with a
variable length scale.

```math
k_{\\ell,\\alpha}(x, y) = k_{\\ell,\\alpha}(\\lVert x - y\\rVert_2^2) =
\\bigg(1 + \\frac{\\lVert x - y\\rVert_2^2}{2\\alpha\\ell^2}\\bigg)^{-\\alpha}
```

External links
* [Rational quadratic covariance function on Wikipedia](https://en.wikipedia.org/wiki/Rational_quadratic_covariance_function)

"""
mutable struct RationalQuadratic{
    F<:AbstractFloat
} <: RadialKernel{SqEuclidean}
    dist::SqEuclidean
    lnℓ::F
    lnα::F

    α::F
    twoαℓ²::F
    function RationalQuadratic(
        dist::SqEuclidean,
        lnℓ::AbstractFloat,
        lnα::AbstractFloat
    )
        lnℓ, lnα = promote(lnℓ, lnα)
        new{typeof(lnℓ)}(
            dist,
            lnℓ,
            lnα,
            exp(lnα),
            exp(2 * lnℓ + lnα + float(logtwo)),
        )
    end
end

RationalQuadratic(lnℓ::AbstractFloat, lnα::AbstractFloat) =
    RationalQuadratic(SqEuclidean(), lnℓ, lnα)

@inline (k::RationalQuadratic)(τ::AbstractFloat) = (1 + τ / k.twoαℓ²)^(-k.α)

numparams(::RationalQuadratic) = (1, 1)
paramtypes(::RationalQuadratic{F}) where F = (F, F)
params(k::RationalQuadratic) = (lnℓ = k.lnℓ, lnα = k.lnα)
function setparams!(k::RationalQuadratic{F}, lnℓ::F, lnα::F) where F
    k.lnℓ = lnℓ
    k.lnα = lnα
    k.α = exp(lnα)
    k.twoαℓ² = exp(2 * lnℓ + lnα + float(logtwo))
end


text:
# Day 0: Hello, World.
# Read a full line of input from stdin and print it out.
input_string = readline()

# Print a string literal saying "Hello, World." to stdout.
println("Hello, World.")

# TODO: Write a line of code here that prints the contents of input_string to stdout.
println(input_string)

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 00.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
@doc """
    CenteredDifference{N}(n, order, step, len, [coeff_func])
    CenteredDifference{N}(n, order, steps, len, [coeff_func])

See also: [`UpwindDifference`](@ref)
"""
CenteredDifference

@doc """
    calculate_weights(n::Int, x₀::Real, x::Vector)

Return a vector `c` such that `c⋅f.(x)` approximates ``f^{(n)}(x₀)`` for a smooth `f`.

The points `x` need not be evenly spaced.

The stencil `c` has the highest approximation order possible given values of `f` at `length(x)` points. More precisely, if `x` has length `m`, there is a function `g` such that ``g(y) = f(y) + O(y-x₀)^{m-n+?}`` and ``c⋅f.(x) = g'(x₀)``.

The algorithm is due to [Fornberg](https://doi.org/10.1090/S0025-5718-1988-0935077-0), with a [modification](http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507) to improve stability.
"""
calculate_weights

@doc """
    DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SVector,T3,F}(<fields>)

Represent a finite-difference derivative operator.

These operators implement the `DiffEqBase.AbstractDiffEqLinearOperator` interface. Therefore, `eltype` returns `T`.

These operators can be contracted over an arbitrary dimension, given by the type parameter `N`.

The `CentredDifference` and `UpwindDifference` types serve as hooks to hang constructors, but no objects are actually constructed with those types. Their constructors return structures of type `DerivativeOperator`, with stencils constructed by `calculate_weights`, and other fields filled out appropriately.

The finite-difference methods are defined for `DerivativeOperator`. In particular, `*`, which actually takes derivatives.

The key data are three stencils of coefficients, stored in `SVector`s. The interior stencil, stored in `stencil_coefs::S1`, is the normal one used in the interior of the grid. The others, `low_boundary_coefs::S2` and `high_boundary_coefs::S2`, are used where the normal stencil would jut out of the grid boundary. These can have a different length than the interior stencil, hence the two types.

When the operator is applied by `mul!`, these stencils are multiplied by the vector in the `coefficients` field. It is set to `coeff_func.(1:len)` in `UpwindDifference`, which seems odd, because `len` is not the stencil length. Scalar multiplication is absorbed into this vector, as required by the `DiffEqBase.AbstractDiffEqLinearOperator` interface.

The `coefficients` field appears to be a more general `DifferentialEquations` thing. There is an `update_coefficients!` method.

The term “left boundary” is used interchangeably with “low boundary”, and “right” with “high”.

# Type Parameters

Presumably, most of these are to force method specialization on:a) the stencil lengths, which are given by the parameters of the `SVector` types `S1` and `S2`, and b) even or arbitrary grid spacing.

- `T<:Real`: Function range type.
- `N`: Contraction dimension.
- `Wind`: Flag `true` for upwind operators, `false` for centered.
- `T2`: Concrete type of `dx` field, which might be `T` or `AbstractVector{T}`.
- `S1`: Concrete type of the interior stencil. (When is this not an `SVector`?)
- `S2<:SVector`: Concrete type of boundary coefficients.
- `T3`: `Vector{T}` if a `coeff_func` was supplied, `Nothing` otherwise.
- `F`: `typeof(coeff_func)` if a `coeff_func` was supplied, `Nothing` otherwise.

# Fields

- `derivative_order::Int`: The `n` in `f⁽ⁿ⁾(x)`.
- `approximation_order::Int`: The degree of polynomial for which the operator is numerically exact.
- `dx::T2`: Grid step.
- `len::Int`: The number of interior points on the grid.
- `stencil_length::Int`: Length of the interior stencil.
- `stencil_coefs::S1`: An `SVector` of the interior stencil coefficients.
- `boundary_stencil_length::Int`:
- `boundary_point_count::Int`:
- `low_boundary_coefs::S2`: An array of `SVector` stencils for the interior points where `stencil_coefs` juts over the left boundary.
- `high_boundary_coefs::S2`: An array of `SVector` stencils for the interior points where `stencil_coefs` juts over the right boundary.
- `coefficients::T3`: Multiplier for the stencil.
- `coeff_func::F`: God only knows what this function is for.

See also: [`AbstractBC`](@ref), [`CenteredDifference`](@ref), [`UpwindDifference`](@ref)
"""
DerivativeOperator

@doc """
    abstract type AbstractBC{T}

`T` is the range type of the discretized function.

Boundary condition operators extrapolate the discretized function
as illustrated in the README, adding a ghost node at each end such
that an interpolated polynomial satisfies the boundary condition.

Some of these are affine operators with a constant term, not linear
operators, so using `*` and `\\` is a minor abuse of notation. For example,
when there is a non-zero Dirichlet boundary value, the value at the
ghost node is constant, and does not scale proportional to the
interior function values. Because boundary condition operators are
not linear, they can not be concretized as matrices.

`GhostDerivativeOperator` somehow concretizes as a pair of matrices.

The implementation of `*(::AffineBC, ::Vector)` returns a
`BoundaryPaddedVector`. This is a structure that stores the ghost
values separately from the original vector. Presumably, this is to
avoid allocations, although the `mul!` method for `GhostDerivativeOperator`
and `Matrix` method allocates a new `BoundaryPaddedArray` for each
column.

In higher-dimensional cases, it might return a
`BoundaryPaddedArray`.

The `*` operator for a derivative multiplying a boundary condition
packages them up as a `GhostDerivativeOperator`. This cannot
return a sparse matrix, because the operator is affine (not linear).
When a vector is multiplied by one of those, it uses `*(::AbstractBC,
::AbstractVector)` to pad with ghost nodes, then
`mul!(::DerivativeOperator, ::AbstractVector)` evaluates the
derivative.

See also: [`DerivativeOperator`](@ref), [`PeriodicBC`](@ref) [`NeumannBC`](@ref), [`DirichletBC`](@ref), [`RobinBC`](@ref).
"""
AbstractBC

@doc """
    PeriodicBC{T}()

q = PeriodicBC{T}()

Qx, Qy, ... = PeriodicBC{T}(size(u)) #When all dimensions are to be extended with
a periodic boundary condition.

Creates a periodic boundary condition, where the lower index end of some u is extended
with the upper index end and vice versa. It is not recommended to concretize this BC
type in to a BandedMatrix, since the vast majority of bands will be all 0s.
SparseMatrix concretization is recommended.

Currently, periodic boundary conditions are only implemented in one
dimension.

# Type parameters

- `T` is the domain type of the discretized function.

See also: [`AbstractBC`](@ref)
"""
PeriodicBC

"""
    q = RobinBC(l, r, dx::T, approximation_order) where T # When this BC extends a dimension with a uniform step size

    q = RobinBC(l, r, dx::Vector{T}, approximation_order) where T # When this BC extends a dimension with a non-uniform step size. dx should be the vector of step sizes for the whole dimension.

`l` and `r` are the BC coefficients, i.e., `(αl, βl, γl)` and `(αl, βl, γl)` (tuples and vectors work)
and correspond to BCs of the form:

`αl * u + βl * u' = γl`
`αr * u + βr * u' = γr`

imposed on the lower (`l`) and higher (`r`) index boundaries, respectively.

`RobinBC` implements a Robin boundary condition operator `Q` that acts on a vector to give an extended
vector as a result (see https://github.com/JuliaDiffEq/DiffEqOperators.jl/files/3267835/ghost_node.pdf).

Write vector b̄₁ as a vertical concatenation with b0 and the rest of the elements of b̄₁, denoted b̄`₁, the same with ū into u0 and ū`. b̄`₁ = b̄`_2 = fill(β/Δx, length(stencil)-1).
Pull out the product of u0 and b0 from the dot product. The stencil used to approximate u` is denoted s. b0 = α+(β/Δx)*s[1].
Rearrange terms to find a general formula for u0:= -b̄`₁̇⋅ū`/b0 + γ/b0, which is dependent on ū` the robin coefficients and Δx.
The non-identity part of Qa is qa:= -b`₁/b0 = -β.*s[2:end]/(α+β*s[1]/Δx). The constant part is Qb = γ/(α+β*s[1]/Δx)
do the same at the other boundary (amounts to a flip of s[2:end], with the other set of boundary coefficients).
"""
RobinBC
@doc (@doc RobinBC) NeumannBC
@doc (@doc RobinBC) Neumann0BC
@doc (@doc RobinBC) DirichletBC
@doc (@doc RobinBC) Dirichlet0BC


text:
using ProgressMeter

include("./plot_ic_reward.jl")
include("./plot_ic_goal_visitation.jl")

# Experiment parameters
log_interval = 100
num_steps = 60000

# Getting baseline TB
folder = "M:/globus/InterestTB/Experiment2_dpi/"
ic = ItemCollection(folder)

# println(diff(ic))
# # asdasd
search_ic = search(ic, Dict("demon_learner" => "SR", "demon_update" => "TB", "behaviour_learner" => "GPI"))
gpi_tb_ic = get_best(search_ic, ["eta",  "alpha_init"], :ttmaze_direct_error)

search_ic = search(ic, Dict("demon_learner" => "SR", "demon_update" => "TB", "behaviour_learner" => "Q"))
esarsa_tb_ic = get_best(search_ic, ["eta",  "alpha_init"], :ttmaze_direct_error)

configs = [
    Dict(
        :search => Dict("demon_learner" => "Q", "demon_update" => "ETB"),
        :label => "Q ETB"
    )
    Dict(
        :search => Dict("demon_learner" => "Q", "demon_update" => "EmphESARSA"),
        :label => "Q EmphESARSA"
    )
    Dict(
        :search => Dict("demon_learner" => "Q", "demon_update" => "InterestTB"),
        :label => "Q InterestTB"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "ETB"),
        :label => "SR ETB"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "EmphESARSA"),
        :label => "SR EmphESARSA"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "InterestTB"),
        :label => "SR InterestTB"
    )
]

folder = "M:/globus/InterestTB/Experiment2_ESARSA/"
ic = ItemCollection(folder)

p = plot()

@showprogress for config in configs
    local search_ic = search(ic, config[:search])
    local best_ic = get_best(search_ic, ["demon_alpha_init", "eta", "behaviour_alpha_init"], :ttmaze_direct_error)
    plot_ic_reward!(p, best_ic, :ttmaze_direct_error, log_interval; smooth_step=1, label=config[:label])
end

plot_ic_reward!(p, esarsa_tb_ic, :ttmaze_direct_error, log_interval; smooth_step=1, label="SR TB")
plot!(p, title="Q")

savefig("plotting/chunlok/generated_plots/Experiment2_rewards_esarsa.svg")

##########################################################################################
# Plotting configs with GPI behaviour
##########################################################################################

folder = "M:/globus/InterestTB/Experiment2_GPI/"
ic = ItemCollection(folder)

# # println(diff(ic))

# Plotting the RMSE plots for all configs with GPI behaviour
p = plot()
@showprogress for config in configs
    local search_ic = search(ic, config[:search])
    local best_ic = get_best(search_ic, ["demon_alpha_init", "eta", "behaviour_alpha_init"], :ttmaze_direct_error)
    plot_ic_reward!(p, best_ic, :ttmaze_direct_error, log_interval; smooth_step=1, label=config[:label])
end

plot_ic_reward!(p, gpi_tb_ic, :ttmaze_direct_error, log_interval; smooth_step=1, label="SR TB")
plot!(p, title="GPI")
savefig("plotting/chunlok/generated_plots/Experiment2_rewards_gpi.svg")

# Plotting goal visitations for each config for GPI.
# for config in configs
#     local p = plot()
#     local search_ic = search(ic, config[:search])
#     local best_ic = get_best(search_ic, ["demon_alpha_init", "eta", "behaviour_alpha_init"], :ttmaze_direct_error)
#     plot_ic_goal_visitation!(p, best_ic)
#     plot!(title=config[:label])
#     savefig("plotting/chunlok/generated_plots/Experiment2_GPI_goal_visitation_$(config[:label])")
# end

# Plotting reward per GVF
# best_ics = []
# for config in configs
#     local search_ic = search(ic, config[:search])
#     local best_ic = get_best(search_ic, ["demon_alpha_init", "eta", "behaviour_alpha_init"], :ttmaze_direct_error)
#     push!(best_ics, best_ic)
# end

# for gvf_i in 1:4
#     p = plot()
#     for (config_i, config) in enumerate(configs)
#         plot_ic_reward_per_gvf!(p, best_ics[config_i], :ttmaze_direct_error, log_interval, num_steps, gvf_i; smooth_step=1, label=config[:label])
#     end
#     plot!(title="GVF $(gvf_i)")
#     savefig("plotting/chunlok/generated_plots/Experiment2_GPI_rewards_gvf_$(gvf_i)")
# end

sr_configs = [
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "ETB"),
        :label => "SR ETB"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "EmphESARSA"),
        :label => "SR EmphESARSA"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "InterestTB"),
        :label => "SR InterestTB"
    )
]

# SR reward plots
# p = plot()
# folder = "M:/globus/InterestTB/Experiment2_GPI/"
# ic = ItemCollection(folder)
# @showprogress for config in sr_configs
#     local search_ic = search(ic, config[:search])
#     local best_ic = get_best(search_ic, ["demon_alpha_init", "eta", "behaviour_alpha_init"], :ttmaze_direct_error)
#     plot_ic_reward!(p, best_ic, :ttmaze_direct_error, log_interval, num_steps; smooth_step=1, label="GPI " * config[:label])
# end

# folder = "M:/globus/InterestTB/Experiment2_ESARSA/"
# ic = ItemCollection(folder)

# @showprogress for config in sr_configs
#     local search_ic = search(ic, config[:search])
#     local best_ic = get_best(search_ic, ["demon_alpha_init", "eta", "behaviour_alpha_init"], :ttmaze_direct_error)
#     plot_ic_reward!(p, best_ic, :ttmaze_direct_error, log_interval, num_steps; smooth_step=1, label="ESARSA " * config[:label])
# end

# savefig("plotting/chunlok/generated_plots/Experiment2_rewards_SR.svg")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
export batchqromp

# Write a batch version of the QR-OMP code

# Should not be considered production-ready

function batchqromp(ψ::AbstractMatrix{T}, u::AbstractVector{T}; invert::Bool=true, verbose::Bool = true, ϵrel::Float64 = 1e-1, batchterms::Int64 = 1, maxterms::Int64=typemax(Int64)) where {T}
    m, n = size(ψ)

    residue = copy(u)
    idxset = Int64[]
    dict = collect(1:n)
    ϵu = norm(u)
    ϵhist = [ϵu]
    ϵrel *= ϵu
    ϕ = zeros(T, n)
    ratio = zeros(T,n)
    cache = zeros(T, m, batchterms)

    @inbounds for j=1:n
        for i=1:m
        ϕ[j] += ψ[i,j]^2
        end
    end

    cache = zeros(T, m, batchterms)
    F = elasticqrfact(cache)
    ek = zeros(T, n, batchterms)
    q = zeros(T, m, batchterms)
    new_idx = zeros(Int64, batchterms)

    @inbounds for k=1:div(n, batchterms)
        # Step  7 choose best candidate

        @fastmath @inbounds for (j, idx) in enumerate(dict)
            ψj = view(ψ,:,idx)
            ratio[idx] = dot(residue, ψj)^2/ϕ[idx]

            # if ratio >= entry
            #     entry = copy(ratio)
            #     new_idx = copy(idx)
            # end
        end
        new_idx .= partialsortperm(view(ratio, dict), 1:batchterms, rev = true)

        # Step 8 Update set of selected basis
        append!(idxset, new_idx)

        # Step 9 Update candidate dictionary
        filter!(x-> x != new_idx, dict)

        # Step 10 Update Q(k-1) and R(k-1) with ψi(k)
        if k>1
            # F = updateqrfactUnblocked!(F, view(ψ,:,new_idx))
            @time updateelasticqrfact!(F, view(ψ,:,new_idx), cache)
            @inbounds for i=1:batchterms
                ek[(k-2)*batchterms+i,i] = 0.0
            end
        else
            # F = qrfactUnblocked(ψ[:,new_idx:new_idx])
            F = elasticqrfact(view(ψ,:,new_idx))
        end

        # Step 11 Update residual
        @inbounds for i=1:batchterms
            ek[(k-1)*batchterms+i,i] = 1.0
        end

        mul!(q, F.Q, view(ek,1:k,1:batchterms))

        @inbounds for (j, idx) in enumerate(dict)
            ψj = view(ψ,:,idx)
            ϕj = dot(q, ψj)
            ϕ[idx] -= ϕj^2
        end
        factor = dot(q, u)
        residue -= factor*q

        # Step 12 Calculate stopping critera
        ϵ = norm(residue)

        if verbose == true
            push!(ϵhist, copy(ϵ))
        end
        if ϵ < ϵrel || k == maxterms #|| size(idxset,1) + batchterms > n
            break
        end
    end

    # Solve the system with the set of indices
    if verbose == true
        if invert == true
            c = zeros(length(idxset))
            ldiv!(c, F, u)
            return idxset, c, ϵhist
        else
            return idxset, ϵhist
        end
    else
        if invert == true
            c = zeros(length(idxset))
            ldiv!(c, F, u)
            return idxset, c
        else
            return idxset
        end
    end
end

max_forks_repo_name: mleprovost/QROMP


text:
module Channel

include("Channel/Gilbert.jl")

end


text:
import QuantumControlBase
using QuantumControlBase: getcontrols, getcontrolderivs, discretize_on_midpoints, evalcontrols
using QuantumPropagators: init_storage
using ConcreteStructs

# Krotov workspace (for internal use)
@concrete terse struct KrotovWrk

    # a copy of the objectives
    objectives

    # the adjoint objectives, containing the adjoint generators for the
    # backward propagation
    adjoint_objectives

    # The kwargs from the control problem
    kwargs

    # Tuple of the original controls (probably functions)
    controls

    # storage for controls discretized on intervals of tlist
    pulses0 :: Vector{Vector{Float64}}

    # second pulse storage: pulses0 and pulses1 alternate in storing the guess
    # pulses and optimized pulses in each iteration
    pulses1 ::  Vector{Vector{Float64}}

    # values of ∫gₐ(t)dt for each pulse
    g_a_int :: Vector{Float64}

    # update shapes S(t) for each pulse, discretized on intervals
    update_shapes ::  Vector{Vector{Float64}}

    lambda_vals :: Vector{Float64}

    is_parametrized :: Vector{Bool}

    parametrization

    # map of controls to options
    pulse_options

    # Result object

    result

    #################################
    # scratch objects, per objective:

    # backward-propagated states
    # note: storage for fw-propagated states is in result.states
    bw_states

    # dynamical generator at a particular point in time
    G

    control_derivs :: Vector{Vector{Union{Function, Nothing}}}

    vals_dict

    fw_storage # forward storage array (per objective)

    fw_storage2 # forward storage array (per objective)

    bw_storage # backward storage array (per objective)

    fw_prop_wrk

    bw_prop_wrk

    use_threads :: Bool

end

function KrotovWrk(problem::QuantumControlBase.ControlProblem)
    use_threads = get(problem.kwargs, :use_threads, false)
    objectives = [obj for obj in problem.objectives]
    adjoint_objectives = [adjoint(obj) for obj in problem.objectives]
    controls = getcontrols(objectives)
    control_derivs = [
        getcontrolderivs(obj.generator, controls) for obj in objectives
    ]
    tlist = problem.tlist
    kwargs = Dict(problem.kwargs)
    # TODO: check that kwargs has all the required parameters
    pulse_options = problem.pulse_options
    update_shapes = [
        discretize_on_midpoints(
                pulse_options[control][:update_shape],
                tlist
        )
        for control in controls
    ]
    lambda_vals = [
        convert(Float64, pulse_options[control][:lambda_a])
        for control in controls
    ]
    is_parametrized = [
        haskey(pulse_options[control], :parametrization)
        for control in controls
    ]
    parametrization = [
        get(pulse_options[control], :parametrization, NoParametrization())
        for control in controls
    ]
    if haskey(kwargs, :continue_from)
        @info "Continuing previous optimization"
        result = kwargs[:continue_from]
        if !(result isa KrotovResult)
            # account for continuing from a different optimization method
            result = convert(KrotovResult, result)
        end
        result.iter_stop = get(problem.kwargs, :iter_stop, 5000)
        result.converged = false
        result.start_local_time = now()
        result.message = "in progress"
        pulses0 = [
            discretize_on_midpoints(control, tlist)
            for control in result.optimized_controls
        ]
    else
        result = KrotovResult(problem)
        pulses0 = [
            discretize_on_midpoints(control, tlist) for control in controls
        ]
    end
    pulses1 = [copy(pulse) for pulse in pulses0]
    g_a_int = zeros(length(pulses0))
    bw_states = [similar(obj.initial_state) for obj in objectives]
    zero_vals = IdDict(control => zero(pulses0[i][1]) for (i, control) in enumerate(controls))
    G = [evalcontrols(obj.generator, zero_vals) for obj in objectives]
    vals_dict = [copy(zero_vals) for _ in objectives]
    # TODO: forward_storage only if g_b != 0
    fw_storage = [init_storage(obj.initial_state, tlist) for obj in objectives]
    # TODO: second forward storage only if second order
    fw_storage2 = [init_storage(obj.initial_state, tlist) for obj in objectives]
    bw_storage = [init_storage(obj.initial_state, tlist) for obj in objectives]
    fw_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj, :fw_prop_method, :prop_method; problem.kwargs...
            )
        )
        for obj in objectives
    ]
    bw_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj, :bw_prop_method, :prop_method; problem.kwargs...
            )
        )
        for obj in objectives
    ]

    fw_prop_wrk = [
        QuantumControlBase.initobjpropwrk(obj, tlist, fw_prop_method[k];
                                            initial_state=obj.initial_state,
                                            kwargs...)
        for (k, obj) in enumerate(objectives)
    ]
    bw_prop_wrk = [
        QuantumControlBase.initobjpropwrk(obj, tlist, bw_prop_method[k];
                                            initial_state=obj.initial_state,
                                            kwargs...)
        for (k, obj) in enumerate(objectives)
    ]
    return KrotovWrk(
        objectives, adjoint_objectives, kwargs, controls, pulses0, pulses1,
        g_a_int, update_shapes, lambda_vals, is_parametrized,
        parametrization, pulse_options, result, bw_states, G,
        control_derivs, vals_dict, fw_storage, fw_storage2,
        bw_storage, fw_prop_wrk, bw_prop_wrk, use_threads
    )
end


text:
# This file is a "Hello, world!" in Julia language for wandbox.

include("test1.jl")

using .Test1

println("Hello, Wandbox!")

Test1.test()

# Julia language references:
#   https://docs.julialang.org/en/v1/

max_stars_repo_name: srz-zumix/wandbox-api

max_issues_repo_name: srz-zumix/wandbox-api

max_forks_repo_name: srz-zumix/wandbox-api


text:
using JuMP
import MathOptInterface
const MOI = MathOptInterface
include("utils.jl")

function reopt(reo_model, data, model_inputs)

    p = build_param(Tech = model_inputs["Tech"],
          Load = model_inputs["Load"],
          TechClass = model_inputs["TechClass"],
          TechIsGrid = model_inputs["TechIsGrid"],
          TechToLoadMatrix = model_inputs["TechToLoadMatrix"],
          TurbineDerate = model_inputs["TurbineDerate"],
          TechToTechClassMatrix = model_inputs["TechToTechClassMatrix"],
          NMILRegime = model_inputs["NMILRegime"],
          r_tax_owner = model_inputs["r_tax_owner"],
          r_tax_offtaker = model_inputs["r_tax_offtaker"],
          pwf_om = model_inputs["pwf_om"],
          pwf_e = model_inputs["pwf_e"],
          pwf_prod_incent = model_inputs["pwf_prod_incent"],
          LevelizationFactor = model_inputs["LevelizationFactor"],
          LevelizationFactorProdIncent = model_inputs["LevelizationFactorProdIncent"],
          StorageCostPerKW = model_inputs["StorageCostPerKW"],
          StorageCostPerKWH = model_inputs["StorageCostPerKWH"],
          OMperUnitSize = model_inputs["OMperUnitSize"],
          CapCostSlope = model_inputs["CapCostSlope"],
          CapCostYInt = model_inputs["CapCostYInt"],
          CapCostX = model_inputs["CapCostX"],
          ProdIncentRate = model_inputs["ProdIncentRate"],
          MaxProdIncent = model_inputs["MaxProdIncent"],
          MaxSizeForProdIncent = model_inputs["MaxSizeForProdIncent"],
          two_party_factor = model_inputs["two_party_factor"],
          analysis_years = model_inputs["analysis_years"],
          AnnualElecLoad = model_inputs["AnnualElecLoad"],
          LoadProfile = model_inputs["LoadProfile"],
          ProdFactor = model_inputs["ProdFactor"],
          StorageMinChargePcent = model_inputs["StorageMinChargePcent"],
          EtaStorIn = model_inputs["EtaStorIn"],
          EtaStorOut = model_inputs["EtaStorOut"],
          InitSOC = model_inputs["InitSOC"],
          MaxSize = model_inputs["MaxSize"],
          MinStorageSizeKW = model_inputs["MinStorageSizeKW"],
          MaxStorageSizeKW = model_inputs["MaxStorageSizeKW"],
          MinStorageSizeKWH = model_inputs["MinStorageSizeKWH"],
          MaxStorageSizeKWH = model_inputs["MaxStorageSizeKWH"],
          TechClassMinSize = model_inputs["TechClassMinSize"],
          MinTurndown = model_inputs["MinTurndown"],
          FuelRate = model_inputs["FuelRate"],
          FuelAvail = model_inputs["FuelAvail"],
          FixedMonthlyCharge = model_inputs["FixedMonthlyCharge"],
          AnnualMinCharge = model_inputs["AnnualMinCharge"],
          MonthlyMinCharge = model_inputs["MonthlyMinCharge"],
          ExportRates = model_inputs["ExportRates"],
          TimeStepRatchetsMonth = model_inputs["TimeStepRatchetsMonth"],
          DemandRatesMonth = model_inputs["DemandRatesMonth"],
          DemandLookbackPercent = model_inputs["DemandLookbackPercent"],
          MaxDemandInTier = model_inputs["MaxDemandInTier"],
          MaxDemandMonthsInTier = model_inputs["MaxDemandMonthsInTier"],
          MaxUsageInTier = model_inputs["MaxUsageInTier"],
          FuelBurnRateM = model_inputs["FuelBurnRateM"],
          FuelBurnRateB = model_inputs["FuelBurnRateB"],
          NMILLimits = model_inputs["NMILLimits"],
          TechToNMILMapping = model_inputs["TechToNMILMapping"],
          DemandRates = model_inputs["DemandRates"],
          TimeStepRatchets = model_inputs["TimeStepRatchets"],
          DemandLookbackMonths = model_inputs["DemandLookbackMonths"],
          CapCostSegCount = model_inputs["CapCostSegCount"],
          FuelBinCount = model_inputs["FuelBinCount"],
          DemandBinCount  = model_inputs["DemandBinCount"],
          DemandMonthsBinCount = model_inputs["DemandMonthsBinCount"],
          TimeStepCount = model_inputs["TimeStepCount"],
          NumRatchets = model_inputs["NumRatchets"],
          TimeStepScaling = model_inputs["TimeStepScaling"],
          OMcostPerUnitProd = model_inputs["OMcostPerUnitProd"])

    MAXTIME = data["inputs"]["Scenario"]["timeout_seconds"]

    return reopt_run(reo_model, MAXTIME, p)
end

function reopt_run(reo_model, MAXTIME::Int64, p::Parameter)

	REopt = reo_model
    Obj = 2  # 1 for minimize LCC, 2 for min LCC AND high mean SOC

    @variables REopt begin
        binNMLorIL[p.NMILRegime], Bin
        binSegChosen[p.Tech, p.Seg], Bin
        dvSystemSize[p.Tech, p.Seg] >= 0
        dvGrid[p.Load, p.TimeStep, p.DemandBin, p.FuelBin, p.DemandMonthsBin] >= 0
        dvRatedProd[p.Tech, p.Load, p.TimeStep, p.Seg, p.FuelBin] >= 0
        dvProdIncent[p.Tech] >= 0
        binProdIncent[p.Tech], Bin
        binSingleBasicTech[p.Tech,p.TechClass], Bin
        dvPeakDemandE[p.Ratchets, p.DemandBin] >= 0
        dvPeakDemandEMonth[p.Month, p.DemandMonthsBin] >= 0
        dvElecToStor[p.Tech, p.TimeStep] >= 0
        dvElecFromStor[p.TimeStep] >= 0
        dvStoredEnergy[p.TimeStepBat] >= 0
        dvStorageSizeKWH >= 0
        dvStorageSizeKW >= 0
        dvMeanSOC >= 0
        dvFuelCost[p.Tech, p.FuelBin]
        dvFuelUsed[p.Tech, p.FuelBin]
        binTechIsOnInTS[p.Tech, p.TimeStep], Bin
        MinChargeAdder >= 0
        binDemandTier[p.Ratchets, p.DemandBin], Bin
        binDemandMonthsTier[p.Month, p.DemandMonthsBin], Bin
        binUsageTier[p.Month, p.FuelBin], Bin
        dvPeakDemandELookback >= 0

    # ADDED due to implied types
        ElecToBatt[p.Tech] >= 0
        UsageInTier[p.Month, p.FuelBin] >= 0
        #TotalTechCapCosts >= 0
        #TotalStorageCapCosts >= 0
        #TotalEnergyCharges >= 0
        #DemandTOUCharges >= 0
        #DemandFlatCharges >= 0
        #TotalDemandCharges >= 0
        #TotalFixedCharges >= 0
        #TotalEnergyExports <= 0
        #TotalProductionIncentive >= 0
        #TotalMinCharge >= 0

    # ADDED due to calculations
        #Year1ElecProd
        #AverageElecProd
        #Year1WindProd
        #AverageWindProd

    # ADDED for modeling
        binMinTurndown[p.Tech, p.TimeStep], Bin
    end

    ##############################################################################
	#############  		Constraints									 #############
	##############################################################################
    #TODO: account for exist formatting
    for t in p.Tech
        if p.MaxSize == 0
            for LD in p.Load
                if p.TechToLoadMatrix[t,LD] == 0
                    for ts in p.TimeStep
                        if p.LoadProfile[LD,ts] == 0
                            for s in p.Seg, fb in p.FuelBin
                            @constraint(REopt, dvRatedProd[t, LD, ts, s, fb] == 0)
                            end
                        end
                    end
                end
            end
        end
    end

    ## Fuel Tracking Constraints
    @constraint(REopt, [t in p.Tech, fb in p.FuelBin],
                sum(p.ProdFactor[t,LD,ts] * p.LevelizationFactor[t] * dvRatedProd[t,LD,ts,s,fb] * p.FuelBurnRateM[t,LD,fb] * p.TimeStepScaling
                    for ts in p.TimeStep, LD in p.Load, s in p.Seg) +
                sum(binTechIsOnInTS[t,ts] * p.FuelBurnRateB[t,LD,fb] * p.TimeStepScaling
                    for ts in p.TimeStep, LD in p.Load) == dvFuelUsed[t,fb])
    @constraint(REopt, FuelCont[t in p.Tech, fb in p.FuelBin],
                dvFuelUsed[t,fb] <= p.FuelAvail[t,fb])
    @constraint(REopt, [t in p.Tech, fb in p.FuelBin],
                sum(p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t] * dvRatedProd[t,LD,ts,s,fb] * p.FuelBurnRateM[t,LD,fb] * p.TimeStepScaling * p.FuelRate[t,fb,ts] * p.pwf_e
                    for ts in p.TimeStep, LD in p.Load, s in p.Seg) +
                sum(binTechIsOnInTS[t,ts] * p.FuelBurnRateB[t,LD,fb] * p.TimeStepScaling * p.FuelRate[t,fb,ts] * p.pwf_e
                    for ts in p.TimeStep, LD in p.Load) == dvFuelCost[t,fb])

    ### Switch Constraints
    @constraint(REopt, [t in p.Tech, ts in p.TimeStep],
                sum(p.ProdFactor[t,LD,ts] * dvRatedProd[t,LD,ts,s,fb] for LD in p.Load, s in p.Seg, fb in p.FuelBin) <=
                p.MaxSize[t] * 100 * binTechIsOnInTS[t,ts])
    @constraint(REopt, [t in p.Tech, ts in p.TimeStep],
                sum(p.MinTurndown[t] * dvSystemSize[t,s] for s in p.Seg) -
      		    sum(dvRatedProd[t,LD,ts,s,fb] for LD in p.Load, s in p.Seg, fb in p.FuelBin) <=
                p.MaxSize[t] * (1 - binTechIsOnInTS[t, ts]))

    ### Boundary Conditions and Size Limits
    @constraint(REopt, dvStoredEnergy[0] == p.InitSOC * dvStorageSizeKWH / p.TimeStepScaling)
    @constraint(REopt, p.MinStorageSizeKWH <= dvStorageSizeKWH)
    @constraint(REopt, dvStorageSizeKWH <=  p.MaxStorageSizeKWH)
    @constraint(REopt, p.MinStorageSizeKW <= dvStorageSizeKW)
    @constraint(REopt, dvStorageSizeKW <=  p.MaxStorageSizeKW)

    ### Battery Operations
    # AC electricity to be stored, generated by each Tech, applied to the storage load, for each timestep
    @constraint(REopt, [t in p.Tech, ts in p.TimeStep],
    	        dvElecToStor[t,ts] == sum(p.ProdFactor[t,"1S",ts] * p.LevelizationFactor[t] * dvRatedProd[t,"1S",ts,s,fb]
                                        for s in p.Seg, fb in p.FuelBin))
    @constraint(REopt, [ts in p.TimeStep],
    	        dvStoredEnergy[ts] == dvStoredEnergy[ts-1] + sum(dvElecToStor[t,ts] * p.EtaStorIn[t,"1S"] for t in p.Tech) - dvElecFromStor[ts] / p.EtaStorOut["1S"])
    @constraint(REopt, [ts in p.TimeStep],
    	        dvElecFromStor[ts] / p.EtaStorOut["1S"] <=  dvStoredEnergy[ts-1])
    @constraint(REopt, [ts in p.TimeStep],
    	        dvStoredEnergy[ts] >=  p.StorageMinChargePcent * dvStorageSizeKWH / p.TimeStepScaling)

    ### Operational
    # Storage inverter is AC rated. Following constrains the energy / timestep throughput of the inverter
    #     to the sum of the energy in and energy out of the battery.
    @constraint(REopt, [ts in p.TimeStep],
    	        dvStorageSizeKW * p.TimeStepScaling >=  dvElecFromStor[ts] + sum(dvElecToStor[t,ts] for t in p.Tech))
    @constraint(REopt, dvMeanSOC == sum(dvStoredEnergy[ts] for ts in p.TimeStep) / p.TimeStepCount)
    @constraint(REopt, [ts in p.TimeStep],
    	        dvStorageSizeKWH >=  dvStoredEnergy[ts] * p.TimeStepScaling)
    @constraint(REopt, [t in p.Tech],
    	        ElecToBatt[t] == sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t,LD,ts] * p.LevelizationFactor[t]
                                     for ts in p.TimeStep, LD in ["1S"], s in p.Seg, fb in p.FuelBin))

    ### Binary Bookkeeping
    @constraint(REopt, [t in p.Tech],
                sum(binSegChosen[t,s] for s in p.Seg) == 1)
    @constraint(REopt, [b in p.TechClass],
                sum(binSingleBasicTech[t,b] for t in p.Tech) <= 1)

    ### Capital Cost Constraints (NEED Had to modify b/c AxisArrays doesn't do well with range 0:20
    @constraint(REopt, [t in p.Tech, s in p.Seg],
                dvSystemSize[t,s] <= p.CapCostX[t,s+1] * binSegChosen[t,s])
    @constraint(REopt, [t in p.Tech, s in p.Seg],
                dvSystemSize[t,s] >= p.CapCostX[t,s] * binSegChosen[t,s])

    ### Production Incentive Cap Module
    @constraint(REopt, [t in p.Tech],
                dvProdIncent[t] <= binProdIncent[t] * p.MaxProdIncent[t] * p.pwf_prod_incent[t])
    @constraint(REopt, [t in p.Tech],
                dvProdIncent[t] <= sum(p.ProdFactor[t, LD, ts] * p.LevelizationFactorProdIncent[t] * dvRatedProd[t,LD,ts,s,fb] *
                                       p.TimeStepScaling * p.ProdIncentRate[t, LD] * p.pwf_prod_incent[t]
                                       for LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))
    @constraint(REopt, [t in p.Tech, LD in p.Load,ts in p.TimeStep],
                sum(dvSystemSize[t,s] for s in p.Seg) <= p.MaxSizeForProdIncent[t] + p.MaxSize[t] * (1 - binProdIncent[t]))

    ### System Size and Production Constraints
    @constraint(REopt, [t in p.Tech, s in p.Seg],
                dvSystemSize[t,s] <=  p.MaxSize[t])
    @constraint(REopt, [b in p.TechClass],
                sum(dvSystemSize[t, s] * p.TechToTechClassMatrix[t,b] for t in p.Tech, s in p.Seg) >= p.TechClassMinSize[b])
    #NEED to tighten bound and check logic

    for t in p.Tech
        if p.MinTurndown[t] > 0
            @constraint(REopt, [LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin],
                        binMinTurndown[t,ts] * p.MinTurndown[t] <= dvRatedProd[t,LD,ts,s,fb])
            @constraint(REopt, [LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin],
                        dvRatedProd[t,LD,ts,s,fb] <= binMinTurndown[t,ts] * p.AnnualElecLoad)
        end
    end

    @constraint(REopt, [t in p.Tech, s in p.Seg, ts in p.TimeStep],
                sum(dvRatedProd[t,LD,ts,s,fb] for LD in p.Load, fb in p.FuelBin) <= dvSystemSize[t, s])

    for LD in p.Load
        if LD != "1R" && LD != "1S"
            @constraint(REopt, [ts in p.TimeStep],
                        sum(p.ProdFactor[t,LD,ts] * p.LevelizationFactor[t] * dvRatedProd[t,LD,ts,s,fb]
                            for t in p.Tech, s in p.Seg, fb in p.FuelBin) <= p.LoadProfile[LD,ts])
        end
    end

    @constraint(REopt, [LD in ["1R"], ts in p.TimeStep],
                sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t,LD,ts] * p.LevelizationFactor[t]
                    for t in p.Tech, s in p.Seg, fb in p.FuelBin) + dvElecFromStor[ts] >= p.LoadProfile[LD,ts])

    ###  Net Meter Module
    @constraint(REopt, sum(binNMLorIL[n] for n in p.NMILRegime) == 1)

    @constraint(REopt, [n in p.NMILRegime],
                sum(p.TechToNMILMapping[t,n] * p.TurbineDerate[t] * dvSystemSize[t,s]
                    for t in p.Tech, s in p.Seg) <= p.NMILLimits[n] * binNMLorIL[n])

    ###  Rate Variable Definitions
    @constraint(REopt, [t in ["UTIL1"], LD in p.Load, fb in p.FuelBin, ts in p.TimeStep],
    	        sum(dvRatedProd[t,LD,ts,s,fb] for s in p.Seg) == sum(dvGrid[LD,ts,db,fb,dbm] for db in p.DemandBin, dbm in p.DemandMonthsBin))
    @constraint(REopt, [t in ["UTIL1"], fb in p.FuelBin, m in p.Month],
    	        UsageInTier[m, fb] ==  sum(dvRatedProd[t,LD,ts,s,fb] * p.TimeStepScaling for LD in p.Load, ts in p.TimeStepRatchetsMonth[m], s in p.Seg))

    ### Fuel Bins
    @constraint(REopt, [fb in p.FuelBin, m in p.Month],
                UsageInTier[m, fb] <= binUsageTier[m, fb] * p.MaxUsageInTier[fb])
    @constraint(REopt, [fb in 2:p.FuelBinCount, m in p.Month],
    	        binUsageTier[m, fb] - binUsageTier[m, fb-1] <= 0)
    @constraint(REopt, [fb in 2:p.FuelBinCount, m in p.Month],
    	        binUsageTier[m, fb] * p.MaxUsageInTier[fb-1] - UsageInTier[m, fb-1] <= 0)

    ### Peak Demand Energy Rates
    for db in p.DemandBin
        @constraint(REopt, [r in p.Ratchets],
                    dvPeakDemandE[r, db] <= binDemandTier[r,db] * p.MaxDemandInTier[db])
        if db >= 2
            @constraint(REopt, [r in p.Ratchets],
                        binDemandTier[r, db] - binDemandTier[r, db-1] <= 0)
            @constraint(REopt, [r in p.Ratchets],
                        binDemandTier[r, db] * p.MaxDemandInTier[db-1] - dvPeakDemandE[r, db-1] <= 0)
        end
    end

    if !isempty(p.TimeStepRatchets)
        @constraint(REopt, [db in p.DemandBin, r in p.Ratchets, ts in p.TimeStepRatchets[r]],
                    dvPeakDemandE[r,db] >= sum(dvGrid[LD,ts,db,fb,dbm] for LD in p.Load, fb in p.FuelBin, dbm in p.DemandMonthsBin))
        @constraint(REopt, [db in p.DemandBin, r in p.Ratchets, ts in p.TimeStepRatchets[r]],
                    dvPeakDemandE[r,db] >= p.DemandLookbackPercent * dvPeakDemandELookback)
    end

    ### Peak Demand Energy Rachets
    for dbm in p.DemandMonthsBin
        @constraint(REopt, [m in p.Month],
        	        dvPeakDemandEMonth[m, dbm] <= binDemandMonthsTier[m,dbm] * p.MaxDemandMonthsInTier[dbm])
        if dbm >= 2
        @constraint(REopt, [m in p.Month],
        	        binDemandMonthsTier[m, dbm] - binDemandMonthsTier[m, dbm-1] <= 0)
        @constraint(REopt, [m in p.Month],
        	        binDemandMonthsTier[m, dbm] * p.MaxDemandMonthsInTier[dbm-1] <= dvPeakDemandEMonth[m, dbm-1])
        end
    end
    @constraint(REopt, [dbm in p.DemandMonthsBin, m in p.Month, ts in p.TimeStepRatchetsMonth[m]],
    	        dvPeakDemandEMonth[m, dbm] >=  sum(dvGrid[LD,ts,db,fb,dbm] for LD in p.Load, db in p.DemandBin, fb in p.FuelBin))
    @constraint(REopt, [LD in p.Load, lbm in p.DemandLookbackMonths],
                dvPeakDemandELookback >= sum(dvPeakDemandEMonth[lbm, dbm] for dbm in p.DemandMonthsBin))

    #= Annual energy produced by Tech's' (except UTIL) used to meet the 1R, 1W, and 1S loads must
    be less than the annual site load (in kWh) =#
    TechIsNotGridSet = filter(t->p.TechIsGrid[t] == 0, p.Tech)
    @constraint(REopt, sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t] *  p.TimeStepScaling
                           for t in TechIsNotGridSet, LD in ["1R", "1W", "1S"],
                           ts in p.TimeStep, s in p.Seg, fb in p.FuelBin) <=  p.AnnualElecLoad)

    ###
    #Added, but has awful bounds
    @constraint(REopt, [t in p.Tech, b in p.TechClass],
                sum(dvSystemSize[t, s] * p.TechToTechClassMatrix[t, b] for s in p.Seg) <= p.MaxSize[t] * binSingleBasicTech[t, b])

    for t in p.Tech
        if p.TechToTechClassMatrix[t, "PV"] == 1 || p.TechToTechClassMatrix[t, "WIND"] == 1
            @constraint(REopt, [ts in p.TimeStep, s in p.Seg],
            	        sum(dvRatedProd[t,LD,ts,s,fb] for fb in p.FuelBin,
                            LD in ["1R", "1W", "1X", "1S"]) ==  dvSystemSize[t, s])
        end
    end

    ### Parts of Objective

    ProdLoad = ["1R", "1W", "1X", "1S"]

    PVTechs = filter(t->p.TechToTechClassMatrix[t, "PV"] == 1, p.Tech)
    @expression(REopt, Year1PvProd, sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.TimeStepScaling
                                            for t in PVTechs, s in p.Seg, fb in p.FuelBin, ts in p.TimeStep, LD in ProdLoad))
    @expression(REopt, AveragePvProd, sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.TimeStepScaling * p.LevelizationFactor[t]
                                              for t in PVTechs, s in p.Seg, fb in p.FuelBin, ts in p.TimeStep, LD in ProdLoad))

    WindTechs = filter(t->p.TechToTechClassMatrix[t, "WIND"] == 1, p.Tech)
    @expression(REopt, Year1WindProd, sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.TimeStepScaling
                                            for t in WindTechs, s in p.Seg, fb in p.FuelBin, ts in p.TimeStep, LD in ProdLoad))
    @expression(REopt, AverageWindProd, sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.TimeStepScaling * p.LevelizationFactor[t]
                                              for t in WindTechs, s in p.Seg, fb in p.FuelBin, ts in p.TimeStep, LD in ProdLoad))

    GeneratorTechs = filter(t->p.TechToTechClassMatrix[t, "GENERATOR"] == 1, p.Tech)
    @expression(REopt, Year1GenProd, sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.TimeStepScaling
                       for t in GeneratorTechs, s in p.Seg, fb in p.FuelBin, ts in p.TimeStep, LD in ProdLoad))
    @expression(REopt, AverageGenProd, sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.TimeStepScaling * p.LevelizationFactor[t]
                       for t in GeneratorTechs, s in p.Seg, fb in p.FuelBin, ts in p.TimeStep, LD in ProdLoad))

    @expression(REopt, TotalTechCapCosts, sum(p.CapCostSlope[t, s] * dvSystemSize[t, s] + p.CapCostYInt[t,s] * binSegChosen[t,s]
                                                for t in p.Tech, s in p.Seg))
    @expression(REopt, TotalStorageCapCosts, dvStorageSizeKWH * p.StorageCostPerKWH + dvStorageSizeKW * p.StorageCostPerKW)
    @expression(REopt, TotalPerUnitSizeOMCosts, sum(p.OMperUnitSize[t] * p.pwf_om * dvSystemSize[t, s] for t in p.Tech, s in p.Seg))

    if !isempty(GeneratorTechs)
        @expression(REopt, TotalPerUnitProdOMCosts, sum(dvRatedProd[t,LD,ts,s,fb] * p.TimeStepScaling * p.ProdFactor[t,LD,ts] * p.OMcostPerUnitProd[t] * p.pwf_om
                                                    for t in GeneratorTechs, LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))
    else
        @expression(REopt, TotalPerUnitProdOMCosts, 0.0)
    end

    @expression(REopt, GenPerUnitSizeOMCosts, sum(p.OMperUnitSize[t] * p.pwf_om * dvSystemSize[t, s] for t in GeneratorTechs, s in p.Seg))
    @expression(REopt, GenPerUnitProdOMCosts, sum(dvRatedProd[t,LD,ts,s,fb] * p.TimeStepScaling * p.ProdFactor[t,LD,ts] * p.OMcostPerUnitProd[t] * p.pwf_om
                                              for t in GeneratorTechs, LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))

    ### Aggregates of definitions
    @expression(REopt, TotalEnergyCharges, sum(dvFuelCost[t,fb]
                                                 for t in p.Tech, fb in p.FuelBin))
    @expression(REopt, TotalEnergyChargesUtil, sum(dvFuelCost["UTIL1",fb] for fb in p.FuelBin))
    @expression(REopt, TotalGenFuelCharges, sum(dvFuelCost[t,fb] for t in GeneratorTechs, fb in p.FuelBin))
    @expression(REopt, DemandTOUCharges, sum(dvPeakDemandE[r, db] * p.DemandRates[r,db] * p.pwf_e
                                               for r in p.Ratchets, db in p.DemandBin))
    @expression(REopt, DemandFlatCharges, sum(dvPeakDemandEMonth[m, dbm] * p.DemandRatesMonth[m, dbm] * p.pwf_e
                                                for m in p.Month, dbm in p.DemandMonthsBin))
    @expression(REopt, TotalDemandCharges, DemandTOUCharges + DemandFlatCharges)
    @expression(REopt, TotalFixedCharges, p.FixedMonthlyCharge * 12 * p.pwf_e)

    ### Utility and Taxable Costs
    @expression(REopt, TotalEnergyExports, sum(dvRatedProd[t,LD,ts,s,fb] * p.TimeStepScaling * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t] *
                                                 p.ExportRates[t,LD,ts] * p.pwf_e for t in p.Tech, LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))
    @expression(REopt, TotalProductionIncentive, sum(dvProdIncent[t] for t in p.Tech))

    ### MinChargeAdder
    if p.AnnualMinCharge > 12 * p.MonthlyMinCharge
        TotalMinCharge = p.AnnualMinCharge * p.pwf_e
    else
        TotalMinCharge = 12 * p.MonthlyMinCharge * p.pwf_e
    end

    if TotalMinCharge > 0
        @constraint(REopt, MinChargeAdder >= TotalMinCharge - (TotalEnergyCharges + TotalDemandCharges + TotalEnergyExports + TotalFixedCharges))
	else
		@constraint(REopt, MinChargeAdder == 0)
    end
    #= Note: 0.999*MinChargeAdder in Obj b/c when TotalMinCharge > (TotalEnergyCharges + TotalDemandCharges + TotalEnergyExports + TotalFixedCharges)
		it is arbitrary where the min charge ends up (eg. could be in TotalDemandCharges or MinChargeAdder).
		0.001*MinChargeAdder is added back into LCC when writing to results.  =#

    # Define Rates
    r_tax_fraction_owner = (1 - p.r_tax_owner)
    r_tax_fraction_offtaker = (1 - p.r_tax_offtaker)

    ### Objective Function
    @expression(REopt, REcosts,
		# Capital Costs
		TotalTechCapCosts + TotalStorageCapCosts +

		# Fixed O&M, tax deductible for owner
		TotalPerUnitSizeOMCosts * r_tax_fraction_owner +

		# Variable O&M, tax deductible for owner
		TotalPerUnitProdOMCosts * r_tax_fraction_owner +

		# Utility Bill, tax deductible for offtaker
		(TotalEnergyCharges + TotalDemandCharges + TotalEnergyExports + TotalFixedCharges + 0.999*MinChargeAdder) * r_tax_fraction_offtaker -

		# Subtract Incentives, which are taxable
		TotalProductionIncentive * r_tax_fraction_owner
	)

    if Obj == 1
		@objective(REopt, Min, REcosts)
	elseif Obj == 2  # Keep SOC high
		@objective(REopt, Min, REcosts - dvMeanSOC)
	end
	optimize!(REopt)

    ##############################################################################
    #############  		Outputs    									 #############
    ##############################################################################
    @expression(REopt, ExportedElecPV,
                sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t] * p.TimeStepScaling
                    for t in PVTechs, LD in ["1W", "1X"], ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))
    @expression(REopt, ExportedElecWIND,
                sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t] * p.TimeStepScaling
                    for t in WindTechs, LD in ["1W", "1X"], ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))
    @expression(REopt, ExportedElecGEN,
                sum(dvRatedProd[t,LD,ts,s,fb] * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t] * p.TimeStepScaling
                    for t in GeneratorTechs, LD in ["1W", "1X"], ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))

    # Needs levelization factor?
    @expression(REopt, ExportBenefitYr1,
                sum(dvRatedProd[t,LD,ts,s,fb] * p.TimeStepScaling * p.ProdFactor[t, LD, ts] * p.ExportRates[t,LD,ts] * p.LevelizationFactor[t]
                    for t in p.Tech, LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))
    @expression(REopt, Year1UtilityEnergy,
                sum(dvRatedProd["UTIL1",LD,ts,s,fb] * p.TimeStepScaling * p.ProdFactor["UTIL1", LD, ts]
                    for LD in p.Load, ts in p.TimeStep, s in p.Seg, fb in p.FuelBin))

    ojv = round(JuMP.objective_value(REopt)+ 0.001*value(MinChargeAdder))
    Year1EnergyCost = TotalEnergyChargesUtil / p.pwf_e
    Year1DemandCost = TotalDemandCharges / p.pwf_e
    Year1DemandTOUCost = DemandTOUCharges / p.pwf_e
    Year1DemandFlatCost = DemandFlatCharges / p.pwf_e
    Year1FixedCharges = TotalFixedCharges / p.pwf_e
    Year1MinCharges = MinChargeAdder / p.pwf_e
    Year1Bill = Year1EnergyCost + Year1DemandCost + Year1FixedCharges + Year1MinCharges

    results = Dict{String, Any}("lcc" => ojv)

    results["batt_kwh"] = value(dvStorageSizeKWH)
    results["batt_kw"] = value(dvStorageSizeKW)

    if results["batt_kwh"] != 0
    	@expression(REopt, soc[ts in p.TimeStep], dvStoredEnergy[ts] / results["batt_kwh"])
        results["year_one_soc_series_pct"] = value.(soc)
    else
        results["year_one_soc_series_pct"] = []
    end

    if !isempty(PVTechs)
        results["PV"] = Dict()
        results["pv_kw"] = round(value(sum(dvSystemSize[t,s] for s in p.Seg, t in PVTechs)), digits=4)
        @expression(REopt, PVtoBatt[t in PVTechs, ts in p.TimeStep],
                    sum(dvRatedProd[t, "1S", ts, s, fb] * p.ProdFactor[t, "1S", ts] * p.LevelizationFactor[t] for s in p.Seg, fb in p.FuelBin))
    end

    WindTechs = filter(t->p.TechToTechClassMatrix[t, "WIND"] == 1, p.Tech)
    if !isempty(WindTechs)
        results["Wind"] = Dict()
        results["wind_kw"] = round(value(sum(dvSystemSize[t,s] for s in p.Seg, t in WindTechs)), digits=4)
        @expression(REopt, WINDtoBatt[t in WindTechs, ts in p.TimeStep],
                    sum(dvRatedProd[t, "1S", ts, s, fb] * p.ProdFactor[t, "1S", ts] * p.LevelizationFactor[t] for s in p.Seg, fb in p.FuelBin))
    end

	results["gen_net_fixed_om_costs"] = 0
	results["gen_net_variable_om_costs"] = 0
	results["gen_total_fuel_cost"] = 0
	results["gen_year_one_fuel_cost"] = 0
	results["gen_year_one_variable_om_costs"] = 0

    GeneratorTechs = filter(t->p.TechToTechClassMatrix[t, "GENERATOR"] == 1, p.Tech)
    if !isempty(GeneratorTechs)
    	if value(sum(dvSystemSize[t,s] for s in p.Seg, t in GeneratorTechs)) > 0
			results["Generator"] = Dict()
            results["generator_kw"] = value(sum(dvSystemSize[t,s] for s in p.Seg, t in GeneratorTechs))
			results["gen_net_fixed_om_costs"] = round(value(GenPerUnitSizeOMCosts) * r_tax_fraction_owner, digits=0)
			results["gen_net_variable_om_costs"] = round(value(GenPerUnitProdOMCosts) * r_tax_fraction_owner, digits=0)
	        results["gen_total_fuel_cost"] = round(value(TotalGenFuelCharges) * r_tax_fraction_offtaker, digits=2)
	        results["gen_year_one_fuel_cost"] = round(value(TotalGenFuelCharges) * r_tax_fraction_offtaker / p.pwf_e, digits=2)
	        results["gen_year_one_variable_om_costs"] = round(value(GenPerUnitProdOMCosts) * r_tax_fraction_owner / p.pwf_om, digits=0)
		end
    end

    net_capital_costs_plus_om = value(TotalTechCapCosts + TotalStorageCapCosts) +
                                value(TotalPerUnitSizeOMCosts + TotalPerUnitProdOMCosts) * r_tax_fraction_owner +
                                value(TotalGenFuelCharges) * r_tax_fraction_offtaker

    push!(results, Dict("year_one_utility_kwh" => round(value(Year1UtilityEnergy), digits=2),
						 "year_one_energy_cost" => round(value(Year1EnergyCost), digits=2),
						 "year_one_demand_cost" => round(value(Year1DemandCost), digits=2),
						 "year_one_demand_tou_cost" => round(value(Year1DemandTOUCost), digits=2),
						 "year_one_demand_flat_cost" => round(value(Year1DemandFlatCost), digits=2),
						 "year_one_export_benefit" => round(value(ExportBenefitYr1), digits=0),
						 "year_one_fixed_cost" => round(value(Year1FixedCharges), digits=0),
						 "year_one_min_charge_adder" => round(value(Year1MinCharges), digits=2),
						 "year_one_bill" => round(value(Year1Bill), digits=2),
						 "year_one_payments_to_third_party_owner" => round(value(TotalDemandCharges) / p.pwf_e, digits=0),
						 "total_energy_cost" => round(value(TotalEnergyCharges) * r_tax_fraction_offtaker, digits=2),
						 "total_demand_cost" => round(value(TotalDemandCharges) * r_tax_fraction_offtaker, digits=2),
						 "total_fixed_cost" => round(value(TotalFixedCharges) * r_tax_fraction_offtaker, digits=2),
						 "total_export_benefit" => round(value(TotalEnergyExports) * r_tax_fraction_offtaker, digits=2),
						 "total_min_charge_adder" => round(value(MinChargeAdder) * r_tax_fraction_offtaker, digits=2),
						 "total_payments_to_third_party_owner" => 0,
						 "net_capital_costs_plus_om" => round(net_capital_costs_plus_om, digits=0),
						 "year_one_energy_produced" => round(value(Year1PvProd), digits=0),
                         "average_yearly_pv_energy_produced" => round(value(AveragePvProd), digits=0),
						 "average_annual_energy_exported_pv" => round(value(ExportedElecPV), digits=0),
						 "year_one_wind_energy_produced" => round(value(Year1WindProd), digits=0),
						 "average_wind_energy_produced" => round(value(AverageWindProd), digits=0),
						 "average_annual_energy_exported_wind" => round(value(ExportedElecWIND), digits=0),
                         "year_one_gen_energy_produced" => round(value(Year1GenProd), digits=0),
                         "average_yearly_gen_energy_produced" => round(value(AverageGenProd), digits=0),
                         "average_annual_energy_exported_gen" => round(value(ExportedElecGEN), digits=0),
						 "net_capital_costs" => round(value(TotalTechCapCosts + TotalStorageCapCosts), digits=2))...)

    @expression(REopt, GeneratorFuelUsed, sum(dvFuelUsed[t, fb] for t in GeneratorTechs, fb in p.FuelBin))
    results["fuel_used_gal"] = round(value(GeneratorFuelUsed), digits=2)

    @expression(REopt, GridToBatt[ts in p.TimeStep],
                sum(dvRatedProd["UTIL1", "1S", ts, s, fb] * p.ProdFactor["UTIL1", "1S", ts] * p.LevelizationFactor["UTIL1"]
					for s in p.Seg, fb in p.FuelBin))
    results["GridToBatt"] = round.(value.(GridToBatt), digits=3)

    @expression(REopt, GridToLoad[ts in p.TimeStep],
                sum(dvRatedProd["UTIL1", "1R", ts, s, fb] * p.ProdFactor["UTIL1", "1R", ts] * p.LevelizationFactor["UTIL1"]
					for s in p.Seg, fb in p.FuelBin))
    results["GridToLoad"] = round.(value.(GridToLoad), digits=3)

	if !isempty(GeneratorTechs)
		@expression(REopt, GENERATORtoBatt[ts in p.TimeStep],
					sum(dvRatedProd[t, "1S", ts, s, fb] * p.ProdFactor[t, "1S", ts] * p.LevelizationFactor[t]
						for t in GeneratorTechs, s in p.Seg, fb in p.FuelBin))
    	results["GENERATORtoBatt"] = round.(value.(GENERATORtoBatt), digits=3)

		@expression(REopt, GENERATORtoGrid[ts in p.TimeStep],
					sum(dvRatedProd[t, LD, ts, s, fb] * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t]
						for t in GeneratorTechs, LD in ["1W", "1X"], s in p.Seg, fb in p.FuelBin))
		results["GENERATORtoGrid"] = round.(value.(GENERATORtoGrid), digits=3)

		@expression(REopt, GENERATORtoLoad[ts in p.TimeStep],
					sum(dvRatedProd[t, "1R", ts, s, fb] * p.ProdFactor[t, "1R", ts] * p.LevelizationFactor[t]
						for t in GeneratorTechs, s in p.Seg, fb in p.FuelBin))
		results["GENERATORtoLoad"] = round.(value.(GENERATORtoLoad), digits=3)
    else
    	results["GENERATORtoBatt"] = []
		results["GENERATORtoGrid"] = []
		results["GENERATORtoLoad"] = []
	end

	if !isempty(PVTechs)
		@expression(REopt, PVtoLoad[ts in p.TimeStep],
					sum(dvRatedProd[t, "1R", ts, s, fb] * p.ProdFactor[t, "1R", ts] * p.LevelizationFactor[t]
						for t in PVTechs, s in p.Seg, fb in p.FuelBin))
    	results["PVtoLoad"] = round.(value.(PVtoLoad), digits=3)

		@expression(REopt, PVtoGrid[ts in p.TimeStep],
					sum(dvRatedProd[t, LD, ts, s, fb] * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t]
						for t in PVTechs, LD in ["1W", "1X"], s in p.Seg, fb in p.FuelBin))
    	results["PVtoGrid"] = round.(value.(PVtoGrid), digits=3)

		@expression(REopt, PVPerUnitSizeOMCosts,
					sum(p.OMperUnitSize[t] * p.pwf_om * dvSystemSize[t, s] for t in PVTechs, s in p.Seg))
		results["pv_net_fixed_om_costs"] = round(value(PVPerUnitSizeOMCosts) * r_tax_fraction_owner, digits=0)
	else
		results["PVtoLoad"] = []
		results["PVtoGrid"] = []
		results["pv_net_fixed_om_costs"] = 0
	end

	if !isempty(WindTechs)
		@expression(REopt, WINDtoLoad[ts in p.TimeStep],
					sum(dvRatedProd[t, "1R", ts, s, fb] * p.ProdFactor[t, "1R", ts] * p.LevelizationFactor[t]
						for t in WindTechs, s in p.Seg, fb in p.FuelBin))
		results["WINDtoLoad"] = round.(value.(WINDtoLoad), digits=3)

		@expression(REopt, WINDtoGrid[ts in p.TimeStep],
					sum(dvRatedProd[t, LD, ts, s, fb] * p.ProdFactor[t, LD, ts] * p.LevelizationFactor[t]
						for t in WindTechs, s in p.Seg, fb in p.FuelBin, LD in ["1W", "1X"]))
		results["WINDtoGrid"] = round.(value.(WINDtoGrid), digits=3)
	else
		results["WINDtoLoad"] = []
    	results["WINDtoGrid"] = []
	end

    results["Load"] = p.LoadProfile["1R", :]

	if termination_status(REopt) == MOI.TIME_LIMIT
		status = "timed-out"
    elseif termination_status(REopt) == MOI.OPTIMAL
        status = "optimal"
    else
        status = "not optimal"
    end

    results["status"] = status
    return results
end


text:
using GLAbstraction, GLPlot, Reactive

window = createdisplay(h=1000,w=1500, eyeposition=Vec3(2,0,0))
text = "whatup internet!?\n#Crusont"

glplot(text, scale=Vec2(1/50), color=Vec4[Vec4(rand(), rand(), 0,1) for i=1:length(text)]) # You can either, supply a texture with colors
#glplot(text, scale=Vec2(1/50), color=Vec4(0,1,0,1)) # or just supply one color
#full api:
#=
gplot(text::AbstractString;
					start=Vec3(0), scale=Vec2(1/500), color=Vec4(0,0,1,1), backgroundcolor=Vec4(0),
					lineheight=Vec3(0,0,0), advance=Vec3(0,0,0), rotation=Quaternion(1f0,0f0,0f0,0f0), textrotation=Quaternion(1f0,0f0,0f0,0f0),
					camera=pcamera
				)
=#
# For all attributes you can either supply values per glyph with a texture, or a scalar value for the whole text.

renderloop(window)


text:
module Unidecode
    using REPL, Unicode
    const emoji_dict = Dict((v => k) for (k, v) in REPL.REPLCompletions.emoji_symbols)
    function decode_emoji(x)
        string((c ∈ keys(emoji_dict) ? emoji_dict[c][2:end] : c for c in graphemes(x))...)
    end
    const latex_dict = Dict((v => k) for (k, v) in REPL.REPLCompletions.latex_symbols)
    function decode_latex(x)
        string((c ∈ keys(latex_dict) ? latex_dict[c][2:end] : c for c in graphemes(x))...)
    end
    function unidecode(x)
        decode_emoji(decode_latex(x))
    end
    export unidecode, decode_emoji, decode_latex
end


text:
module CalibrationTests

using CalibrationErrors
using HypothesisTests
using Parameters
using StatsBase
using StatsFuns

using LinearAlgebra
using Random

using CalibrationErrors: CalibrationErrorEstimator, SKCE, skce_kernel, skce_result_type,
    MatrixKernel

export ConsistencyTest
export DistributionFreeTest, AsymptoticLinearTest, AsymptoticQuadraticTest

# re-export
export pvalue

# defaults
abstract type CalibrationTest <: HypothesisTests.HypothesisTest end

estimator(test::CalibrationTest) = test.estimator
HypothesisTests.default_tail(test::CalibrationTest) = :right

include("consistency.jl")

include("skce/utils.jl")
include("skce/distribution_free.jl")
include("skce/asymptotic_linear.jl")
include("skce/asymptotic_quadratic.jl")

end # module


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

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
# This file was generated, do not modify it. # hide
I, D = searchbatch(G, Q, 10; parallel=true)

Threads.@threads for i in eachindex(Q)
    res, cost = search(G, Q[i], KnnResult(10))
    println(res.id) # do something with `res`
end

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
using CairoMakie, Random, Distributions
Random.seed!(13)
n = 3000
colors = to_colormap(:spring, 8)[3:end]

set_theme!(theme_dark())
fig = Figure(resolution = (600, 400))
ax = Axis(fig[1,1]; palette = (; patchcolor = colors),
    xticks = (1:7, ["cat 1", "A", "B", "C", "D", "E", "F"]),
    yticks = ([-5], ["cat 2"]), yticklabelrotation = π/2)
boxplot!(ax, fill(-5,n), rand(Normal(0, 0.5), n); orientation=:horizontal,
    whiskerwidth = 1, width = 2, color = (:orange, 0.95),
    whiskercolor = :red, mediancolor = :yellow, markersize = 8,
    strokecolor = :black, strokewidth = 1, label = "horizontal")
boxplot!(ax, fill(1,n), rand(Normal(1,  3), n); whiskerwidth = 1,
    width = 0.5, color = :dodgerblue, whiskercolor = :dodgerblue,
    mediancolor = :white, markersize = 5, strokecolor = :white,
    strokewidth = 1, label = "vertical")
for i in 2:7
    boxplot!(ax, fill(i,n), rand(Normal(rand(-2:5), 2*rand() + 0.3), n);
        whiskerwidth = 1, width = 0.35)
end
axislegend(ax, position = :lt)
save("boxplot.png", fig, px_per_unit = 2)
display(fig)


text:
# =============
# weight matrix
# =====================

using SparseArrays: AbstractSparseMatrix

import LinearAlgebra:ishermitian, issymmetric, transpose, adjoint

struct AdjacencyMatrix{G <: AbstractGraph} <: AbstractSparseMatrix{Bool, Int}
    graph::G
end

getindex(matrix::AdjacencyMatrix, s, d) = has_edge(matrix.graph, s, d)

function size(matrix::AdjacencyMatrix)
    nvg = Int(nv(matrix.graph))
    return (nvg, nvg)
end

ishermitian(::AdjacencyMatrix{<:EdgeValGraph}) = true

adjoint(matrix::AdjacencyMatrix{<:EdgeValGraph}) = matrix
transpose(matrix::AdjacencyMatrix{<:EdgeValGraph}) = transpose

struct ValMatrix{K, Tv, G <: AbstractEdgeValGraph} <:
        AbstractSparseMatrix{Tv, Int} # TODO maybe not int

    graph::G
    zero_value::Tv
end

# TODO constructor

function ValMatrix(g::AbstractEdgeValGraph{V, E_VALS}, key;
        zero_value= zero(E_VAL_for_key(E_VALS, key))) where {V, E_VALS}

    T = E_VAL_for_key(E_VALS, key)
    Z = typeof(zero_value)
    Tv = Union{T, Z}

    return ValMatrix{key, Tv, typeof(g)}(g, zero_value)
end

function ValMatrix(g::OneEdgeValGraph{V, E_VAL};
        zero_value= zero(E_VAL)) where {V, E_VAL}

    key = edgeval_keys(g)[1] # TODO implement keys

    return ValMatrix(g, key; zero_value=zero_value)
end

function size(matrix::ValMatrix)
    nvg = Int(nv(matrix.graph))
    return (nvg, nvg)
end

function getindex(matrix::ValMatrix{K, Tv}, s, d) where {K, Tv}
    return get_edgeval_or(matrix.graph, s, d, matrix.zero_value, key=K)
end

issymmetric(::ValMatrix{<:EdgeValGraph}) = true

transpose(matrix::ValMatrix{<:EdgeValGraph}) = transpose

### weights

function weights(g::AbstractEdgeValGraph; key=nokey)

    if g isa ZeroEdgeValGraph && key == nokey
        return LightGraphs.DefaultDistance(nv(g))
    elseif g isa OneEdgeValGraph && key == nokey
        return ValMatrix(g)
    end

    return ValMatrix(g, key)
end


text:
function (parareal::Parareal)(solution::PararealSolution, problem::AbstractInitialValueProblem; mode::String="SERIAL")
    cache = PararealCache(problem, parareal)
    if nprocs() == 1 || mode == "SERIAL"
        parareal_serial!(cache, solution, problem, parareal)
    elseif nprocs() > 1 && mode == "DISTRIBUTED"
        parareal_distributed!(cache, solution, problem, parareal)
    elseif nprocs() > 1 && mode == "MPI"
        parareal_mpi!(cache, solution, problem, parareal)
    else
        error("$mode is not available. Please choose between SERIAL, DISTRIBUTED, and MPI.")
    end
    return solution
end

function (parareal::Parareal)(problem::AbstractInitialValueProblem; mode::String="SERIAL", saveiterates::Bool=false)
    solution = PararealSolution(problem, parareal; saveiterates=saveiterates)
    parareal(solution, problem; mode=mode)
    return solution
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   International Geomagnetic Field Model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
#   [2] https://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f
#   [3] https://www.mathworks.com/matlabcentral/fileexchange/34388-international-geomagnetic-reference-field--igrf--model
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export igrf, igrfd

################################################################################
#                                  Functions
################################################################################

"""
    igrfd(date::Number, [r,h]::Number, λ::Number, Ω::Number, T[, P, dP]; show_warns = true)

**IGRF Model**

*Current version: v13*

Compute the geomagnetic field vector [nT] at the date `date` [Year A.D.] and
position (`r`, `λ`, `Ω`).

The position representation is defined by `T`. If `T` is `Val(:geocentric)`,
then the input must be **geocentric** coordinates:

1. Distance from the Earth center `r` [m];
1. Geocentric latitude `λ` (-90°, +90°); and
2. Geocentric longitude `Ω` (-180°, +180°).

If `T` is `Val(:geodetic)`, then the input must be **geodetic** coordinates:

1 Altitude above the reference ellipsoid `h` (WGS-84) \\[m];
2. Geodetic latitude `λ` (-90°, +90°); and
3. Geodetic longitude `Ω` (-180°, +180°).

If `T` is omitted, then it defaults to `Val(:geocentric)`.

Notice that the output vector will be represented in the same reference system
selected by the parameter `T` (geocentric or geodetic). The Y-axis of the output
reference system always points East. In case of **geocentric coordinates**, the
Z-axis points toward the center of Earth and the X-axis completes a right-handed
coordinate system. In case of **geodetic coordinates**, the X-axis is tangent to
the ellipsoid at the selected location and points toward North, whereas the
Z-axis completes a right-hand coordinate system.

The optional arguments `P` and `dP` must be two matrices with at least 14x14
real numbers. If they are present, then they will be used to store the Legendre
coefficients and their derivatives. In this case, no allocation will be
performed when computing the magnetic field. If they are not present, then 2
allocations will happen to create them.

# Keywords

* `show_warns`: Show warnings about the data (**Default** = `true`).

# Remarks

The `date` must be greater or equal to 1900 and less than or equal 2030. Notice
that a warning message is printed for dates greater than 2025.

# Disclaimer

This function is an independent implementation of the IGRF model. It contains a
more readable code than the original one in FORTRAN, because it uses features
available in Julia language.

"""
@inline igrfd(date::Number, r::Number, λ::Number, Ω::Number; show_warns = true) =
    igrfd(date, r, λ, Ω, Val(:geocentric); show_warns = show_warns)

@inline igrfd(date::Number, r::Number, λ::Number, Ω::Number,
              P::AbstractMatrix{T}, dP::AbstractMatrix{T};
              show_warns = true) where T<:Real =
    igrfd(date, r, λ, Ω, Val(:geocentric), P, dP; show_warns = show_warns)

@inline function igrfd(date::Number, rh::Number, λ::Number, Ω::Number, R::T;
                       show_warns = true) where {T<:Union{Val{:geocentric},
                                                          Val{:geodetic}}}

    P  = Matrix{Float64}(undef, 14, 14)
    dP = Matrix{Float64}(undef, 14, 14)

    return igrfd(date, rh, λ, Ω, R, P, dP; show_warns = show_warns)
end

@inline function igrfd(date::Number, rh::Number, λ::Number, Ω::Number, R::T,
                       P::AbstractMatrix{S}, dP::AbstractMatrix{S};
                       show_warns = true) where {T<:Union{Val{:geocentric},
                                                          Val{:geodetic}},
                                                 S<:Real}

    # Check if the latitude and longitude are valid.
    ( (λ < -90) || (λ > 90) ) &&
    error("The latitude must be between -90° and +90° rad.")

    ( (Ω < -180) || (Ω > 180) ) &&
    error("The longitude must be between -180° and +180° rad.")

    return igrf(date, rh, deg2rad(λ), deg2rad(Ω), R, P, dP; show_warns = show_warns)
end

"""
    igrf(date::Number, [r,h]::Number, λ::Number, Ω::Number, T[, P, dP]; show_warns = true)

**IGRF Model**

*Current version: v13*

Compute the geomagnetic field vector [nT] at the date `date` [Year A.D.] and
position (`r`, `λ`, `Ω`).

The position representation is defined by `T`. If `T` is `Val(:geocentric)`,
then the input must be **geocentric** coordinates:

1. Distance from the Earth center `r` [m];
1. Geocentric latitude `λ` (-π/2, +π/2) \\[rad]; and
2. Geocentric longitude `Ω` (-π, +π) \\[rad].

If `T` is `Val(:geodetic)`, then the input must be **geodetic** coordinates:

1 Altitude above the reference ellipsoid `h` (WGS-84) \\[m];
2. Geodetic latitude `λ` (-π/2, +π/2) \\[rad]; and
3. Geodetic longitude `Ω` (-π, +π) \\[rad].

If `T` is omitted, then it defaults to `Val(:geocentric)`.

Notice that the output vector will be represented in the same reference system
selected by the parameter `T` (geocentric or geodetic). The Y-axis of the output
reference system always points East. In case of **geocentric coordinates**, the
Z-axis points toward the center of Earth and the X-axis completes a right-handed
coordinate system. In case of **geodetic coordinates**, the X-axis is tangent to
the ellipsoid at the selected location and points toward North, whereas the
Z-axis completes a right-hand coordinate system.

The optional arguments `P` and `dP` must be two matrices with at least 14x14
real numbers. If they are present, then they will be used to store the Legendre
coefficients and their derivatives. In this case, no allocation will be
performed when computing the magnetic field. If they are not present, then 2
allocations will happen to create them.

# Keywords

* `show_warns`: Show warnings about the data (**Default** = `true`).

# Remarks

The `date` must be greater or equal to 1900 and less than or equal 2030. Notice
that a warning message is printed for dates greater than 2025.

# Disclaimer

This function is an independent implementation of the IGRF model. It contains a
more readable code than the original one in FORTRAN, because it uses features
available in Julia language.

"""
igrf(date::Number, r::Number, λ::Number, Ω::Number; show_warns = true) =
    igrf(date, r, λ, Ω, Val(:geocentric); show_warns = show_warns)

igrf(date::Number, r::Number, λ::Number, Ω::Number, P::AbstractMatrix{T},
     dP::AbstractMatrix{T}; show_warns = true) where T<:Real =
    igrf(date, r, λ, Ω, Val(:geocentric), P, dP; show_warns = show_warns)

function igrf(date::Number, r::Number, λ::Number, Ω::Number, R::T;
              show_warns::Bool = true) where {T<:Union{Val{:geocentric},
                                                       Val{:geodetic}}}

    P  = Matrix{Float64}(undef, 14, 14)
    dP = Matrix{Float64}(undef, 14, 14)

    return igrf(date, r, λ, Ω, R, P, dP; show_warns = show_warns)
end

function igrf(date::Number, r::Number, λ::Number, Ω::Number, ::Val{:geocentric},
              P::AbstractMatrix{T}, dP::AbstractMatrix{T};
              show_warns::Bool = true) where T<:Real

    # Model data
    # ==========

    max_year = 2030    # Maximum year in which we can compute.
    rel_year = 2025    # Maximum year with reliable data.
    dat_year = 2020    # Last year in which real measurements are available.

    # Input verification
    # ==================

    # Check the data, since this model is valid for years between 1900 and
    # `max_year`.
    ( (date < 1900) || (date > max_year) ) &&
    error("This IGRF version will not work for years outside the interval [1900, $max_year).")

    # Check if the latitude and longitude are valid.
    ( (λ < -pi/2) || (λ > pi/2) ) &&
    error("The latitude must be between -π/2 and +π/2 rad.")

    ( (Ω < -pi) || (Ω > pi) ) &&
    error("The longitude must be between -π and +π rad.")

    # Warn the user that for dates after the year `rel_year` the accuracy maybe
    # reduced.
    show_warns && (date > rel_year) &&
    @warn("The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than $rel_year.")

    # Input variables conversion
    # ==========================

    # Convert latitude / longitude to colatitude and east-longitude.
    θ = pi/2 - λ
    ϕ = (Ω >= 0) ? Ω : 2pi + Ω

    # The input variable `r` is in [m], but all the algorithm requires it to be
    # in [km].
    r /= 1000

    # Preliminary setup
    # =================

    # Compute the epoch that will be used to obtain the coefficients. This is
    # necessary because the IGRF provides coefficients every 5 years. Between
    # two epochs, those coefficients must be interpolated.
    idx   = floor(Int, clamp((date-1900)*0.2+1, 0, (rel_year-1900)/5))
    epoch = 1900 + (idx-1)*5

    # Compute the fraction of time from the epoch of the coefficient selected by
    # `idx`.
    Δt = date - epoch

    # Compute the maximum spherical harmonic degree for the selected date.
    n_max = (epoch < 1995) ? 10 : 13

    # Make sure we have the required amount of space in the matrices.
    (rows, cols) = size(P)
    ((rows < n_max+1) || (cols < n_max+1)) && error("Matrix `P` must have at least $(n_max+1) rows and columns.")

    (rows, cols) = size(dP)
    ((rows < n_max+1) || (cols < n_max+1)) && error("Matrix `dP` must have at least $(n_max+1) rows and columns.")

    # Compute the Schmidt quasi-normalized associated Legendre functions and
    # their first order derivative, neglecting the phase term.
    legendre!(Val(:schmidt), P, θ, false, n_max, n_max)
    dlegendre!(Val(:schmidt), dP, θ, P, false, n_max, n_max)

    # Parameters and auxiliary variables
    # ==================================

    # Auxiliary variables to select the IGRF coefficients.
    H = H_igrf
    G = G_igrf

    # Get the last index of the matrices G and H, which must be the same.
    ~, gend = size(G)
    ~, hend = size(H)

    (gend != hend) &&
    error("Internal error: the matrices G and H must have the same number of columns.")

    # Reference radius [km].
    a = 6371.2

    # Auxiliary variables to decrease the computational burden.
    sin_ϕ,  cos_ϕ  = sincos(1ϕ)
    sin_2ϕ, cos_2ϕ = sincos(2ϕ)
    ratio   = a/r
    fact    = ratio

    # Initialization of variables
    # ===========================

    dVr = 0.0   # Derivative of the Geomagnetic potential w.r.t. r.
    dVθ = 0.0   # Derivative of the Geomagnetic potential w.r.t. θ.
    dVϕ = 0.0   # Derivative of the Geomagnetic potential w.r.t. ϕ.
    ΔG  = 0.0   # Auxiliary variable to interpolate the G coefficients.
    ΔH  = 0.0   # Auxiliary variable to interpolate the H coefficients.
    kg  = 1     # Index to obtain the values of the matrix `G`.
    kh  = 1     # Index to obtain the values of the matrix `H`.

    # Geomagnetic potential
    # =====================

    @inbounds for n = 1:n_max
        aux_dVr = 0.0
        aux_dVθ = 0.0
        aux_dVϕ = 0.0

        # Compute the contributions when `m = 0`
        # ======================================

        # Get the coefficients in the epoch and interpolate to the desired
        # time.
        Gnm_e0 = G[kg,idx+2]

        if date < dat_year
            Gnm_e1 = G[kg,idx+3]
            ΔG     = (Gnm_e1-Gnm_e0)/5
        else
            ΔG     = G[kg,gend]
        end

        Gnm  = Gnm_e0 + ΔG*Δt
        kg  += 1

        aux_dVr += -(n+1)/r*Gnm*P[n+1,1]
        aux_dVθ += Gnm*dP[n+1,1]

        # Sine and cosine with m = 1
        # ==========================
        #
        # This values will be used to update recursively `sin(m*ϕ)` and
        # `cos(m*ϕ)`, reducing the computational burden.
        sin_mϕ   = +sin_ϕ    # sin( 1*λ_gc)
        sin_m_1ϕ = 0.0       # sin( 0*λ_gc)
        sin_m_2ϕ = -sin_ϕ    # sin(-1*λ_gc)
        cos_mϕ   = +cos_ϕ    # cos( 1*λ_gc)
        cos_m_1ϕ = 1.0       # cos( 0*λ_gc)
        cos_m_2ϕ = +cos_ϕ    # cos(-2*λ_gc)

        # Other auxiliary variables that depend only on `n`
        # =================================================

        fact_dVr = (n+1)/r

        # Compute the contributions when `m ∈ [1,n]`
        # ==========================================

        for m = 1:n
            # Compute recursively `sin(m*ϕ)` and `cos(m*ϕ)`.
            sin_mϕ = 2cos_ϕ*sin_m_1ϕ-sin_m_2ϕ
            cos_mϕ = 2cos_ϕ*cos_m_1ϕ-cos_m_2ϕ

            # Compute the coefficients `G_nm` and `H_nm`
            # ==========================================

            # Get the coefficients in the epoch and interpolate to the
            # desired time.
            Gnm_e0 = G[kg,idx+2]
            Hnm_e0 = H[kh,idx+2]

            if date < dat_year
                Gnm_e1 = G[kg,idx+3]
                Hnm_e1 = H[kh,idx+3]
                ΔG     = (Gnm_e1-Gnm_e0)/5
                ΔH     = (Hnm_e1-Hnm_e0)/5
            else
                ΔG     = G[kg,gend]
                ΔH     = H[kh,hend]
            end

            Gnm    = Gnm_e0 + ΔG*Δt
            Hnm    = Hnm_e0 + ΔH*Δt
            kg    += 1
            kh    += 1

            GcHs_nm = Gnm*cos_mϕ + Hnm*sin_mϕ
            GsHc_nm = Gnm*sin_mϕ - Hnm*cos_mϕ

            # Compute the contributions for `m`
            # =================================

            aux_dVr += -fact_dVr*GcHs_nm*P[n+1,m+1]
            aux_dVθ += GcHs_nm*dP[n+1,m+1]
            aux_dVϕ += (θ == 0) ? -m*GsHc_nm*dP[n+1,m+1] : -m*GsHc_nm*P[n+1,m+1]

            # Update the values for the next step
            # ===================================

            sin_m_2ϕ = sin_m_1ϕ
            sin_m_1ϕ = sin_mϕ
            cos_m_2ϕ = cos_m_1ϕ
            cos_m_1ϕ = cos_mϕ
        end

        # Perform final computations related to the summation in `n`
        # ==========================================================

        # fact = (a/r)^(n+1)
        fact    *= ratio

        # aux_<> *= (a/r)^(n+1)
        aux_dVr *= fact
        aux_dVϕ *= fact
        aux_dVθ *= fact

        dVr += aux_dVr
        dVϕ += aux_dVϕ
        dVθ += aux_dVθ
    end

    dVr *= a
    dVϕ *= a
    dVθ *= a

    # Compute the Geomagnetic field vector in the geocentric reference frame
    # ======================================================================

    x = +1/r*dVθ
    y = (θ == 0) ? -1/r*dVϕ : -1/(r*sin(θ))*dVϕ
    z = dVr

    B_gc = SVector{3,Float64}(x,y,z)

    return B_gc
end

function igrf(date::Number, h::Number, λ::Number, Ω::Number, ::Val{:geodetic},
              P::AbstractMatrix{T}, dP::AbstractMatrix{T};
              show_warns = true) where T<:Real

    # TODO: This method has a small error (≈ 0.01 nT) compared with the
    # `igrf12syn`. However, the result is exactly the same as the MATLAB
    # function in [3]. Hence, this does not seem to be an error in the
    # conversion from geodetic to geocentric coordinates. This is probably
    # caused by a numerical error. Further verification is necessary.

    # Convert the geodetic coordinates to geocentric coordinates.
    (λ_gc, r) = geodetic_to_geocentric(λ, h)

    # Compute the geomagnetic field in geocentric coordinates.
    B_gc = igrf(date, r, λ_gc, Ω, Val(:geocentric), P, dP; show_warns = show_warns)

    # Convert to geodetic coordinates.
    D_gd_gc = create_rotation_matrix(λ_gc - λ,:Y)
    B_gd    = D_gd_gc*B_gc

    return B_gd
end