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

"""
    Scenario(d::Dict)

Constructor for Scenario struct, where `d` has upper-case keys:
- Site (required)
- ElectricTariff (required)
- ElectricLoad (required)
- PV (optional, can be Array)
- Wind (optional)
- Storage (optional)
- ElectricUtility (optional)
- Financial (optional)
- Generator (optional)
- DomesticHotWaterLoad (optional)
- SpaceHeatingLoad (optional)
- ExistingBoiler (optional)

All values of `d` are expected to be `Dicts` except for `PV`, which can be either a `Dict` or `Dict[]`.
```
struct Scenario
    settings::Settings
    site::Site
    pvs::Array{PV, 1}
    wind::Wind
    storage::Storage
    electric_tariff::ElectricTariff
    electric_load::ElectricLoad
    electric_utility::ElectricUtility
    financial::Financial
    generator::Generator
    dhw_load::DomesticHotWaterLoad
    space_heating_load::SpaceHeatingLoad
    existing_boiler::ExistingBoiler
end
```
"""
function Scenario(d::Dict)
    if haskey(d, "Settings")
        settings = Settings(;dictkeys_tosymbols(d["Settings"])...)
    else
        settings = Settings()
    end

    site = Site(;dictkeys_tosymbols(d["Site"])...)

    pvs = PV[]
    if haskey(d, "PV")
        if typeof(d["PV"]) <: AbstractArray
            for (i, pv) in enumerate(d["PV"])
                check_pv_tilt!(pv, site)
                if !(haskey(pv, "name"))
                    pv["name"] = string("PV", i)
                end
                push!(pvs, PV(;dictkeys_tosymbols(pv)...))
            end
        elseif typeof(d["PV"]) <: AbstractDict
            check_pv_tilt!(d["PV"], site)
            push!(pvs, PV(;dictkeys_tosymbols(d["PV"])...))
        else
            error("PV input must be Dict or Dict[].")
        end
    end

    if haskey(d, "Financial")
        financial = Financial(; dictkeys_tosymbols(d["Financial"])...)
    else
        financial = Financial()
    end

    if haskey(d, "ElectricUtility")
        electric_utility = ElectricUtility(; dictkeys_tosymbols(d["ElectricUtility"])...)
    else
        electric_utility = ElectricUtility()
    end

    if haskey(d, "Storage")
        # only modeling electrochemical storage so far
        storage_dict = Dict(:elec => dictkeys_tosymbols(d["Storage"]))
        storage = Storage(storage_dict, financial)
    else
        storage_dict = Dict(:elec => Dict(:max_kw => 0))
        storage = Storage(storage_dict, financial)
    end

    electric_load = ElectricLoad(; dictkeys_tosymbols(d["ElectricLoad"])...,
                                   latitude=site.latitude, longitude=site.longitude,
                                   time_steps_per_hour=settings.time_steps_per_hour
                                )

    electric_tariff = ElectricTariff(; dictkeys_tosymbols(d["ElectricTariff"])...,
                                       year=electric_load.year,
                                       NEM=electric_utility.net_metering_limit_kw > 0,
                                       time_steps_per_hour=settings.time_steps_per_hour
                                    )

    if haskey(d, "Wind")
        wind = Wind(; dictkeys_tosymbols(d["Wind"])...,
                    average_elec_load=sum(electric_load.loads_kw) / length(electric_load.loads_kw))
    else
        wind = Wind(; max_kw=0)
    end

    if haskey(d, "Generator")
        generator = Generator(; dictkeys_tosymbols(d["Generator"])...)
    else
        generator = Generator(; max_kw=0)
    end

    max_heat_demand_kw = 0.0
    if haskey(d, "DomesticHotWaterLoad")
        add_doe_reference_names_from_elec_to_thermal_loads(d["ElectricLoad"], d["DomesticHotWaterLoad"])
        dhw_load = DomesticHotWaterLoad(; dictkeys_tosymbols(d["DomesticHotWaterLoad"])...,
                                          latitude=site.latitude, longitude=site.longitude,
                                          time_steps_per_hour=settings.time_steps_per_hour
                                        )
        max_heat_demand_kw = maximum(dhw_load.loads_kw)
    else
        dhw_load = DomesticHotWaterLoad(; fuel_loads_mmbtu_per_hour=repeat([0.0], 8760))
    end

    if haskey(d, "SpaceHeatingLoad")
        add_doe_reference_names_from_elec_to_thermal_loads(d["ElectricLoad"], d["SpaceHeatingLoad"])
        space_heating_load = SpaceHeatingLoad(; dictkeys_tosymbols(d["SpaceHeatingLoad"])...,
                                                latitude=site.latitude, longitude=site.longitude,
                                                time_steps_per_hour=settings.time_steps_per_hour
                                              )

        max_heat_demand_kw = maximum(space_heating_load.loads_kw .+ max_heat_demand_kw)
    else
        space_heating_load = SpaceHeatingLoad(; fuel_loads_mmbtu_per_hour=repeat([0.0], 8760))
    end

    if max_heat_demand_kw > 0
        vals = Dict{Symbol, Any}()
        vals[:max_heat_demand_kw] = max_heat_demand_kw
        vals[:time_steps_per_hour] = settings.time_steps_per_hour
        if haskey(d, "ExistingBoiler")
            vals = merge(vals, dictkeys_tosymbols(d["ExistingBoiler"]))
        end
        existing_boiler = ExistingBoiler(; vals...)
    else
        existing_boiler = ExistingBoiler(0.0, 0.0, Real[])
    end

    return Scenario(
        settings,
        site,
        pvs,
        wind,
        storage,
        electric_tariff,
        electric_load,
        electric_utility,
        financial,
        generator,
        dhw_load,
        space_heating_load,
        existing_boiler
    )
end

"""
    Scenario(fp::String)

Consruct Scenario from filepath `fp` to JSON with keys aligned with the `Scenario(d::Dict)` method.
"""
function Scenario(fp::String)
    Scenario(JSON.parsefile(fp))
end

function check_pv_tilt!(pv::Dict, site::Site)
    if !(haskey(pv, "tilt"))
        pv["tilt"] = site.latitude
    end
end

function add_doe_reference_names_from_elec_to_thermal_loads(elec::Dict, thermal::Dict)
    string_keys = [
        "doe_reference_name",
        "blended_doe_reference_names",
        "blended_doe_reference_percents",
    ]
    for k in string_keys
        if !(k in keys(thermal)) && k in keys(elec)
            thermal[k] = elec[k]
        end
    end
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
using libm
srand(1)

X = exp(randn(10_000_000))

macro testsum(name,fn,X)
    quote
        function test(X)
            s = 0.0
            for x in X
                s += $fn(x)
            end
            s
        end

        test($X)
        test($X)

        println($name)
        @time test($X)
        println()
    end
end

syslog(x::Float64) = Core.Intrinsics.nan_dom_err(ccall((:log,:libm),Float64,(Float64,),x), x)
syslog(x::Float32) = Core.Intrinsics.nan_dom_err(ccall((:logf,:libm),Float32,(Float32,),x), x)

# println("Float64")
# println("-------------------")
# @testsum "Openlibm:" log X
# @testsum "System:" syslog X

# @testsum "Julia fdlibm:" libm.log X
# @testsum "Julia tang:" libm.log_tang X
# # @testsum "muladd:" libm.log_muladd X

# println()
# println()

X = float32(X)

println("Float32")
println("-------------------")
@testsum "Openlibm:" log X
@testsum "System:" syslog X

#@testsum "Julia fdlibm:" libm.log X
@testsum "Julia tang:" libm.log_tang X

ext: jl


text:
"""
    rasterize(
        res::Tuple{Int64, Int64},
        points::Matrix{Float64},
        colors::Matrix{Float64};
        xlim = (minimum(points[:,1]), maximum(points[:,1])),
        ylim = (minimum(points[:,2]), maximum(points[:,2]))
        )::Array{Float64, 3}

Rasterize the `points` represented as coordinates in 2-row matrix colored by
RGBA float colors in the corresponding 4-row matrix into a raster of size
`res`. `xlim` and `ylim` are the point-coordinates of raster edges.

The raster is returned as a 3D RGBA array with coordinates as in
`[colorChannel, x, y]`.
"""
function rasterize(
    res::Tuple{Int64,Int64},
    points::Matrix{Float64},
    colors::Matrix{Float64};
    xlim = (minimum(points[1, :]), maximum(points[1, :])),
    ylim = (minimum(points[2, :]), maximum(points[2, :])),
)::Raster

    # Color blending is computed on integers because floats are just slower.
    # `precision` is where the value of float 1.0 is mapped.
    # Note that because there's multiplication in the process, the value of
    # precision should not exceed the square root of maximum integer precision,
    # in this case around 2^31. 100k as used here is MORE THAN ENOUGH for any
    # thinkable use in this universe.
    precision = Int64(100000)
    # Most notably, having precision > 256 dodges the ugly low-alpha color
    # artifacts that are pretty common with many graphics libraries.

    # integer raster
    ra = zeros(Int64, 4, res[1], res[2])

    # minima
    mins = [xlim[1], ylim[1]]
    # actual raster point sizes
    iszs = [res[1] / (xlim[2] - xlim[1]), res[2] / (ylim[2] - ylim[1])]

    if size(points,1) != 2 || size(colors) != (4, size(points,2))
        throw(ArgumentError("wrong input matrix sizes"))
    end

    @inbounds for i = 1:size(points, 2)
        # convert the input coordinates to raster indexes
        posx =  trunc(Int64, (points[1, i] - mins[1]) * iszs[1]) + 1
        posy =  trunc(Int64, (points[2, i] - mins[2]) * iszs[2]) + 1

        # skip if it's off limits
        if posx < 1 || posx > res[1] || posy < 1 || posy > res[2]
            continue
        end

        # get the src color
        srcr = trunc(Int64, precision * colors[1, i])
        srcg = trunc(Int64, precision * colors[2, i])
        srcb = trunc(Int64, precision * colors[3, i])
        srca = trunc(Int64, precision * colors[4, i])

        # premultiply alpha
        srcr = (srcr * srca) ÷ precision
        srcg = (srcg * srca) ÷ precision
        srcb = (srcb * srca) ÷ precision

        # get the dst color
        dstr = ra[1, posx, posy]
        dstg = ra[2, posx, posy]
        dstb = ra[3, posx, posy]
        dsta = ra[4, posx, posy]

        # aaaand it blends!
        ra[1, posx, posy] = srcr + dstr - ((srca * dstr) ÷ precision)
        ra[2, posx, posy] = srcg + dstg - ((srca * dstg) ÷ precision)
        ra[3, posx, posy] = srcb + dstb - ((srca * dstb) ÷ precision)
        ra[4, posx, posy] = srca + dsta - ((srca * dsta) ÷ precision)
    end

    #unmultiply alpha and convert back to floats in range 0..1
    raf = zeros(Float64, 4, res[1], res[2])

    for i = 1:res[1]
        for j = 1:res[2]
            # get the alpha element
            raa = ra[4, i, j]
            if raa == 0 #transparent black
                continue
            end
            # remove alpha (which is still in precision-scale, so this
            # effectively removes precision multiplication from r,g,b channels
            raf[1, i, j] = ra[1, i, j] / raa
            raf[2, i, j] = ra[2, i, j] / raa
            raf[3, i, j] = ra[3, i, j] / raa
            # and scale alpha back down to 0..1
            raf[4, i, j] = ra[4, i, j] / precision
        end
    end

    raf
end


text:
import Bridge: ProcessOrCoefficients

"""
    invSolve!(::EulerMaruyama, Y, W::SamplePath, P::ProcessOrCoefficients)

Compute Wiener process `W` that would have been needed to be simulated in order
to obtain path `Y` under law `P` when using Euler-Maruyama numerical scheme
"""
function inv_solve!(::EulerMaruyama, Y, W::SamplePath, P::ProcessOrCoefficients
                  ) where {T}
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = W.tt

    yy = Y.yy
    #tt[:] = Y.tt
    ww[.., N] = zero(ww[.., N])

    for i in N-1:-1:1
        yᵢ₊₁ = yy[.., i+1]
        yᵢ = yy[.., i]

        ww[.., i] = ww[.., i+1] - σ(tt[i], yᵢ, P)\(yᵢ₊₁ - yᵢ - _b((i,tt[i]), yᵢ, P)*(tt[i+1]-tt[i]))
    end
    for i in N:-1:1
        ww[.., i] = ww[.., i] - ww[.., 1]
    end
    W
end


text:
# Script performing value iteration with tau process
#
# This is a fixed / extended version of the VI module with correct(ed) boundary
# conditions

module VI2

import Roots: fzero, ConvergenceFailed
import Distributions: Gamma, pdf, quantile
import Base: size
using NLopt

using CIRProc

export FDTask, RTTask, getrr,  vi, gettaus, getgs, pdf,
    max_rr_particle,
    XTauProcSpace, GTauProcSpace

# discretised parameter space for VI
abstract AbstractProcSpace
# PDE solver to compute the expectation
abstract AbstractPDE2DSolver

# ----------------------------------------------------------------------------
# Task definition
# ----------------------------------------------------------------------------

abstract AbstractTask

# fixed duration task (infinite duration)
immutable FDTask <: AbstractTask
    p::CIRProcess
    c::Float64
    rcorr::Float64
    rincorr::Float64

    function FDTask(p::CIRProcess, c::Real, rcorr::Real, rincorr::Real)
        new(p, float(c), float(rcorr), float(rincorr))
    end
    FDTask(p, c) = FDTask(p, c, 1, 0)
    FDTask(p) = FDTask(p, 0.1)
end
getrr(t::FDTask, pc::Real, dt::Real) =
    float(pc) * t.rcorr + (1.0 - float(pc)) * t.rincorr - float(dt) * t.c

# reaction time task
immutable RTTask <: AbstractTask
    p::CIRProcess
    c::Float64
    rcorr::Float64
    rincorr::Float64
    ti::Float64
    tp::Float64

    function RTTask(p::CIRProcess,
                    c::Real, rcorr::Real, rincorr::Real, ti::Real, tp::Real)
        ti, tp = float(ti), float(tp)
        ti >= 0 || error("ti needs to be non-negative")
        tp >= 0 || error("tp needs to be non-negative")
        new(p, float(c), float(rcorr), float(rincorr), ti, tp)
    end
    RTTask(p, c, ti) = RTTask(p, c, 1, 0, ti, 0)
    RTTask(p, c) = RTTask(p, c, 1)
    RTTask(p) = RTTask(p)
end
getrr(t::RTTask, pc::Real, dt::Real) =
    (float(pc) * t.rcorr + (1.0 - float(pc)) * t.rincorr - float(dt) * t.c) /
    (dt + t.ti + (1.0 - float(pc)) * t.tp)

# ----------------------------------------------------------------------------
# Heuristic bound optimization
# ----------------------------------------------------------------------------

include("rr_optim.jl")

# ----------------------------------------------------------------------------
# Dynamic programming
# ----------------------------------------------------------------------------

# Performs value iteration to find the optimal bound
#
# t is the task description
# r is the nx x ny discretized process parameter space
# dt is the time discretisation
#
# The function returns (bound, v, er), where bound is an nx-element
# array, v is the nx x ny value function, and er is the expected reward
function vi(t::FDTask, r::AbstractProcSpace, dt::Float64)
    rdiff = t.rcorr - t.rincorr
    tnorm = FDTask(t.p, t.c / rdiff)
    bound, v = vinorm(t, r, dt)
    v = rdiff * v .+ t.rincorr
    er = avgreward(v, gettaus(r), t.p)
    return bound, v, er
end

# Performs value iteration to find the optimal bound
#
# The function assumes
#
# t is the task description
# r is the nx x ny discretized process parameter space
# dt is the time discretisation
#
# The function returns (bound, v, rho), where bound is an ny-element
# array, v is the nx x ny value function, and rho is the reward rate
function vi(t::RTTask, r::AbstractProcSpace, dt::Float64; verbose::Bool=false)
    const maxiter, roottol = 20, 1e-12
    rdiff = t.rcorr - t.rincorr
    bestrho, bestv0 = 0.0, Inf

    function rr_initial_value(rho::Real, rc::AbstractProcSpace)
        k = rdiff + rho * t.tp
        tnorm = FDTask(t.p, (t.c + rho) / k)
        bound, v = vinorm(tnorm, rc, dt)
        v0 = avgreward(v * k .+ t.rincorr .- rho * (t.ti + t.tp), gettaus(rc), t.p)
        !verbose || println("rho = $rho, <v(1/2,tau)> = $v0")
        if abs(v0) < bestv0; bestrho, bestv0 = rho, abs(v0); end
        v0
    end

    # reward rate bounds, for DT from 0 to Inf, R from rincorr to rcorr
    rrs = [t.rcorr / t.ti, (t.rcorr + t.rincorr) / (2 * t.ti + t.tp), -t.c]
    rrmax, rrmin = maximum(rrs), minimum(rrs)
    rrini = rrmin + 0.8 * (rrmax - rrmin)

    # find reward rate that makes v(1/2) = 0
    const nx, ny = size(r)
    const xdiv, ydiv = 5, 5
    rcoarse = coarser(r, xdiv, ydiv)
    const nxcoarse, nycoarse = size(rcoarse)
    !verbose || println("-- Root finding to find reward rate")
    !verbose || println("-- Coarse grid, nx = $nxcoarse, ny = $nycoarse")
    # make sure to return last rho, even if convergence failed
    rho = try
        fzero(x -> rr_initial_value(x, rcoarse),
              rrini, [rrmin, rrmax], max_iter=maxiter, tol=roottol)
    catch y
        if isa(y, ConvergenceFailed)
            bestrho
        else
            rethrow(y)
        end
    end
    !verbose || println("-- Fine grid, nx = ", nx, ", ny = ", ny)
    bestv0 = Inf
    rho = try
        fzero(x -> rr_initial_value(x, r),
              rho, rho .+ 0.1 * (rrmax - rrmin) * [-1, 1],
              max_iter=maxiter, tol=roottol)
    catch y
        if isa(y, ConvergenceFailed)
            bestrho
        else
            rethrow(y)
        end
    end
    !verbose || println("-- Reward rate rho = $rho")

    k = rdiff + rho * t.tp;
    tnorm = FDTask(t.p, (t.c + rho) / k)
    bound, v = vinorm(tnorm, r, dt)

    return bound, v * k .+ t.rincorr .- rho * (t.ti + t.tp), rho
end

# returns the tau discretisation
#
# ntau is the number of tau steps
# maxtau is the largest tau to represent
# t is the task, in which case tau is discretised up to 2 * 99th prctile
#
# The function returns a vector of size ntau
gettaus(ntau::Int, maxtau::Float64) = linspace(0.0, maxtau, ntau)
gettaus(ntau::Int, p::CIRProcess) = gettaus(ntau, 2*quantile(p, 0.99))

# returns the g discretisation
#
# ng is the number of g steps
#
# The function returns a vector of size ng
getgs(ng::Int) = linspace(0.0, 1.0, ng)

# performs value iteration for fixed duration tasks and normalised costs
#
# t is the task description. Only the process and cost are considered
# r is the nx x ny discretized process parameter space
# dt is the time discretisation
#
# The function returns (v, bound), where v is the nx x ny final value function,
# and bound is the bound found over ny-space.
function vinorm(t::FDTask, r::AbstractProcSpace, dt::Float64)
    p = PDE2DADISolver(r, dt)
    vini, vd, gs = getvini(r), getvd(r), getgs(r)
    vinormiter(p, vini, vd, gs, t.c, dt)
end

# returns value function and policy resulting from value iteration
#
# p is the solver to find the expected value function
# vini is the K x J matrix specifying the initial value function
# vd is the K x J matrix speficying the value function for immediate deicsions
# gs is the K-vector of beliefs
# c is the cost per unit time
# dt is the time discretisation (needs to match p)
# maxiter specifies the maximum number of VI steps
# maxvdiff is the largest acceptable abs(v - vprev).dt to consider convergence
#
# The function returns (v, bound), where v is the K x J final value function,
# and bound is the found bound over J-space. This bound is computed by
# interpolating over values in gs. It is, in fact, the only point where gs is
# being used.
#
# vini and vd need to correspond to each other at the boundaries in k, but
# not in j.
function vinormiter{T <: Real}(p::AbstractPDE2DSolver,
    vini::Matrix{T}, vd::Matrix{T}, gs::Vector{T}, c::Real, dt::Real,
    maxiter::Int=0, maxvdiff::Real=0.001)

    const K, J = size(vini, 1), size(vini, 2)
    @assert K == size(vd, 1) && J == size(vd, 2)
    @assert K == length(gs)
    @assert dt > zero(dt)
    if maxiter == 0
        maxiter = iceil(50.0/dt)
    end

    const dc = c * dt
    converged = false
    v = vini
    vnext = Array(Float64, K, J)
    for n = 2:maxiter
        vprev = copy(v)
        vnext = expectedu(p, v)

        # bounds except for j = 1 need to remain unchanged
        @assert all([isapprox(vnext[k,end], vprev[k,end]) for k = 1:K])
        @assert all([isapprox(vnext[1,j], vprev[1,j]) for j = 1:J])
        @assert all([isapprox(vnext[end,j], vprev[end,j]) for j = 1:J])
        # boundaries for k = 1 and k = K need to be decision boundaries
        @assert all([isapprox(vnext[1,j], vd[1,j]) for j = 1:(J-1)])
        @assert all([isapprox(vnext[end,j], vd[end,j]) for j=  1:(J-1)])

        # based on vnext, update value function, using unrolled
        # v = max(vnext .- dc, vd). v[:,end] is unchanged to preserve
        # boundary condition at j=J
        for j = 1:(J-1), k = 1:K
            v[k,j] = max(vnext[k,j] - dc, vd[k,j])
        end

        # bounds except for j = 1 need to remain unchanged
        @assert all([isapprox(v[k,end], vprev[k,end]) for k = 1:K])
        @assert all([isapprox(v[1,j], vprev[1,j]) for j = 1:J])
        @assert all([isapprox(v[end,j], vprev[end,j]) for j = 1:J])

        # compute value function difference
        if maxvdiff > zero(maxvdiff)
            vdiff = -Inf
            for j=1:J, k=1:K
                const vdiffjk = abs(v[k,j] - vprev[k,j])
                if vdiff < vdiffjk; vdiff = vdiffjk; end
            end
            if vdiff/dt <= maxvdiff; converged = true;  break;  end
        end
    end
    if maxvdiff > zero(maxvdiff) && !converged
        warn("value iteration did not converge within $maxiter iterations")
    end

    interpolatebound(vnext .- dc, vd, gs), v
end

# Returns the decision value function for (1,0) reward
#
# gs is ng-element array of reals
#
# function returns ng x n array of the decision value function
function getvd(gs::Vector, n::Int)
    repmat(max(gs, 1.-gs), 1, n)
end

# Returns initial value function proposal over g's for given bound
#
# bound defines the bound location
# gs is ng-element array of reals
#
# function returns initial value function over g's that corresponds to bound
# at given location.
function getvini{T <: Real}(bound::T, gs::Vector{T})
    vini = 1.0/(2*bound-1)^3 * (bound^3*(3*bound-2) .-
        (1+6*bound*(bound-1))*gs .+ 6*bound*(bound-1)*gs.^2 .+ 2*gs.^3 .- gs.^4)
    max(vini, max(gs, 1.-gs))
end

# vnext is ng x n_nu array
# vd is ng x n_nu array
# gs is ng-element array of reals
#
# function returns bound on intersecting vnext and vd, using linear interp.
#
# the returned bound is a vector is size n_nu
function interpolatebound{T <: Real}(vnext::Matrix{T}, vd::Matrix{T}, gs::Vector{T})
    # argument properties
    ng, n_nu = size(vnext)
    dg = gs[2] - gs[1]

    # compute bound for each nu separately
    bound = fill(NaN, n_nu)
    g_lo = div(ng,2)+1
    for n = 1:n_nu
        g_idx = find(vd[g_lo:end,n] .>= vnext[g_lo:end,n])
        if isempty(g_idx)
            # vd >= vnext nowhere -> bound is at g = 1
            bound[n] = 1.0
        elseif g_idx[1] == 1
            # vd >= vnext everywhere -> bound is at g = 0.5
            bound[n] = 0.5
        else
            # adjust g_idx to be gs index before/below boundary
            g_idx = g_idx[1]+g_lo-2
            # find boundary by linear interpolation between g_idx and g_idx+1
            bound[n] = gs[g_idx] + dg/(
                1+(vd[g_idx+1,n]-vnext[g_idx+1,n])/(vnext[g_idx,n]-vd[g_idx,n]))
        end
    end

    bound
end

# Returns the expected reward at trial start
#
# The expected return is V(g,tau) at g=1/2 with the expectation over tau,
# assuming that tau follows the prior p(tau).
#
# v is the ng x ntau value function
# taus is the tau discretisation
# p are the CIR process parameters
#
# The function returns the scalar expected reward.
function avgreward{T <: Real}(v::Matrix{T}, taus::Vector{T}, p::CIRProcess)
    # argument properties
    ng, ntau = size(v)
    @assert length(taus) == ntau

    # find values at belief g=0.5
    vg2 = fill(NaN, ntau)
    if mod(ng,2) == 0
        g2_idx = div(ng,2)
        vg2 = squeeze(mean(v[g2_idx:(g2_idx+1),:],1),1)
    else
        vg2 = v[div(ng,2)+1,:]
    end

    empmean(p, taus, vg2)
end

# ----------------------------------------------------------------------------
# Process spaces
# ----------------------------------------------------------------------------

# Process space over (X, tau)
immutable XTauProcSpace <: AbstractProcSpace
    xs::Vector{Float64}
    gs::Vector{Float64}
    taus::Vector{Float64}
    p::CIRProcess

    function XTauProcSpace(p::CIRProcess, nx::Int, ntau::Int, maxtau::Real)
        maxtau > zero(maxtau) || error("maxtau needs to be positive")
        const xrange = [-1.0, 1] * 0.5log(2.0nx - 1.0) # 1/2nx distance from [0,1]
        const xs = linspace(xrange[1], xrange[2], nx)
        new(xs, 1.0 ./ (1.0 .+ exp(-2.0xs)), gettaus(ntau, maxtau), p)
    end
end
XTauProcSpace(p::CIRProcess, nx::Int, ntau::Int) =
    XTauProcSpace(p, nx, ntau, 2.0quantile(p, 0.99))
size(r::XTauProcSpace) = length(r.xs), length(r.taus)
getgs(r::XTauProcSpace) = r.gs
gettaus(r::XTauProcSpace) = r.taus
getvd(r::XTauProcSpace) = getvd(r.gs, length(r.taus))
function getvini(r::XTauProcSpace)
    const nx, ntau = size(r)
    const gmax = max(r.gs[1], r.gs[end])
    vini = Array(Float64, nx, ntau)
    for j = 1:ntau
        # upper rim of v for large |x| and tau is at gmax
        vbase = getvini(0.5+(gmax-0.5)*(j-1)/(ntau-1), r.gs)
        vini[:,j] = (vbase .- gmax)*(ntau-j)/(ntau-1) .+ gmax;
    end
    vini
end
function moments(r::XTauProcSpace, dt::Real)
    const nx, ntau = size(r)
    const dx = (2.0r.gs .- 1.0) .* r.taus'             # (2g(X)-1) tau
    const dx2 = repmat(r.taus', nx, 1)                 # tau
    const dtau = repmat(r.p.theta * (r.p.mu .- r.taus'), nx, 1) # theta (mu-tau)
    const dtau2 = repmat(2.0r.p.theta * r.p.sig^2 / r.p.mu * r.taus',
                         nx, 1)                        # 2 theta sig^2 / mu tau
    dx, dx2, dtau, dtau2, zeros(nx, ntau), r.xs[2]-r.xs[1], r.taus[2]-r.taus[1]
end
coarser(r::XTauProcSpace, xdiv::Real, ydiv::Real) =
    XTauProcSpace(r.p, int(div(length(r.xs),xdiv)), int(div(length(r.taus),ydiv)),
                  maximum(r.taus))

# Process space over (g, tau)
immutable GTauProcSpace <: AbstractProcSpace
    gs::Vector{Float64}
    taus::Vector{Float64}
    p::CIRProcess

    function GTauProcSpace(p::CIRProcess, ng::Int, ntau::Int, maxtau::Real)
        maxtau > zero(maxtau) || error("maxtau needs to be positive")
        new(linspace(0.0, 1.0, ng), gettaus(ntau, maxtau), p)
    end
end
GTauProcSpace(p::CIRProcess, ng::Int, ntau::Int) =
    GTauProcSpace(p, ng, ntau, 2.0quantile(p, 0.99))
size(r::GTauProcSpace) = length(r.gs), length(r.taus)
getgs(r::GTauProcSpace) = r.gs
gettaus(r::GTauProcSpace) = r.taus
getvd(r::GTauProcSpace) = getvd(r.gs, length(r.taus))
function getvini(r::GTauProcSpace)
    const ng, ntau = size(r)
    vini = Array(Float64, ng, ntau)
    for j = 1:ntau
        vbase = getvini(0.5+0.5*(j-1)/(ntau-1), r.gs)
        vini[:,j] = (vbase .- 1.0)*(ntau-j)/(ntau-1) .+ 1.0;
    end
    vini
end
function moments(r::GTauProcSpace, dt::Real)
    const ng, ntau = size(r)
    const dg2 = (4.0(1.0 .- r.gs).^2 .* r.gs.^2) .* r.taus' # 4(1-g)^2 g^2 tau
    const dtau = repmat(r.p.theta * (r.p.mu .- r.taus'), ng, 1) # theta (mu-tau)
    const dtau2 = repmat(2.0r.p.theta * r.p.sig^2 / r.p.mu * r.taus',
                         ng, 1)                            # theta sig^2 / mu tau
    zeros(ng, ntau), dg2, dtau, dtau2, zeros(ng, ntau),
    r.gs[2] - r.gs[1], r.taus[2] - r.taus[1]
end
coarser(r::GTauProcSpace, xdiv::Real, ydiv::Real) =
    GTauProcSpace(r.p, int(div(length(r.gs),xdiv)), int(div(length(r.taus),ydiv)),
                  maximum(r.taus))

# ----------------------------------------------------------------------------
# PDE solvers to compute the expected value
# ----------------------------------------------------------------------------

# all solvers operate over a 2D space, equally spaced in each of the two
# dimensions (but spacing does not need to be equal across dimensions).
# It assumes the boundary conditions
# - constant boundaries in x
# - constant upper boundary in y
# - lose lower boundary in y, assuming <dx> = <dx2> = <dy2> = 0 and <dy> > 0

# PDE solver using single matrix inversion over collapsed 2D space
immutable PDE2DSolver <: AbstractPDE2DSolver
    ukj::Matrix{Float64}
    ukpj::Matrix{Float64}
    uknj::Matrix{Float64}
    ukjp::Matrix{Float64}
    ukjn::Matrix{Float64}
    ukdjd::Matrix{Float64}
    L::SparseMatrixCSC{Float64,Int}

    function PDE2DSolver{T <: Real}(dx::Matrix{T}, dx2::Matrix{T},
        dy::Matrix{T}, dy2::Matrix{T}, dxy::Matrix{T},
        Dx::Real, Dy::Real, dt::Real)
        Dx, Dy, dt = float(Dx), float(Dy), float(dt)
        Dx > 0.0 || error("Dx needs to be positive")
        Dy > 0.0 || error("Dy needs to be positive")
        dt > 0.0 || error("dt needs to be positive")
        const K, J = size(dx, 1), size(dx, 2)
        K == size(dy, 1) && J == size(dy, 2) || error("dx and dy need to be of same size")
        K == size(dx2, 1) && J == size(dx2, 2) || error("dx and dx2 need to be of same size")
        K == size(dy2, 1) && J == size(dy2, 2) || error("dx and dy2 need to be of same size")
        K == size(dxy, 1) && J == size(dxy, 2) || error("dx and dxy need to be of same size")
        const dtDx, dtDy, dtDx2, dtDy2 = dt/Dx, dt/Dy, dt/(Dx*Dx), dt/(Dy*Dy)
        const ukj = 0.5(dtDx2*float(dx2) .+ dtDy2*float(dy2))
        const ukpj = 0.25(dtDx*float(dx) .+ dtDx2*float(dx2))
        const uknj = 0.25(dtDx2*float(dx2) .- dtDx*float(dx))
        const ukjp = 0.25(dtDy*float(dy) .+ dtDy2*float(dy2))
        const ukjn = 0.25(dtDy2*float(dy2) .- dtDy*float(dy))
        const ukdjd = dt/(8.0Dx*Dy)*float(dxy)
        new(ukj, ukpj, uknj, ukjp, ukjn, ukdjd,
            _get2DSolverL(ukj, ukpj, uknj, ukjp, ukjn, ukdjd))
    end
end
function PDE2DSolver(r::AbstractProcSpace, dt::Real)
    dt > zero(dt) || error("dt needs to be positive")
    const dx, dx2, dy, dy2, dxy, Dx, Dy = moments(r, dt)
    PDE2DSolver(dx, dx2, dy, dy2, dxy, Dx, Dy, dt)
end
size(p::PDE2DSolver) = size(p.ukj)
expectedu{T <: Real}(p::PDE2DSolver, u::Matrix{T}) =
    reshape(p.L \ _getb(p, float(u)), size(u, 1), size(u, 2))

# returns L matrix for PDE2DSolver
#
# The given matrices are of size K x J, containing various pre-computed
# constants.
#
# The returned sparse L is of size KJ x KJ
function _get2DSolverL(ukj, ukpj, uknj, ukjp, ukjn, ukdjd)
    const K, J = size(ukj)
    const maxels = 3K - 2 + 9(J-2)*(K-2) + 2(J-2)
    row, col, v = Array(Int, maxels), Array(Int, maxels), Array(Float64, maxels)
    # fill in-between block elements
    addnonzero(ri, ci, vi, i) = v == 0.0 ? i : begin   # only add non-zero elements
                                                   row[i], col[i], v[i] = ri, ci, vi
                                                   i + 1
                                               end
    # diagonals of first two block in top block row determine j=1 boundary
    row[1], col[1], v[1] = 1, 1, 1.0
    for k = 2:(K-1)
        const ki = 2k - 2
        row[ki], col[ki], v[ki] = k, k, 1.0 + 2.0ukjp[k,1]
        row[ki+1], col[ki+1], v[ki+1] = k, k+K, -2.0ukjp[k,1]
    end
    i = 2(K-1) + 1
    row[i-1], col[i-1], v[i-1] = K, K, 1.0
    # all j > 1 blocks
    for j = 2:(J-1)
        # first row has 1 along diagonal
        i = addnonzero((j-1)K + 1, (j-1)K + 1, 1.0, i)
        # other rows
        for k = 2:(K-1)
            const ri = (j-1)K + k  # row of L
            const ukdjdkj = ukdjd[k,j]
            # lower tri-diagonal block tri-diagonal elements
            i = addnonzero(ri, ri-K-1, -ukdjdkj, i)
            i = addnonzero(ri, ri-K, -ukjn[k,j], i)
            i = addnonzero(ri, ri-K+1, ukdjdkj, i)
            # diagonal block tri-diagonal elements
            i = addnonzero(ri, ri-1, -uknj[k,j], i)
            i = addnonzero(ri, ri, 1.0+ukj[k,j], i)
            i = addnonzero(ri, ri+1, -ukpj[k,j], i)
            # upper tri-diagonal block tri-diagonal elements
            i = addnonzero(ri, ri+K-1, -ukdjdkj, i)
            i = addnonzero(ri, ri+K, -ukjp[k,j], i)
            i = addnonzero(ri, ri+K+1, ukdjdkj, i)
        end
        # last row has 1 along diagonal
        i = addnonzero(j*K, j*K, 1.0, i)
    end
    # last block diagonal is identity
    const lastbase = (J-1)K
    for k = 1:K
        row[i+k-1], col[i+k-1], v[i+k-1] = lastbase + k, lastbase + k, 1.0
    end
    i += K - 1
    # return sparse L
    sparse(row[1:i], col[1:i], v[1:i], K*J, K*J)
end

# returns b vector for PDE2DSolver
#
# The given u is of size K x J, with x along rows, and y along columns.
#
# The returned b has KJ elements.
function _getb(p::PDE2DSolver, u::Matrix{Float64})
    const K, J = size(p)
    @assert K == size(u, 1) && J == size(u, 2)
    # fill b as b[k,j] and then reshape
    b = Array(Float64, K, J)
    # b[:,1] = u[:,1]
    # first block for j = 1 is different due to boundary conditions
    b[1,1] = u[1,1]
    for k = 2:(K-1)
        b[k,1] = (1.0 - 2.0p.ukjp[k,1])*u[k,1] + 2.0p.ukjp[k,1]*u[k,2]
    end
    b[K,1] = u[K,1]
    # for all j > 1
    for j = 2:(J-1)
        b[1,j] = u[1,j]
        for k = 2:(K-1)
            b[k,j] = (1.0-p.ukj[k,j])*u[k,j] +
                p.ukpj[k,j]*u[k+1,j] + p.uknj[k,j]*u[k-1,j] +
                p.ukjp[k,j]*u[k,j+1] + p.ukjn[k,j]*u[k,j-1] +
                p.ukdjd[k,j]*(u[k+1,j+1]+u[k-1,j-1]-u[k+1,j-1]-u[k-1,j+1])
        end
        b[K,j] = u[K,j]
    end
    b[:,J] = u[:,J]
    reshape(b, K*J)
end

# PDE solver using single matrix inversion over collapse 2D space
# This solver uses the ADI method and does not support <dx dy> != 0
immutable PDE2DADISolver <: AbstractPDE2DSolver
    ukjx::Matrix{Float64}
    ukjy::Matrix{Float64}
    ukpj::Matrix{Float64}
    uknj::Matrix{Float64}
    ukjp::Matrix{Float64}
    ukjn::Matrix{Float64}
    L1::Vector{Tridiagonal{Float64}}
    L2::Vector{Tridiagonal{Float64}}

    function PDE2DADISolver{T <: Real}(dx::Matrix{T}, dx2::Matrix{T},
        dy::Matrix{T}, dy2::Matrix{T}, Dx::Real, Dy::Real, dt::Real)
        Dx, Dy, dt = float(Dx), float(Dy), float(dt)
        Dx > 0.0 || error("Dx needs to be positive")
        Dy > 0.0 || error("Dy needs to be positive")
        dt > 0.0 || error("dt needs to be positive")
        const K, J = size(dx, 1), size(dx, 2)
        K == size(dy, 1) && J == size(dy, 2) || error("dx and dy need to be of same size")
        K == size(dx2, 1) && J == size(dx2, 2) || error("dx and dx2 need to be of same size")
        K == size(dy2, 1) && J == size(dy2, 2) || error("dx and dy2 need to be of same size")
        const dtDx, dtDy, dtDx2, dtDy2 = dt/Dx, dt/Dy, dt/(Dx*Dx), dt/(Dy*Dy)
        const ukjx = 0.5dtDx2*float(dx2)
        const ukjy = 0.5dtDy2*float(dy2)
        const ukpj = 0.25(dtDx*float(dx) .+ dtDx2*float(dx2))
        const uknj = 0.25(dtDx2*float(dx2) .- dtDx*float(dx))
        const ukjp = 0.25(dtDy*float(dy) .+ dtDy2*float(dy2))
        const ukjn = 0.25(dtDy2*float(dy2) .- dtDy*float(dy))
        const L1, L2 = _get2DADISolverL(ukjx, ukjy, ukpj, uknj, ukjp, ukjn)
        new(ukjx, ukjy, ukpj, uknj, ukjp, ukjn, L1, L2)
    end
end
size(p::PDE2DADISolver) = size(p.ukj)
function PDE2DADISolver(r::AbstractProcSpace, dt::Real)
    dt > zero(dt) || error("dt needs to be positive")
    const dx, dx2, dy, dy2, dxy, Dx, Dy = moments(r, dt)
    @assert all([isapprox(dxyi, 0.0) for dxyi in vec(dxy)])
    PDE2DADISolver(dx, dx2, dy, dy2, Dx, Dy, dt)
end

function expectedu{T <: Real}(p::PDE2DADISolver, u::Matrix{T})
    const K, J = size(u, 1), size(u, 2)
    @assert K == size(p.ukjx, 1) && J == size(p.ukjx, 2)
    # step 1 from n+1 to n+1/2
    u1 = Array(Float64, K, J)
    # boundary condition for j = 1
    u1[1,1] = u[1,1]
    for k = 2:(K-1)
        u1[k,1] = (1.0-2.0p.ukjp[k,1])*u[k,1] + 2.0p.ukjp[k,1]*u[k,2]
    end
    u1[K,1] = u[K,1]
    # 2 <= j <= J - 1
    bj = Array(Float64, K)
    for j = 2:(J-1)
        bj[1] = u[1,j]
        for k = 2:(K-1)
            @inbounds begin
                bj[k] = (1.0-p.ukjy[k,j])*u[k,j] +
                        p.ukjp[k,j]*u[k,j+1] + p.ukjn[k,j]*u[k,j-1]
            end
        end
        bj[K] = u[K,j]
        u1[:,j] = p.L1[j-1] \ bj
    end
    # boundary condition for j = J
    u1[:,J] = u[:,J]

    # step 2 from n+1/2 to n
    u2 = Array(Float64, K, J)
    # boundary condition for k = 1
    u2[1,:] = u1[1,:]
    # 2 <= k <= K - 1
    bk = Array(Float64, J)
    for k = 2:(K-1)
        bk[1] = u1[k,1]
        for j = 2:(J-1)
            @inbounds begin
                bk[j] = (1.0-p.ukjx[k,j])*u1[k,j] +
                        p.ukpj[k,j]*u1[k+1,j] + p.uknj[k,j]*u1[k-1,j]
            end
        end
        bk[J] = u1[k,J]
        u2[k,:] = p.L2[k-1] \ bk
    end
    # boundary condition for k = K
    u2[K,:] = u1[K,:]

    u2
end

function _get2DADISolverL(ukjx, ukjy, ukpj, uknj, ukjp, ukjn)
    const K, J = size(ukjx, 1), size(ukjx, 2)
    L1 = Array(Tridiagonal{Float64}, J-2)
    # L1 matrices for 2 <= j <= J-1
    for j = 2:(J-1)
        L1[j-1] = Tridiagonal([-uknj[2:(K-1),j], 0.0],
                              [1.0, 1.0.+ukjx[2:(K-1),j], 1.0],
                              [0.0, -ukpj[2:(K-1),j]])
    end
    L2 = Array(Tridiagonal{Float64}, K-2)
    # L2 matrices for 2 <= k <= K-1
    for k = 2:(K-1)
        L2[k-1] = Tridiagonal([-vec(ukjn[k,2:(J-1)]), 0.0],
                              [1.0+2.0ukjp[k,1], 1.0.+vec(ukjy[k,2:(J-1)]), 1.0],
                              [-2.0ukjp[k,1], -vec(ukjp[k,2:(J-1)])])
    end
    L1, L2
end

end # module


text:
# ------------------------------------------------------------------------------
# hyperelasticity example
# ------------------------------------------------------------------------------

using LinearAlgebra

# setup
mesh = Mesh3D(1.0, 1.0, 1.0)
finep = 2
coarsep = 1
numbercomponents = 3
dimension = 3
finebasis = TensorH1LagrangeBasis(finep + 1, finep + 1, numbercomponents, dimension)
coarsebasis = TensorH1LagrangeBasis(coarsep + 1, finep + 1, numbercomponents, dimension)
ctofbasis = TensorH1LagrangeBasis(
    coarsep + 1,
    finep + 1,
    numbercomponents,
    dimension,
    lagrangequadrature = true,
)

# constants
E = 1E6                             # Young's modulus
ν = 0.3                             # Poisson's ratio
K = E / (3 * (1 - 2 * ν))           # bulk modulus
λ = E * ν / ((1 + ν) * (1 - 2 * ν)) # Lamé parameters
μ = E / (2 * (1 + ν))

# state
gradu = [1; 2; 3] * ones(1, 3);

function neohookeanweakform(deltadu::Array{Float64}, w::Array{Float64})
    # dP = dF S + F dS

    # deformation gradient
    F = gradu + I
    J = det(F)
    # Green-Lagrange strain tensor
    E = (gradu * gradu' + gradu' * gradu) / 2
    # right Cauchy-Green tensor
    C = 2 * E + I
    C_inv = C^-1
    # second Piola-Kirchhoff
    S = λ * log(J) * C_inv + 2 * μ * C_inv * E

    # delta du
    deltadu = deltadu'
    # dF
    dF = deltadu + I
    # deltaE
    deltaE = (deltadu * deltadu' + deltadu' * deltadu) / 2
    # dS
    dS = λ * sum(C_inv .* deltaE) * C_inv + 2 * (μ - λ * log(J)) * C_inv * deltaE * C_inv
    # dP
    dP = (dF * S + F * dS) * w[1]

    return [dP']
end

# linearized Neo-Hookean operators
function makeoperator(basis::TensorBasis)
    inputs = [
        OperatorField(basis, [EvaluationMode.gradient], "gradent of deformation"),
        OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights"),
    ]
    outputs = [
        OperatorField(
            basis,
            [EvaluationMode.gradient],
            "test function gradient of deformation",
        ),
    ]
    return Operator(neohookeanweakform, mesh, inputs, outputs)
end
fineoperator = makeoperator(finebasis)
coarseoperator = makeoperator(coarsebasis)

# Chebyshev smoother
chebyshev = Chebyshev(fineoperator)

# p-multigrid preconditioner
multigrid = PMultigrid(fineoperator, coarseoperator, chebyshev, [ctofbasis])

# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π, π])
eigenvalues = real(eigvals(A))

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


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Source <: StaticInjection
        name::String
        available::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        X_th::Float64
        services::Vector{Service}
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

This struct acts as an infinity bus.

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `activepower::Float64`
- `reactivepower::Float64`
- `X_th::Float64`: Source Thevenin impedance, validation range: (0, nothing)
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Source <: StaticInjection
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    "Source Thevenin impedance"
    X_th::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Source(name, available, bus, activepower, reactivepower, X_th, services=Device[], ext=Dict{String, Any}(), )
    Source(name, available, bus, activepower, reactivepower, X_th, services, ext, InfrastructureSystemsInternal(), )
end

function Source(; name, available, bus, activepower, reactivepower, X_th, services=Device[], ext=Dict{String, Any}(), )
    Source(name, available, bus, activepower, reactivepower, X_th, services, ext, )
end

# Constructor for demo purposes; non-functional.
function Source(::Nothing)
    Source(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        X_th=0,
        services=Device[],
        ext=Dict{String, Any}(),
    )
end

InfrastructureSystems.get_name(value::Source) = value.name
"""Get Source available."""
get_available(value::Source) = value.available
"""Get Source bus."""
get_bus(value::Source) = value.bus
"""Get Source activepower."""
get_activepower(value::Source) = value.activepower
"""Get Source reactivepower."""
get_reactivepower(value::Source) = value.reactivepower
"""Get Source X_th."""
get_X_th(value::Source) = value.X_th
"""Get Source services."""
get_services(value::Source) = value.services
"""Get Source ext."""
get_ext(value::Source) = value.ext
"""Get Source internal."""
get_internal(value::Source) = value.internal

InfrastructureSystems.set_name!(value::Source, val::String) = value.name = val
"""Set Source available."""
set_available!(value::Source, val::Bool) = value.available = val
"""Set Source bus."""
set_bus!(value::Source, val::Bus) = value.bus = val
"""Set Source activepower."""
set_activepower!(value::Source, val::Float64) = value.activepower = val
"""Set Source reactivepower."""
set_reactivepower!(value::Source, val::Float64) = value.reactivepower = val
"""Set Source X_th."""
set_X_th!(value::Source, val::Float64) = value.X_th = val
"""Set Source services."""
set_services!(value::Source, val::Vector{Service}) = value.services = val
"""Set Source ext."""
set_ext!(value::Source, val::Dict{String, Any}) = value.ext = val
"""Set Source internal."""
set_internal!(value::Source, val::InfrastructureSystemsInternal) = value.internal = val


text:
using Plots
using FFTW
using NonlinearSchrodinger

function T1A_K(ψᵢ, b)
    γ = 1
    p = 1
    q = 1

    n = 1
    μ = 0

    #ψ = ifft(im*b.ω.*fftshift(fft(ψᵢ)))
    #ψₒ = exp.(-γ*ψᵢ.^(n-1).*b.dx.*ψ).*ψᵢ
    #ψₒ = fftshift(fft(ψₒ))
    #ψₒ = ifft(exp.(-b.dx*(im*b.ω*μ .+ im*b.ω.^3*p .+ im*b.ω.^5*q)).*ψₒ)

    ψₒ = fftshift(fft(ψᵢ))
    ψₒ = ifft(exp.(-b.dx*(im*b.ω*μ .+ im*b.ω.^3*p .+ im*b.ω.^5*q)).*ψₒ)
    ψₒ = exp.(-γ*ψᵢ.^(n-1).*b.dx.*ifft(+im*b.ω.*fftshift(fft(ψₒ)))).*ψₒ
end

xᵣ = 0=>50
box = Box(xᵣ, 200, dx=5e-2, Nₜ = 512, n_periods =3)

ψ = zeros(box.Nₓ, box.Nₜ) .+ im

ψ[1, :] = 105/169*sech.(1/2/sqrt(13)*box.t).^4
for i = 1:box.Nₓ-1
   ψ[i+1, :] .=  T1A_K(ψ[i,:], box)
end
heatmap(abs.(ψ))


text:
#!/usr/bin/env julia
#
# @license Apache-2.0
#
# Copyright (c) 2018 The Stdlib Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import JSON

"""
    gen( x, name )

Generate fixture data and write to file.

# Arguments

* `x`: domain
* `name::AbstractString`: output filename

# Examples

``` julia
julia> x = linspace( -1000, 1000, 2001 );
julia> gen( x, \"data.json\" );
```
"""
function gen( x, name )
	y = Array( Float64, length( x ) );
	for i in eachindex(x)
		y[i] = cbrt( x[i] );
	end

	# Store data to be written to file as a collection:
	data = Dict([
		("x", x),
		("expected", y)
	]);

	# Based on the script directory, create an output filepath:
	filepath = joinpath( dir, name );

	# Write the data to the output filepath as JSON:
	outfile = open( filepath, "w" );
	write( outfile, JSON.json(data) );
	close( outfile );
end

# Get the filename:
file = @__FILE__;

# Extract the directory in which this file resides:
dir = dirname( file );

# Subnormal values:
x = linspace( 1e-309, 1e-324, 503 );
gen( x, "subnormal.json" );

# Negative tiny values:
x = linspace( -1e-300, -1e-308, 503 );
gen( x, "tiny_negative.json" );

# Positive tiny values:
x = linspace( 1e-300, 1e-308, 503 );
gen( x, "tiny_positive.json" );

# Small(er) values:
x = linspace( -0.8, 0.8, 503 );
gen( x, "smaller.json" );

# Negative small values:
x = linspace( -0.8, -3.0, 503 );
gen( x, "small_negative.json" );

# Positive small values:
x = linspace( 0.8, 3.0, 503 );
gen( x, "small_positive.json" );

# Negative medium values:
x = linspace( -3.0, -20.0, 503 );
gen( x, "medium_negative.json" );

# Positive medium values:
x = linspace( 3.0, 20.0, 503 );
gen( x, "medium_positive.json" );

# Large negative values:
x = linspace( -50, -20, 503 );
gen( x, "large_negative.json" );

# Large positive values:
x = linspace( 20, 50, 503 );
gen( x, "large_positive.json" );

# Very large negative values:
x = linspace( -50, -500, 503 );
gen( x, "very_large_negative.json" );

# Very large positive values:
x = linspace( 50, 500, 503 );
gen( x, "very_large_positive.json" );

# Huge negative values:
x = linspace( -1e200, -1e208, 503 );
gen( x, "huge_negative.json" );

# Huge positive values:
x = linspace( 1e300, 1e308, 503 );
gen( x, "huge_positive.json" );

max_stars_repo_name: andbmme/stdlib

max_issues_repo_name: andbmme/stdlib

max_forks_repo_name: andbmme/stdlib


text:
# This file contains methods to generate a data set of instances
include("io.jl")

"""
Generate an n*n grid

Argument
- n: size of the grid
"""
function generateInstance(n::Int64)
    if n<=3
        println("size should be greater than 3")
        return
    end
	b=n # number of black cells
    directions=[1 0;-1 0;0 1;0 -1;1 1;1 -1;-1 -1;-1 1]
    grid=trues(n,n)
    t=zeros(Int,n,n)
    # put black cells and keep the white cells continuous
    for r=1:n
        canPutBlackCell=false
        while !canPutBlackCell
            c=rand(1:n)
            canPutBlackCell=true
            for i=1:size(directions,1)
                if r+directions[i,1]>=1 && r+directions[i,1]<=n && c+directions[i,2]>=1 && c+directions[i,2]<=n && grid[r+directions[i,1],c+directions[i,2]]==0
                    canPutBlackCell=false
                end
            end
            if canPutBlackCell==true
                grid[r,c]=0
            end
        end
    end
    for r=1:n
        canPutNumberCell=false
        while !canPutNumberCell
            c=rand(1:n)
            if grid[r,c]!=0
                visNum=visiable(grid,r,c)
                t[r,c]=visNum
                canPutNumberCell=true
            end
        end
    end
    # println(grid)
    return t
end
"""
count the visiable cells from position (l,c)
"""
function visiable(x_vals::BitArray{2},l::Int,c::Int)
    rowNum=size(x_vals,1)
    colNum=size(x_vals,2)

    up=0 # the number of white cells in top direction until meet black cell
    right=0 # the number of white cells in right direction until meet black cell
    down=0 # the number of white cells in down direction until meet black cell
    left=0 # the number of white cells in left direction until meet black cell
    for i1 in 1:l-1
        if  x_vals[l-i1,c]== 0  #meet the black cell
            down=i1-1
            break
        else
            down=l-1
        end
    end
    for j1 in 1:rowNum-l
        if  x_vals[l+j1,c]== 0
            up=j1-1
            break
        else
            up=rowNum-l
        end
    end
    for i2 in 1:c-1
        if  x_vals[l,c-i2]== 0
            left=i2-1
            break
        else
            left=c-1
        end
    end
    for j2 in 1:colNum-c
        if  x_vals[l,c+j2]== 0
            right=j2-1
            break
        else
            right=colNum-c
        end
    end
    return up+down+left+right+1
end
"""
Generate all the instances

Remark: a grid is generated only if the corresponding output file does not already exist
"""
function generateDataSet()
	# different grid size
	for n in [4,5,6]
		for m in [4,5,6]
			#generate 2 instances for each size
			for instance=1:10
				fileName="../data/instance_$(n)_$(m)_$(instance).txt"
				if !isfile(fileName)
					println("-- Generating file " * fileName)
	                saveInstance(generateInstance(n), fileName)
	            end
	        end
	    end
    end
end


text:
using Printf
# the header implemented based on LAS 1.2
# TODO: check compatibility other LAS versions

#=
COMPATIBILITY WITH LAS 1.2:
One unavoidable change has been made to the Public Header Block; Start of Waveform
Data Packet Record. This long, long has been added to the end of the block and thus
little or no change will be needed in LAS 1.2 readers that do not need waveform data.
There are no changes to Point Data Record types 0 through 3. The waveform encoded
data types have been added as Point Data Record types 4 and 5.

THE ADDITIONS OF LAS 1.4 INCLUDE:
Backward compatibility with LAS 1.1 – LAS 1.3 when payloads consist of only legacy
content
=#

mutable struct LasHeader
    file_source_id::UInt16
    global_encoding::UInt16
    guid_1::UInt32
    guid_2::UInt16
    guid_3::UInt16
    guid_4::AbstractString
    version_major::UInt8
    version_minor::UInt8
    system_id::AbstractString
    software_id::AbstractString
    creation_dayofyear::UInt16
    creation_year::UInt16
    header_size::UInt16
    data_offset::UInt32
    n_vlr::UInt32
    data_format_id::UInt8
    data_record_length::UInt16
    records_count::UInt32
    point_return_count::Vector{UInt32}
    x_scale::Float64
    y_scale::Float64
    z_scale::Float64
    x_offset::Float64
    y_offset::Float64
    z_offset::Float64
    x_max::Float64
    x_min::Float64
    y_max::Float64
    y_min::Float64
    z_max::Float64
    z_min::Float64
    start_of_waveform_data_packet_record::UInt64
    start_of_first_extended_variable_length_record::UInt64
    number_of_extended_variable_length_records::UInt32
    extended_number_of_point_records::UInt64
    extended_number_of_points_by_return::Vector{UInt64}
    variable_length_records::Vector{LasVariableLengthRecord}
    user_defined_bytes::Vector{UInt8}
end

function Base.getproperty(h::LasHeader, s::Symbol)
    if s == :creation_doy
        Base.depwarn("LasHeader.creation_doy is deprecated, use LasHeader.creation_dayofyear", :getproperty)
        return getfield(h, :creation_dayofyear)
    else
        return getfield(h, s)
    end
end

function Base.show(io::IO, header::LasHeader)
    n = header.extended_number_of_point_records > 0 ? Int(header.extended_number_of_point_records) : Int(header.records_count)
    println(io, "LasHeader with $(n) points.")
    println(io, string("\tfile_source_id = ", header.file_source_id))
    println(io, string("\tglobal_encoding = ", header.global_encoding))
    println(io, string("\tguid_1 = ", header.guid_1))
    println(io, string("\tguid_2 = ", header.guid_2))
    println(io, string("\tguid_3 = ", header.guid_3))
    println(io, string("\tguid_4 = ", header.guid_4))
    println(io, string("\tversion_major = ", header.version_major))
    println(io, string("\tversion_minor = ", header.version_minor))
    println(io, string("\tsystem_id = ", header.system_id))
    println(io, string("\tsoftware_id = ", header.software_id))
    println(io, string("\tcreation_dayofyear = ", header.creation_dayofyear))
    println(io, string("\tcreation_year = ", header.creation_year))
    println(io, string("\theader_size = ", header.header_size))
    println(io, string("\tdata_offset = ", header.data_offset))
    println(io, string("\tn_vlr = ", header.n_vlr))
    println(io, string("\tdata_format_id = ", header.data_format_id))
    println(io, string("\tdata_record_length = ", header.data_record_length))
    println(io, string("\trecords_count = ", header.records_count))
    println(io, string("\tpoint_return_count = ", [Int(i) for i in header.point_return_count]))
    println(io, string("\tx_scale = ", header.x_scale))
    println(io, string("\ty_scale = ", header.y_scale))
    println(io, string("\tz_scale = ", header.z_scale))
    println(io, string("\tx_offset = ", header.x_offset))
    println(io, string("\ty_offset = ", header.y_offset))
    println(io, string("\tz_offset = ", header.z_offset))
    println(io, @sprintf "\tx_max = %.7f" header.x_max)
    println(io, @sprintf "\tx_min = %.7f" header.x_min)
    println(io, @sprintf "\ty_max = %.7f" header.y_max)
    println(io, @sprintf "\ty_min = %.7f" header.y_min)
    println(io, @sprintf "\tz_max = %.7f" header.z_max)
    println(io, @sprintf "\tz_min = %.7f" header.z_min)

    if header.version_minor > 3
        println(io, string("\tnumber_of_extended_variable_length_records = ", header.number_of_extended_variable_length_records))
        println(io, string("\textended_number_of_point_records = ", header.extended_number_of_point_records))
        println(io, string("\textended_number_of_points_by_return = ", [Int(i) for i in header.extended_number_of_points_by_return]))
    end
    if !isempty(header.variable_length_records)
        nrecords = min(10, size(header.variable_length_records, 1))

        println(io, string("\tvariable_length_records (showing first $(nrecords)/$(size(header.variable_length_records, 1))) = "))
        for vlr in header.variable_length_records[1:nrecords]
            println(io, "\t\t($(vlr.user_id), $(vlr.record_id)) => ($(vlr.description), $(sizeof(vlr.data)) bytes...)")
        end
        println("\t\t...")
    end
end

function readstring(io, nb::Integer)
    bytes = read(io, nb)
    # strip possible null bytes
    lastchar = findlast(bytes .!= 0)
    if lastchar == nothing
        return ""
    else
        return String(bytes[1:lastchar])
    end
end

function writestring(io, str::AbstractString, nb::Integer)
    n = length(str)
    npad = nb - n
    if npad < 0
        error("string too long")
    elseif npad == 0
        write(io, str)
    else
        writestr = string(str * "\0"^npad)
        write(io, writestr)
    end
end

function Base.read(io::IO, ::Type{LasHeader})
    file_source_id = read(io, UInt16)
    global_encoding = read(io, UInt16)
    guid_1 = read(io, UInt32)
    guid_2 = read(io, UInt16)
    guid_3 = read(io, UInt16)
    guid_4 = readstring(io, 8)
    version_major = read(io, UInt8)
    version_minor = read(io, UInt8)
    system_id = readstring(io, 32)
    software_id = readstring(io, 32)
    creation_dayofyear = read(io, UInt16)
    creation_year = read(io, UInt16)
    header_size = read(io, UInt16)
    data_offset = read(io, UInt32)
    n_vlr = read(io, UInt32)
    data_format_id = read(io, UInt8)
    data_record_length = read(io, UInt16)
    records_count = read(io, UInt32)
    point_return_count = read!(io, Vector{UInt32}(undef, 5))
    x_scale = read(io, Float64)
    y_scale = read(io, Float64)
    z_scale = read(io, Float64)
    x_offset = read(io, Float64)
    y_offset = read(io, Float64)
    z_offset = read(io, Float64)
    x_max = read(io, Float64)
    x_min = read(io, Float64)
    y_max = read(io, Float64)
    y_min = read(io, Float64)
    z_max = read(io, Float64)
    z_min = read(io, Float64)
    lasversion = VersionNumber(version_major, version_minor)
    start_of_waveform_data_packet_record = UInt64(0)
    start_of_first_extended_variable_length_record = UInt64(0)
    number_of_extended_variable_length_records = UInt32(0)
    extended_number_of_point_records = UInt64(0)
    extended_number_of_points_by_return = fill(UInt32(0), 15)

    if lasversion >= v"1.3"
        # start of waveform data record (unsupported)
        start_of_waveform_data_packet_record = read(io, UInt64)
        start_of_first_extended_variable_length_record = read(io, UInt64)
        number_of_extended_variable_length_records = read(io, UInt32)
        extended_number_of_point_records = read(io, UInt64) # this is a UInt64 which is different from records_count which is UInt32
        extended_number_of_points_by_return = read!(io, Vector{UInt64}(undef, 15))
    end

    vlrs = [read(io, LasVariableLengthRecord, false) for i=1:n_vlr]

    # From here until the data_offset everything is read in
    # as user_defined_bytes. To avoid a seek that we cannot do on STDIN,
    # we calculate how much to read in.
    vlrlength = n_vlr == 0 ? 0 : sum(sizeof, vlrs)
    pos = header_size + vlrlength
    user_defined_bytes = read(io, data_offset - pos)

    # put it all in a type
    header = LasHeader(
        file_source_id,
        global_encoding,
        guid_1,
        guid_2,
        guid_3,
        guid_4,
        version_major,
        version_minor,
        system_id,
        software_id,
        creation_dayofyear,
        creation_year,
        header_size,
        data_offset,
        n_vlr,
        data_format_id,
        data_record_length,
        records_count,
        point_return_count,
        x_scale,
        y_scale,
        z_scale,
        x_offset,
        y_offset,
        z_offset,
        x_max,
        x_min,
        y_max,
        y_min,
        z_max,
        z_min,
        start_of_waveform_data_packet_record,
        start_of_first_extended_variable_length_record,
        number_of_extended_variable_length_records,
        extended_number_of_point_records,
        extended_number_of_points_by_return,
        vlrs,
        user_defined_bytes,
    )
end

function Base.write(io::IO, h::LasHeader)
    write(io, h.file_source_id)
    write(io, h.global_encoding)
    write(io, h.guid_1)
    write(io, h.guid_2)
    write(io, h.guid_3)
    writestring(io, h.guid_4, 8)
    write(io, h.version_major)
    write(io, h.version_minor)
    writestring(io, h.system_id, 32)
    writestring(io, h.software_id, 32)
    write(io, h.creation_dayofyear)
    write(io, h.creation_year)
    write(io, h.header_size)
    write(io, h.data_offset)
    @assert length(h.variable_length_records) == h.n_vlr
    write(io, h.n_vlr)
    write(io, h.data_format_id)
    write(io, h.data_record_length)
    write(io, h.records_count)
    @assert length(h.point_return_count) == 5
    write(io, h.point_return_count)
    write(io, h.x_scale)
    write(io, h.y_scale)
    write(io, h.z_scale)
    write(io, h.x_offset)
    write(io, h.y_offset)
    write(io, h.z_offset)
    write(io, h.x_max)
    write(io, h.x_min)
    write(io, h.y_max)
    write(io, h.y_min)
    write(io, h.z_max)
    write(io, h.z_min)
    lasversion = VersionNumber(h.version_major, h.version_minor)
    if lasversion >= v"1.3"
        # start of waveform data record (unsupported)
        write(io, h.start_of_waveform_data_packet_record) # start_of_waveform_data_packet_record
        write(io, h.start_of_first_extended_variable_length_record) # start_of_first_extended_variable_length_record
        write(io, h.number_of_extended_variable_length_records) # number_of_extended_variable_length_records
        write(io, h.extended_number_of_point_records)  # extended_number_of_point_records
        write(io, h.extended_number_of_points_by_return) # extended_number_of_points_by_return
    end
    for i in 1:h.n_vlr
        write(io, h.variable_length_records[i])
    end
    write(io, h.user_defined_bytes)
    # note that for LAS 1.4 a few new parts need to be written
    # possibly introduce typed headers like the points
    nothing
end

"""If true, GPS Time is standard GPS Time (satellite GPS Time) minus 1e9.
If false, GPS Time is GPS Week Time.

Note that not all software sets this encoding correctly."""
is_standard_gps(h::LasHeader) = isodd(h.global_encoding)

"Check if the projection information is in WKT format (true) or GeoTIFF (false)"
function is_wkt(h::LasHeader)
    wkit_bit = Bool((h.global_encoding & 0x0010) >> 4)
    if !wkit_bit && h.data_format_id > 5
        throw(DomainError("WKT bit must be true for point types higher than 5"))
    end
    wkit_bit
end

"Retrieve the bounding box from a LasHeader as a NamedTuple"
boundingbox(h::LasHeader) = (xmin = h.x_min, ymin = h.y_min, zmin = h.z_min,
    xmax = h.x_max, ymax = h.y_max, zmax = h.z_max)

"Get the boundingbox in the same stored integer form as the raw point coordinates"
boundingbox_unscaled(h::LasHeader) = (
    xmin = round(Int32, (h.x_min - h.x_offset) / h.x_scale),
    ymin = round(Int32, (h.y_min - h.y_offset) / h.y_scale),
    zmin = round(Int32, (h.z_min - h.z_offset) / h.z_scale),
    xmax = round(Int32, (h.x_max - h.x_offset) / h.x_scale),
    ymax = round(Int32, (h.y_max - h.y_offset) / h.y_scale),
    zmax = round(Int32, (h.z_max - h.z_offset) / h.z_scale)
)


text:
module GoogleLanguageService

using Logging
using HTTP
include("language_service_pb.jl")

export Document, AnalyzeEntitiesRequest, AnalyzeEntitiesResponse, Document_Type, API_KEY, analyze_entities

API_KEY = get(ENV, "API_KEY", nothing)

function analyze_entities(text::String)

    if API_KEY == nothing
        throw("API_KEY not found")
    else
        @info "Using API_KEY from current environment"
    end

    doc = Document(
        _type = Document_Type.PLAIN_TEXT,
        language = "en",
        content = text
    )

    # analyse
    ent = AnalyzeEntitiesRequest(
        document = doc
    )

    # create bytebuffer
    iob = IOBuffer()
    writeproto(iob, ent)
    write_proto_bytes = take!(iob)

    # Make API call
    headers = Dict(
        "X-Goog-Api-Key" => API_KEY,
        "User-Agent" => "myapp/0.1",
        "Content-Type" => "application/x-protobuf",
    )

    @info "Making HTTP request to gRPC endpoint"
    data_response_bytes = HTTP.request(
        "POST",
        "https://language.googleapis.com/\$rpc/google.cloud.language.v1.LanguageService/AnalyzeEntities",
        headers,
        write_proto_bytes,
    ).body # this is byte array

    api_response = AnalyzeEntitiesResponse()
    response_buffer = IOBuffer(data_response_bytes)
    readproto(IOBuffer(data_response_bytes), api_response)

    return api_response
end

end # module

max_stars_repo_name: sivakon/grpc-julia-gcp-example

max_issues_repo_name: sivakon/grpc-julia-gcp-example

max_forks_repo_name: sivakon/grpc-julia-gcp-example


text:
using Crux, Flux, POMDPs, POMDPGym, Distributions, POMDPPolicies, BSON
init_mujoco_render() # Required for visualization

# Construct the Mujoco environment
mdp = GymPOMDP(:HalfCheetah, version = :v3)
S = state_space(mdp)
adim = length(POMDPs.actions(mdp)[1])
amin = -1*ones(Float32, adim)
amax = 1*ones(Float32, adim)
rand_policy = FunctionPolicy((s) -> Float32.(rand.(Uniform.(amin, amax))))

## (Load from bson) get expert trajectories
expert_trajectories = BSON.load("examples/il/expert_data/half_cheetah_mujoco.bson")[:data]
μ_s = mean(expert_trajectories[:s], dims=2)[:]
σ_s = std(expert_trajectories[:s], dims=2)[:] .+ 1f-3
S = state_space(mdp, μ=μ_s, σ=σ_s)

# Initializations that match the default PyTorch initializations
Winit(out, in) = Float32.(rand(Uniform(-sqrt(1/in), sqrt(1/in)), out, in))
binit(in) = (out) -> Float32.(rand(Uniform(-sqrt(1/in), sqrt(1/in)), out))

## Behavioral Cloning
μ() = ContinuousNetwork(Chain(Dense(S.dims[1], 64, tanh), Dense(64, 32, tanh), Dense(32, adim)))
log_std() = -0.5f0*ones(Float32, adim)

𝒮_bc = BC(π=μ(),
          𝒟_demo=expert_trajectories,
          S=S,
          opt=(epochs=100000, batch_size=1024),
          log=(period=500,),
          max_steps=1000)
solve(𝒮_bc, mdp)

## AdVIL
A() = ContinuousNetwork(Chain(Dense(S.dims[1], 256, relu, init=Flux.orthogonal),
            Dense(256, 256, relu, init=Flux.orthogonal),
            Dense(256, 6, init=Flux.orthogonal)))
D() = ContinuousNetwork(Chain(Dense(S.dims[1] + adim, 256, relu, init=Flux.orthogonal),
            Dense(256, 256, relu, init=Flux.orthogonal),
            Dense(256, 1, init=Flux.orthogonal)))

𝒮_advil = AdVIL(π=ActorCritic(A(),D()), 𝒟_demo=expert_trajectories, S=S, a_opt=(epochs=100000, optimizer=ADAM(8f-6), batch_size=1024), c_opt=(optimizer=ADAM(8e-4),), max_steps=1000)
solve(𝒮_advil, mdp)


text:
# # fast algorithm from https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon
#
# struct Pt
#     x::Float64
#     y::Float64
# end
#
# # takes two points and generates a function giving the line. All points giving 0 from that function will be on the line
# _points_to_line(pt1::Pt, pt2::Pt) = (pt::Pt) -> (pt2.y - pt1.y)*pt.x + (pt2.x - pt1.x)*pt.y + (pt2.x*pt1.y - pt1.x*pt2.y)
# function _lineintersects(pt1_1, pt1_2, pt2_1, pt2_2,
#                         linefun1 = _points_to_line(pt1_1, pt1_2), linefun2 = _points_to_line(pt2_1, pt2_2))
#     (pt1_1 == pt1_2 || pt2_1 == pt2_2) && return 0 # line has length 0
#     tmp1, tmp2 = linefun1(pt2_1), linefun1(pt2_2)
#     sign(tmp1) == sign(tmp2) && return 0
#     tmp3, tmp4 = linefun2(pt1_1), linefun2(pt1_2)
#     sign(tmp3) == sign(tmp4) && return 0
#     return 1
# end
#
#
# function point_in_polygon(pt::Pt, polygon::AbstractVector{<:Pt},
#                         bbox = (extrema(p.x for p in polygon)..., extrema(p.y for p in polygon)...))
#
#     (pt.x < bbox[1] || pt.x > bbox[2] || pt.y < bbox[3] || pt.y > bbox[4] || bbox[1] == bbox[2] || bbox[3] == bbox[4]) && return false
#
#     raypt = Pt(bbox[1] - (bbox[2] - bbox[1]), pt.y)
#     rayline = _points_to_line(raypt, pt)
#     intersections = 0
#
#     @inbounds for i in eachindex(polygon)[2:end]
#         intersections += _lineintersects(raypt, pt, polygon[i-1], polygon[i], rayline)
#     end
#     intersections += _lineintersects(raypt, pt, last(polygon), first(polygon), rayline)
#     Bool(intersections%2)
# end


text:
function draw(lcmgl::LCMGLClient, mesh::AbstractMesh; face_color=[.7, .7, .2, .4], draw_normals=false)
    color(lcmgl, face_color...)
    begin_mode(lcmgl, LCMGL.TRIANGLES)
    verts = vertices(mesh)
    for face in faces(mesh)
        for point in verts[face][end:-1:1]
            vertex(lcmgl, point[1], point[2], point[3])
        end
    end
    end_mode(lcmgl)
    color(lcmgl, 0., 0., 0., 1.0)
    begin_mode(lcmgl, LCMGL.LINES)
    for face in faces(mesh)
        for i = 1:length(verts[face])
            vertex(lcmgl, verts[face][i]...)
            j = i + 1
            if i == length(verts[face])
                j = 1
            end
            vertex(lcmgl, verts[face][j]...)
        end
    end
    if draw_normals
      normals = decompose(Normal{3, Float64}, mesh)
      for i = 1:length(verts)
        vertex(lcmgl, verts[i]...)
        vertex(lcmgl, (verts[i] + Point{3, Float64}(normals[i]) * 0.05)...)
      end
    end
    end_mode(lcmgl)
end

function draw(mesh::AbstractMesh)
    LCMGLClient("soft_robot") do lcmgl
      draw(lcmgl, mesh)
      switch_buffer(lcmgl)
    end
end

function draw(robot::SoftRobot, state::SoftRobotState)
  LCMGLClient("soft_robot") do lcmgl
    draw(lcmgl, robot, state)
    switch_buffer(lcmgl)
  end
end

function draw(lcmgl::LCMGLClient, robot::SoftRobot, state::SoftRobotState)
  body_mesh = HomogenousMesh(state.positions, robot.faces)
  draw(lcmgl, body_mesh)

  LCMGLClient("$(lcmgl.name)_barrier") do lcmgl
    draw(lcmgl, barrier_mesh(robot, state), face_color=[.8, .1, .1, .3])
    switch_buffer(lcmgl)
  end
end


text:
const TimeStamp = typeof(time_ns())

struct Metadata <: Object
    source::LineNumberNode
    _module::Module
    location::UUID
    threadid::Int
    timestamp::TimeStamp
end

Metadata(source, _module, location = UUID(0)) =
    Metadata(source, _module, location, Threads.threadid(), time_ns())

metadata_expr(__source__, __module__) =
    :($Metadata($(QuoteNode(__source__)), $(QuoteNode(__module__)), $(QuoteNode(uuid4()))))

function _summary(io, metadata::Metadata)
    dt = roughly_as_datetime(metadata.timestamp)
    file = string(metadata.source.file)
    line = metadata.source.line
    print(io, metadata._module, " ", file, ":", line, " @", dt)
end

function Base.show(io::IO, ::MIME"text/plain", metadata::Metadata)
    dt = roughly_as_datetime(metadata.timestamp)
    file = string(metadata.source.file)
    line = metadata.source.line
    print(io, "<")
    _summary(io, metadata)
    print(io, ">")
end

roughly_as_datetime(t::TimeStamp) = now() + Second(floor(Int, (time_ns() - t) * 1e-9))


text:
# broken code i need to reference

function p2(lines)
  t = 0 #p1 ans
  scores = [] #p2
  for (i,line) in enumerate(lines)
    cor = false
    Q = Queue{Char}()
    for c in collect(line)
      # println("$c")
      if c ∈ ['(', '[', '{', '<']
        enqueue!(Q, c)
        # @show Q
      elseif c == ')'
        if isempty(Q) || last(Q) != '('
          t += 3
          cor = true
          break
        else
          !isempty(Q) && dequeue!(Q)
        end
      elseif c == ']'
      if isempty(Q) || last(Q) != '['
          t += 57
          cor = true
          break
        else
          !isempty(Q) && dequeue!(Q)
        end
      elseif c == '}'
        if isempty(Q) || last(Q) != '{'
          t += 1197
          cor = true
          break
        else
          !isempty(Q) && dequeue!(Q)
        end
      elseif c == '>'
        if isempty(Q) || last(Q) != '<'
          t += 25137
          cor = true
          break
        else
          !isempty(Q) && dequeue!(Q)
        end
      end
    end
    if !cor
      score  = 0
      map = Dict('(' => 1, '[' => 2, '{' => 3, '<' => 4)
      while !isempty(Q)
        val = dequeue!(Q)
        score *= 5
        score += map[val]
      end
      push!(scores, score)
    end
  end
  # @show t
end


text:
printstyled("  dataless SEED\n", color=:light_green)
redirect_stdout(out) do
  metafile = path*"/SampleFiles/SEED/jones.hood.dataless"
  S = read_meta("dataless", metafile, v=3,
                s="2008-01-01T00:00:00",
                t="2008-02-01T00:00:00",
                units=true)
  S2 = read_dataless( metafile, v=3,
                      s=DateTime("2008-01-01T00:00:00"),
                      t=DateTime("2008-02-01T00:00:00"),
                      units=true)
  @test S == S2
  files = ls(path*"/SampleFiles/SEED/*.dataless")
  for i in files
    println("Reading file ", i)
    S = read_meta("dataless", i, v=0, units=false)
    S = read_meta("dataless", i, v=1, units=false)
    S = read_meta("dataless", i, v=2, units=false)
    S = read_meta("dataless", i, v=3, units=false)
    S = read_meta("dataless", i, v=3, units=true)
  end
end


text:
"""
    CommandException <: Exception

Abstract type for command exceptions. All command exceptions
should contain an `exitcode` field.
"""
abstract type CommandException <: Exception end

struct CommandExit <: CommandException
    exitcode::Int
end

struct CommandError <: CommandException
    msg::String
    exitcode::Int
end

cmd_error(msg::String, code::Int = 1) = throw(CommandError(msg, code))
cmd_exit(code::Int = 0) = throw(CommandExit(code))

function Base.showerror(io::IO, e::CommandError)
    print(io, "CommandError: ", e.msg, " exit with ", e.exitcode)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.1
#   kernelspec:
#     display_name: Julia 1.1.0
#     language: julia
#     name: julia-1.1
# ---

# # Spatial diffusion / smoothing example
#
# The unit testing of `MeshArrays.jl` uses the `smooth()` function. Its effect on a random noise field is to "smooth it out". The average smoother efficiency is predictable and can be set via a `scale parameter`. The method as defined inside `MeshArrays.demo2` is applicable for any regular, `cube-sphere`, or `lat-lon-cap` Earth grid.
#
# For more about how this diffusion based method works and behaves, please refer to **Weaver and Courtier, 2001.** Correlation modelling on the sphere using a generalized diffusion equation. https://doi.org/10.1002/qj.49712757518

# #### Load  `MeshArrays` and `Plots`  modules

using MeshArrays, Plots
include(joinpath(dirname(pathof(MeshArrays)),"Plots.jl"))

# Grid the `6` faces of a cube with `16*16` points on each face. Distances, areas, etc. are all set to `1.0` for simplicity.

GridVariables=GridOfOnes("cs",6,16);

# #### Diffuse a noisy field over the surface of a cube

DemoVariables=MeshArrays.demo2(GridVariables);

# #### Display the 6 cube faces using  `heatmap`
#
# The orignal noise field is `DemoVariables[1]` while the smoothed one
# is `DemoVariables[2]`.

heatmap(DemoVariables[1],title="initial noise",clims=(-0.5,0.5))

# After `smooth()` has been applied via `demo2()`, the noise is smoother and muted.

heatmap(DemoVariables[2],title="smoothed noise",clims=(-0.5,0.5))

# The cost of applying `smooth()` predictably  increases with the decorrelation scale. For more about how this works, please refer to **Weaver and Courtier, 2001.** Correlation modelling on the sphere using a generalized diffusion equation. https://doi.org/10.1002/qj.49712757518

Rini=DemoVariables[1]
DXCsm=DemoVariables[3]
DYCsm=DemoVariables[4]
@time Rend=smooth(Rini,DXCsm,DYCsm,GridVariables);
@time Rend=smooth(Rini,2DXCsm,2DYCsm,GridVariables);

max_stars_repo_name: lmilechin/MeshArrayNotebooks

max_issues_repo_name: lmilechin/MeshArrayNotebooks

max_forks_repo_name: lmilechin/MeshArrayNotebooks


text:
### LINEAR MAP CONSTRUCTION

# for a given function of function-like object A, which acts upon data of type `input`
# and returns data of type `output`
# return a LinearMap that acts upon a vector form of u
linear_map(A,input,output;eltype=Float64) = _linear_map(A,input,output,eltype,Val(length(input)),Val(length(output)))

linear_map(A,input;eltype=Float64) = _linear_map(A,input,eltype,Val(length(input)))

linear_map(A::AbstractMatrix{T},input::AbstractVector{T};eltype=Float64) where {T} = LinearMap{eltype}(A)

linear_map(A::SaddleSystem{T,Ns,Nc},::Any;eltype=Float64) where {T,Ns,Nc} = LinearMap{T}(x->A*x,Ns+Nc)

linear_map(A::UniformScaling,input;eltype=Float64) = LinearMap{eltype}(x->A*x,length(input))

function linear_inverse_map(A,input;eltype=Float64)
    hasmethod(\,Tuple{typeof(A),typeof(input)}) || error("No such backslash operator exists")
    return LinearMap{eltype}(_create_vec_backslash(A,input),length(input))
end
linear_inverse_map(A::SaddleSystem{T,Ns,Nc},::Any;eltype=Float64) where {T,Ns,Nc} = LinearMap{T}(x->A\x,Ns+Nc)

# Square operators. input of zero length
_linear_map(A,input,eltype,::Val{0}) =
      LinearMap{eltype}(x -> (),0,0)

# Square operators. input of non-zero length
_linear_map(A,input,eltype,::Val{M}) where {M} =
      LinearMap{eltype}(_create_fcn(A,input),length(input))

# input and output have zero lengths
_linear_map(A,input,output,eltype,::Val{0},::Val{0}) =
      LinearMap{eltype}(x -> (),0,0)

# input is 0 length, output is not
_linear_map(A,input,output,eltype,::Val{0},::Val{M}) where {M} =
      LinearMap{eltype}(x -> _unwrap_vec(zero(output)),length(output),0)

# output is 0 length, input is not
_linear_map(A,input,output,eltype,::Val{N},::Val{0}) where {N} =
      LinearMap{eltype}(x -> (),0,length(input))

# non-zero lengths of input and output
_linear_map(A,input,output,eltype,::Val{N},::Val{M}) where {N,M} =
      LinearMap{eltype}(_create_fcn(A,input),length(output),length(input))

# Create a function for operator A that can act upon an input of type AbstractVector
# and return an output of type AbstractVector. It should wrap the input vector
# in the input data type associated with A and it should then unwrap its
# output back into vector form
function _create_fcn(A,input)
    # if A has an associated * operation, then use this
    if hasmethod(*,Tuple{typeof(A),typeof(input)})
        fcn = _create_vec_multiplication(A,input)
    # or if A is a function or function-like object, then use this
    elseif hasmethod(A,Tuple{typeof(input)})
        fcn = _create_vec_function(A,input)
    # or just quit
    else
        error("No function exists for this operator to act upon this type of data")
    end
    return fcn
end

_create_vec_multiplication(A,u::TU) where {TU} = (x -> _unwrap_vec(A*_wrap_vec(x,u)))
_create_vec_function(A,u::TU) where {TU} = (x -> _unwrap_vec(A(_wrap_vec(x,u))))
_create_vec_backslash(A,u::TU) where {TU} = (x -> _unwrap_vec(A\_wrap_vec(x,u)))

#### WRAPPERS ####
# wrap the vector x in type u, unless u is already a subtype of AbstractVector,
# in which case it just keeps it as is.
_wrap_vec(x::AbstractVector{T},u::TU) where {T,TU} = TU(reshape(x,size(u)...))

# if the vector x is simply a reshaped form of type u, then just get the
# parent of x
_wrap_vec(x::Base.ReshapedArray,u::TU) where {TU} = parent(x)

#### UNWRAPPERS ####
# Usually vec suffices
_unwrap_vec(x) = vec(x)
_unwrap_vec(x::Tuple) = x


text:
# ---
# title: 1209. Remove All Adjacent Duplicates in String II
# id: problem1209
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/remove-all-adjacent-duplicates-in-string-ii/description/>
# hidden: true
# ---
#
# Given a string `s`, a _k_   _duplicate removal_  consists of choosing `k`
# adjacent and equal letters from `s` and removing them causing the left and the
# right side of the deleted substring to concatenate together.
#
# We repeatedly make `k` duplicate removals on `s` until we no longer can.
#
# Return the final string after all such duplicate removals have been made.
#
# It is guaranteed that the answer is unique.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abcd", k = 2
#     Output: "abcd"
#     Explanation: There's nothing to delete.
#
# **Example 2:**
#
#
#
#     Input: s = "deeedbbcccbdaa", k = 3
#     Output: "aa"
#     Explanation: First delete "eee" and "ccc", get "ddbbbdaa"
#     Then delete "bbb", get "dddaa"
#     Finally delete "ddd", get "aa"
#
# **Example 3:**
#
#
#
#     Input: s = "pbbcggttciiippooaais", k = 2
#     Output: "ps"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `2 <= k <= 10^4`
#   * `s` only contains lower case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using NeuroGraphs
using Documenter

DocMeta.setdocmeta!(NeuroGraphs, :DocTestSetup, :(using NeuroGraphs); recursive=true)

makedocs(;
    modules=[NeuroGraphs],
    authors="Zachary P. Christensen <zchristensen7@gmail.com> and contributors",
    repo="https://github.com/JuliaNeuroscience/NeuroGraphs.jl/blob/{commit}{path}#{line}",
    sitename="NeuroGraphs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaNeuroscience.github.io/NeuroGraphs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaNeuroscience/NeuroGraphs.jl",
)


text:
module Variance

require("Corpora.jl")
require("CommonTools.jl")
require("CVMAbstract.jl")
require("DataTypes.jl")

## Imports
using DataTypes
using Base.LinAlg.BLAS.axpy!
using Corpora.Token
using Corpora.text
using CommonTools: tanh, square!, l2sq, l2sq_distance, mult_by_0!
using CommonTools: add!, subtract_vec_matrix!
using Base.LinAlg.BLAS.scal!
using CommonTools.CVMIntermediates

using CVMAbstractModel
import CVMAbstractModel: compose!, store_gradients!

## Exports
export VarianceModel
export compose
export compose!
export store_gradients!
export VarianceParameters

immutable VarianceModel <: AbstractModel
    wordvecdim::Int64;
    margin::Float;
    regularizer::Float;
    lr_init::Float;
    lr_min::Float;
    niter::Int64;
    gamma::Float;
    nnoise::Int64;
    batchsize::Int64;
end

type VarianceParameters <: AbstractParameters
    wordvecs::Matrix{Float};
    gradients::Matrix{Float};
end

VarianceParameters() = VarianceParameters(Array(Float,0,0), Array(Float,0,0))

function l2sq!(v1::DenseArray{Float,1}, v2::DenseArray{Float,1},
                    o::DenseArray{Float,1})
    @simd for i=1:size(v1,1)
        o[i] += (v1[i]-v2[i])^2
    end
end

function compose!(model::VarianceModel, local1::Vector{Float},
                    s::Vector{Token}, o::Vector{Float})
    # Calculate the mean
    m = local1
    mult_by_0!(m)
    #= m = zeros(Float, model.wordvecdim) =#
    for i=1:size(s,1)
        #= Base.LinAlg.BLAS.axpy!(const_one, s[i].vector, m) =#
        add!(s[i].vector, m)
        #@simd for j=1:model.wordvecdim
        #    m[j] += s[i].vector[j]
        #end
    end
    scal!(model.wordvecdim, 1.0/size(s,1), m, 1)
    mult_by_0!(o);
    tmp = 0.0
    for i=1:size(s,1)
        l2sq!(s[i].vector, m, o)
        #= @simd for j=1:model.wordvecdim =#
        #=     tmp = s[i].vector[j] - m[j] =#
        #=     o[j] += tmp * tmp =#
        #=     #o[j] += (s[i].vector[j] - m[j])^2 =#
        #= end =#
    end
    scal!(model.wordvecdim, 1.0/size(s,1), o, 1)

    return o
end

function compose!(model::VarianceModel, s::Vector{Token},
                                intermediate::CVMIntermediates)
    return compose!(model, intermediate.local1, s, intermediate.result)
end

function compose(model::VarianceModel, s::Vector{Token})
    # Calculate the mean
    m = zeros(Float, model.wordvecdim)
    for i=1:size(s,1)
        add!(s[i].vector,m)
    end

    # Calculate the output
    o = zeros(Float, model.wordvecdim)
    for i=1:size(s,1)
        @simd for j=1:model.wordvecdim
            o[j] += (s[i].vector[j] .- m[j])^2
        end
    end
    return o
end

function multstore!(scalar::Float, vec::DenseArray{Float,1},
                    idx, m::DenseArray{Float,2})
    @simd for j=1:size(vec,1)
        m[j,idx] = scalar
        m[j,idx] *= vec[j]
    end
end

function store_gradients!(model::VarianceModel,
                    pre_grad::Vector{Float},
                    s::Vector{Token}, intermediate::CVMIntermediates;
                    intermediates_precomputed::Bool=false)

    #= result = Array(Float64,(model.wordvecdim,size(s,1))) =#
    result = intermediate.ofull
    tmp_scaling = 2.0/size(s,1)
    for i=1:size(s,1)
        multstore!(tmp_scaling, s[i].vector, i, result)
    end
    m = zeros(model.wordvecdim)
    for i=1:size(s,1)
        #= Base.LinAlg.BLAS.axpy!(const_one, s[i].vector, m) =#
        add!(s[i].vector,m)
        #@simd for j=1:model.wordvecdim
        #    m[j] += s[i].vector[j]
        #end
    end
    factor = 2.0/(size(s,1)^2)
    scal!(model.wordvecdim, factor, m,1)
    for i=1:size(s,1)
        subtract_vec_matrix!(m, i, result)
        #= @simd for j=1:model.wordvecdim =#
        #=     result[j,i] -= m[j] =#
        #= end =#
    end

    tmp = intermediate.local1
    #= tmp = Array(Float,model.wordvecdim) =#
    for i=1:size(s,1)
        @simd for j=1:model.wordvecdim
            tmp[j] = result[j,i] * pre_grad[j]
        end
        add!(tmp, s[i].gradient)
    end
end

end

max_stars_repo_name: ogh/binclusion

max_issues_repo_name: ogh/binclusion

max_forks_repo_name: ogh/binclusion


text:
using ..ChemistryFeaturization.AbstractType:
    AbstractFeatureDescriptor, AbstractFeaturization
using ..ChemistryFeaturization.FeatureDescriptor
import ..ChemistryFeaturization.FeatureDescriptor: output_shape
using ..ChemistryFeaturization.Atoms: AtomGraph
using ..ChemistryFeaturization.Utils.ElementFeatureUtils

export GraphNodeFeaturization, encode
export encodable_elements, decode, chunk_vec, output_shape

using DataFrames

"""
    GraphNodeFeaturization(encoded_atom_features::Vector{AtomFeatureDescriptor})
    GraphNodeFeaturization(feature_names, lookup_table, nbins, logspaced, categorical)

A featurization for AtomGraph objects that encodes features associated with each node. Contains a collection of `AbstractAtomFeatureDescriptor` objects, and can be initialized by passing those, or by passing parameters for constructing them.

## Required Arguments
- `feature_names::Vector{String}`: Names of each feature
- `lookup_table::DataFrame` (optional): source of data for features, if not included in built-in `atom_data_df`

## Keyword Arguments
- `nbins::Union{Vector{Integer},Integer}`: Number of bins to use for one-hot encoding of continuous-valued features. Will be ignored for categorical features.
- `logspaced::Union{Vector{Bool},Bool}`: Whether to logarithmically space the bins
- `categorical::Union{Vector{Bool},Bool}`: Whether each feature is categorical or continous-valued.
"""
struct GraphNodeFeaturization <: AbstractFeaturization
    features::Vector{<:AbstractAtomFeatureDescriptor} # TODO: this might need another name
end

# NB: this constructor only works if every feature is an ElementFeature
function GraphNodeFeaturization(
    feature_names::Vector{String},
    lookup_table::Union{DataFrame,Nothing} = nothing;
    nbins::Union{Vector{<:Integer},Integer,Nothing} = nothing,
    logspaced::Union{Vector{Bool},Bool,Nothing} = nothing,
    categorical::Union{Vector{Bool},Bool,Nothing} = nothing,
)
    num_features = length(feature_names)
    local lookup_table_here, logspaced_here, categorical_here, nbins_here
    if isnothing(lookup_table)
        lookup_table_here = atom_data_df
    else
        # need to merge them in case some data is in one place and some the other
        lookup_table_here = outerjoin(atom_data_df, lookup_table, on = :Symbol)
    end

    if isnothing(logspaced)
        logspaced_here = map(fn -> default_log(fn, lookup_table_here), feature_names)
    else
        logspaced_here = get_param_vec(logspaced, num_features)
    end

    if isnothing(categorical)
        categorical_here =
            map(fn -> default_categorical(fn, lookup_table_here), feature_names)
    else
        categorical_here = get_param_vec(categorical, num_features)
    end

    if isnothing(nbins)
        nbins_here = [default_nbins for i = 1:num_features]
    else
        nbins_here = get_param_vec(nbins, num_features, pad_val = default_nbins)
    end

    afs = map(zip(feature_names, nbins_here, logspaced_here, categorical_here)) do args
        ElementFeatureDescriptor(
            args[1],
            lookup_table_here,
            nbins = args[2],
            logspaced = args[3],
            categorical = args[4],
        )
    end
    GraphNodeFeaturization(afs)
end

# pretty printing, short
function Base.show(io::IO, fzn::GraphNodeFeaturization)
    st = "GraphNodeFeaturization encoding $(length(fzn.features)) features"
    print(io, st)
end

# pretty printing, long
function Base.show(io::IO, ::MIME"text/plain", fzn::GraphNodeFeaturization)
    st = "GraphNodeFeaturization encoding $(length(fzn.features)) features:\n"
    for feature in fzn.features
        st = string(st, "\t", feature, "\n")
    end
    st = st[1:end-1]
    print(io, st)
end

output_shape(fzn::GraphNodeFeaturization) = sum(output_shape.(fzn.features))

encodable_elements(fzn::GraphNodeFeaturization) =
    intersect([encodable_elements(f) for f in fzn.features]...)

"""
    chunk_vec(vec, nbins)

Helper function that divides up an already-constructed feature vector into "chunks" (one for each feature) of lengths specified by the vector nbins.

Sum of nbins should be equal to the length of vec.

# Examples
```jldoctest
julia> chunk_vec([1,0,0,1,0], [3,2])
2-element Array{Array{Bool,1},1}:
 [1, 0, 0]
 [1, 0]
 ```
"""
function chunk_vec(vec::Vector{<:Real}, nbins::Vector{<:Integer})
    chunks = fill(Bool[], size(nbins, 1))
    @assert length(vec) == sum(nbins) "Total number of bins doesn't match length of feature vector."
    for i = 1:size(nbins, 1)
        if i == 1
            start_ind = 1
        else
            start_ind = sum(nbins[1:i-1]) + 1
        end
        chunks[i] = vec[start_ind:start_ind+nbins[i]-1]
    end
    return chunks
end

function encode(fzn::GraphNodeFeaturization, ag::AtomGraph)
    encoded = reduce(vcat, map((x) -> encode(x, ag), fzn.features))
    return encoded
end

function decode(fzn::GraphNodeFeaturization, encoded::Matrix{<:Real})
    num_atoms = size(encoded, 2)
    nbins = [output_shape(f) for f in fzn.features]
    local decoded = Dict{Integer,Dict{String,Any}}()
    for i = 1:num_atoms
        #println("atom $i")
        chunks = chunk_vec(encoded[:, i], nbins)
        decoded[i] = Dict{String,Any}()
        for (chunk, f) in zip(chunks, fzn.features)
            #println("    $(f.name): $(decode(f, chunk))")
            decoded[i][f.name] = decode(f, chunk)
        end
    end
    return decoded
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Rotations from an Earth-Centered, Earth-Fixed (ECEF) reference frame to
#   another ECEF reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export r_ecef_to_ecef

"""
    r_ecef_to_ecef([T,] ECEFo, ECEFf, JD_UTC::Number, eop_data)

Compute the rotation from an Earth-Centered, Earth-Fixed (`ECEF`) reference
frame to another ECEF reference frame at the Julian Day [UTC] `JD_UTC`. The
rotation description that will be used is given by `T`, which can be `DCM` or
`Quaternion`. The origin ECEF frame is selected by the input `ECEFo` and the
destination ECEF frame is selected by the input `ECEFf`. The model used to
compute the rotation is specified by the selection of the origin and destination
frames. Currently, there are two models supported: IAU-76/FK5 and IAU-2006 with
2010 conventions.

# Rotation description

The rotations that aligns the origin ECEF frame with the destination ECEF frame
can be described by Direction Cosine Matrices or Quaternions. This is selected
by the parameter `T`.

The possible values are:

- `DCM`: The rotation will be described by a Direction Cosine Matrix.
- `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported.**

# ECEF Frame

The supported ECEF frames for both origin `ECEFo` and destination `ECEFf` are:

- `ITRF()`: ECEF will be selected as the International Terrestrial Reference
    Frame (ITRF).
- `PEF()`: ECEF will be selected as the Pseudo-Earth Fixed (PEF) reference
    frame.
- `TIRS()`: ECEF will be selected as the Terrestrial Intermediate Reference
    System (TIRS).

# EOP Data

The conversion between the supported ECEF frames **always** depends on EOP Data
(see [`get_iers_eop`](@ref) and [`read_iers_eop`](@ref)). If IAU-76/FK5 model is
used, then the type of `eop_data` must be [`EOPData_IAU1980`](@ref). Otherwise,
if IAU-2006/2010 model is used, then the type of `eop_data` must be
[`EOPData_IAU2000A`](@ref).

# Returns

The rotation description represented by `T` that rotates the ECEF reference
frame into alignment with the ECI reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop();

julia> r_ecef_to_ecef(PEF(), ITRF(), date_to_jd(1986, 6, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  1.0          0.0         -4.34677e-7
 -6.29476e-13  1.0         -1.44815e-6
  4.34677e-7   1.44815e-6   1.0

julia> r_ecef_to_ecef(Quaternion, PEF(), ITRF(), date_to_jd(1986, 6, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 1.0 - 7.24073e-7⋅i + 2.17339e-7⋅j + 2.17339e-7⋅k

julia> eop_IAU2000A = get_iers_eop(Val(:IAU2000A));

julia> r_ecef_to_ecef(TIRS(), ITRF(), date_to_jd(1986, 6, 19, 21, 35, 0), eop_IAU2000A)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  1.0          3.08408e-11  -4.34677e-7
 -3.14703e-11  1.0          -1.44815e-6
  4.34677e-7   1.44815e-6    1.0

julia> r_ecef_to_ecef(Quaternion, TIRS(), ITRF(), date_to_jd(1986, 6, 19, 21, 35, 0), eop_IAU2000A)
Quaternion{Float64}:
  + 1.0 - 7.24073e-7⋅i + 2.17339e-7⋅j + 2.17339e-7⋅k
```
"""
@inline function r_ecef_to_ecef(
    T_ECEFo::T_ECEFs,
    T_ECEFf::T_ECEFs,
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return r_ecef_to_ecef(DCM, T_ECEFo, T_ECEFf, JD_UTC, eop_data)
end

@inline function r_ecef_to_ecef(
    T_ECEFo::T_ECEFs_IAU_2006,
    T_ECEFf::T_ECEFs_IAU_2006,
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_ecef_to_ecef(DCM, T_ECEFo, T_ECEFf, JD_UTC, eop_data)
end

################################################################################
#                                  IAU-76/FK5
################################################################################

#                                 ITRF <=> PEF
# ==============================================================================

function r_ecef_to_ecef(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:PEF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Return the rotation.
    return r_itrf_to_pef_fk5(T, x_p, y_p)
end

@inline function r_ecef_to_ecef(
    T::T_ROT,
    T_ECEFo::Val{:PEF},
    T_ECEFf::Val{:ITRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_ecef_to_ecef(T, T_ECEFf, T_ECEFo, JD_UTC, eop_data))
end

################################################################################
#                                IAU-2006/2010
################################################################################

#                                ITRF <=> TIRS
# ==============================================================================

function r_ecef_to_ecef(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:TIRS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Return the rotation.
    return r_itrf_to_tirs_iau2006(T, JD_UTC, x_p, y_p)
end

@inline function r_ecef_to_ecef(
    T::T_ROT,
    T_ECEFo::Val{:TIRS},
    T_ECEFf::Val{:ITRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return inv_rotation(r_ecef_to_ecef(T, T_ECEFf, T_ECEFo, JD_UTC, eop_data))
end


text:
function payoff(S::AbstractMatrix{num}, euPayoff::EuropeanPayoff, rfCurve::abstractZeroRateCurve, T1::num2 = maturity(euPayoff)) where {abstractZeroRateCurve <: AbstractZeroRateCurve, num <: Number, num2 <: Number}
    r = rfCurve.r
    T = euPayoff.T
    NStep = size(S, 2) - 1
    index1 = round(UInt, T / T1 * NStep) + 1
    @views ST = S[:, index1]
    phi(Sti::numtype_) where {numtype_ <: Number} = payout(Sti, euPayoff)
    payoff2 = phi.(ST)

    return payoff2 * exp(-integral(r, T))
end