text:
friction_coefficients(model::Rocket{T}) where T = T[]

function RoboDojo.indices_z(model::Rocket)
    nq = model.nq
    q = collect(1:nq)
    γ = collect(1:0)
    sγ = collect(1:0)
    ψ = collect(1:0)
    b = collect(1:0)
    sψ = collect(1:0)
    sb = collect(1:0)
    IndicesZ(q, γ, sγ, ψ, b, sψ, sb)
end

function RoboDojo.indices_θ(model::Rocket; nf=0)
    nq = model.nq
    nu = model.nu
    nw = model.nw

    q1 = collect(1:nq)
    q2 = collect(1:0)
    u = collect(nq .+ (1:nu))
    w = collect(1:0)
    f = collect(1:0)
    h = collect(nq + nu .+ (1:1))

    Indicesθ(q1, q2, u, w, f, h)
end

RoboDojo.num_var(model::Rocket) = model.nq

RoboDojo.nominal_configuration(model::Rocket) = zeros(model.nq)

function RoboDojo.indices_optimization(model::Rocket)
    nq = model.nq
    nz = nq
    IndicesOptimization(
        nz,
        nz,
        [collect(1:0), collect(1:0)],
        [collect(1:0), collect(1:0)],
        Vector{Vector{Vector{Int}}}(),
        Vector{Vector{Vector{Int}}}(),
        collect(1:(nq)),
        collect(1:0),
        collect(1:0),
        Vector{Vector{Int}}(),
        collect(1:0))
end

function RoboDojo.initialize_z!(z, model::Rocket, idx::IndicesZ, q)
    z[idx.q] .= q
end


text:
using Knet, KnetModules
import KnetModules.forward

# -- ImageNet --
const basic_channels = [64, 128, 256, 512]
const bneck_channels = 4basic_channels

resnet18(;o...) =
    ResNet(BasicBlock, basic_channels, [2, 2, 2, 2]; o...)

resnet34(;o...) =
    ResNet(BasicBlock, basic_channels, [3, 4, 6, 3]; o...)

resnet50(;o...) =
    ResNet(Bottleneck, bneck_channels, [3, 4, 6, 3]; o...)

resnet101(;o...) =
    ResNet(Bottleneck, bneck_channels, [3, 4, 23, 3]; o...)

resnet152(;o...) =
    ResNet(Bottleneck, bneck_channels, [3, 8, 36, 3]; o...)

# -- CIFAR --
for (m, d) in zip([:resnet20, :resnet32, :resnet56, :resnet110],
                  [20, 32, 56, 110])
    eval(:(
        $(m)(;o...) = ResNetCifar($(d); o...)
    ))
end

type BasicBlock <: KnetModule
    conv1; bn1
    conv2; bn2
    downsample
end

function BasicBlock(input::Int, output::Int;
                    stride=1, downsample=nothing)
    conv1 = Conv(3,3,input, output;
                  padding=1, stride=stride,
                  bias=false)
    bn1 = BatchNorm(output)
    conv2 = Conv(3, 3, output, output;
                  padding=1, stride=1,
                  bias=false)
    bn2 = BatchNorm(output)
    return BasicBlock(conv1, bn1, conv2, bn2, downsample)
end

function forward(ctx, b::BasicBlock, x)
    o = relu.(@mc b.bn1(@mc b.conv1(x)))
    o = @mc b.bn2(@mc b.conv2(o))
    if b.downsample !== nothing
        x = @mc b.downsample(x)
    end
    return relu.(o .+ x)
end

type Bottleneck <: KnetModule
    conv1; bn1
    conv2; bn2
    conv3; bn3
    downsample
end

let s = true
    global stride3!, stride1!, stride3
    stride3!() = (s=true)
    stride1!() = (s=false)
    stride3() = s
end

function Bottleneck(input::Int, output::Int;
                    stride=1,downsample=nothing)
    # middle dim.
    planes = div(output, 4)
    # layers
    conv1 = Conv(1, 1, input, planes;
                 bias=false,
                 stride=(stride3() ? 1 : stride))
    bn1 = BatchNorm(planes)
    conv2 = Conv(3, 3, planes, planes;
                 bias=false,
                 stride=(stride3() ? stride: 1),
                 padding=1)
    bn2 = BatchNorm(planes)
    conv3 = Conv(1, 1, planes, output;
                  bias=false)
    bn3 = BatchNorm(output)
    # return
    return Bottleneck(conv1, bn1, conv2, bn2,
                      conv3, bn3, downsample)
end

function forward(ctx, b::Bottleneck, x)
    o = relu.(@mc b.bn1(@mc b.conv1(x)))
    o = relu.(@mc b.bn2(@mc b.conv2(o)))
    o = @mc b.bn3(@mc b.conv3(o))
    if b.downsample !== nothing
        x = @mc b.downsample(x)
    end
    return relu.(o .+ x)
end

function _make_layer(block::Type, input::Int, output::Int,
                     nblocks::Int; stride=1)
    downsample = nothing
    if stride !== 1 || input !== output
        downsample = Sequential(
            Conv(1, 1, input, output;
                stride=stride, bias=false),
            BatchNorm(output))
    end
    modules = Sequential()
    add!(modules, block(input, output;
                        stride=stride,
                        downsample=downsample))
    for i = 2:nblocks
        add!(modules, block(output, output))
    end
    return modules
end

type ResNet <: KnetModule
    conv1::Conv
    bn1::BatchNorm
    layers::Array{Sequential, 1}
    output::Union{Void, Linear}
    dataset::Symbol
    stage::Int
    avg_pool::Bool
end

function ResNet(block::Type,
                nchannels::Array{Int, 1},
                nrepeats::Array{Int, 1};
                nclasses=1000,
                stage=0)
    conv1 = Conv(7, 7, 3, 64; stride=2, padding=3, bias=false)
    bn1 = BatchNorm(64)
    layers = []
    nchannels = [64, nchannels...]
    for i = 1:4
        push!(layers, _make_layer(
            block, nchannels[i],
            nchannels[i+1], nrepeats[i];
            stride = 1+Int(i>1)))
        (i == stage || (stage==5 && i==4)) &&
            return ResNet(conv1, bn1, layers, nothing, stage==5)
    end

    output = Linear(nclasses, nchannels[end])
    return ResNet(conv1, bn1, layers, output, :imagenet, stage, true)
end

function ResNetCifar(depth::Int; nclasses=10)
    conv1 = Conv(3, 3, 3, 16; bias=false, padding=1)
    bn1 = BatchNorm(16)
    n = Int((depth-2)/3)
    rep = div(n, 2)
    layers = [
        _make_layer(BasicBlock, 16, 16, rep),
        _make_layer(BasicBlock, 16, 32, rep; stride=2),
        _make_layer(BasicBlock, 32, 64, rep; stride=2),
    ]
    output = Linear(nclasses, 64)
    return ResNet(conv1, bn1, layers, output, :cifar, true)
end

function forward(ctx, r::ResNet, x)
    # initial conv
    if r.dataset == :cifar
        o = relu.(@mc r.bn1(@mc r.conv1(x)))
    else
        o = pool(relu.(@mc r.bn1(@mc r.conv1(x)));
                 window=3, stride=2, padding=Int(stride3()))
    end
    # residual blocks
    for (i,l) in enumerate(r.layers)
        o = @mc l(o)
        r.stage == i && return o
    end
    # output layer
    if r.avg_pool
        o = pool(o; mode=2, window=size(o,1,2))
        r.stage == 5 && return o
    end
    if r.output !== nothing
        o = @mc r.output(o)
    end
    return o
end

ext: jl


text:
export  Rectifier,
        square, dsquare, d2square,
        softplus, dsoftplus, d2softplus, invsoftplus,
        explinearunit, dexplinearunit, d2explinearunit, invexplinearunit,
        inverse!, inverse, vinverse,
        grad_x!, grad_x, vgrad_x,
        grad_x_logeval, grad_x_logeval!, vgrad_x_logeval,
        hess_x_logeval, hess_x_logeval!, vhess_x_logeval,
        hess_x!, hess_x, vhess_x,
        evaluate!, vevaluate

"""
$(TYPEDEF)

This structure defines a continuous rectifier g,
i.e. a positive and monotonically increasing function
(e.g. square function, exponential , softplus, explinearunit).

## Fields

$(TYPEDFIELDS)

"""
struct Rectifier
    T::String
end

square(x) = x^2
dsquare(x) = 2.0*x
d2square(x) = 2.0

# Softplus tools
softplus(x) = (log(1.0 + exp(-abs(log(2.0)*x))) + max(log(2.0)*x, 0.0))/log(2.0)
dsoftplus(x) = 1/(1 + exp(-log(2.0)*x))
d2softplus(x) = log(2.0)/(2.0*(1.0 + cosh(log(2.0)*x)))
invsoftplus(x) = min(log(exp(log(2.0)*x) - 1.0)/log(2.0), x)

explinearunit(x) = x < 0.0 ? exp(x) : x + 1.0
dexplinearunit(x) = x < 0.0 ? exp(x) : 1.0
d2explinearunit(x) = x < 0.0 ? exp(x) : 0.0
invexplinearunit(x) = x < 1.0 ? log(x) : x - 1.0

# Type of the rectifier should be in the following list:
# "squared", "exponential", "softplus", "explinearunit"

Rectifier() = Rectifier("softplus")

function (g::Rectifier)(x)
    if g.T=="squared"
        return square(x)
    elseif g.T=="exponential"
        return exp(x)
    elseif g.T=="softplus"
        return softplus(x)
    elseif g.T=="explinearunit"
        return explinearunit(x)
    end
end

function evaluate!(result, g::Rectifier, x)
    @assert size(result,1) == size(x,1) "Dimension of result and x don't match"
    if g.T=="squared"
        vmap!(square, result, x)
        return result
    elseif g.T=="exponential"
        vmap!(exp, result, x)
        return result
    elseif g.T=="softplus"
        vmap!(softplus, result, x)
        return result
    elseif g.T=="explinearunit"
        vmap!(explinearunit, result, x)
        return result
    end
end

vevaluate(g::Rectifier, x) = evaluate!(zero(x), g, x)

function inverse(g::Rectifier, x)
    @assert x>=0 "Input to rectifier is negative"
    if g.T=="squared"
        error("squared rectifier is not invertible")
    elseif g.T=="exponential"
        return log(x)
    elseif g.T=="softplus"
        return invsoftplus(x)
    elseif g.T=="explinearunit"
        return invexplinearunit(x)
    end
end

function inverse!(result, g::Rectifier, x)
    @assert all(x .> 0) "Input to rectifier is negative"
    @assert size(result,1) == size(x,1) "Dimension of result and x don't match"
    if g.T=="squared"
        error("squared rectifier is not invertible")
    elseif g.T=="exponential"
        vmap!(log, result, x)
        return result
    elseif g.T=="softplus"
        vmap!(invsoftplus, result, x)
        return result
    elseif g.T=="explinearunit"
        vmap!(invexplinearunit, result, x)
        return result
    end
end

vinverse(g::Rectifier, x)  = inverse!(zero(x), g, x)

function grad_x(g::Rectifier, x)
    if g.T=="squared"
        return dsquare(x)
    elseif g.T=="exponential"
        return exp(x)
    elseif g.T=="softplus"
        return dsoftplus(x)
    elseif g.T=="explinearunit"
        return dexplinearunit(x)
    end
end

function grad_x!(result, g::Rectifier, x)
    @assert size(result,1) == size(x,1) "Dimension of result and x don't match"
    if g.T=="squared"
        vmap!(dsquare, result, x)
        return result
    elseif g.T=="exponential"
        vmap!(exp, result, x)
        return result
    elseif g.T=="softplus"
        vmap!(dsoftplus, result, x)
        return result
    elseif g.T=="explinearunit"
        vmap!(dexplinearunit, result, x)
        return result
    end
end

vgrad_x(g::Rectifier, x) = grad_x!(zero(x), g, x)

""" Compute g′(x)/g(x) i.e d/dx log(g(x))"""
function grad_x_logeval(g::Rectifier, x::T) where {T <: Real}
    if g.T=="squared"
        return dsquare(x)/square(x)
    elseif g.T=="exponential"
        return 1.0
    elseif g.T=="softplus"
        return dsoftplus(x)/softplus(x)
    elseif g.T=="explinearunit"
        return dexplinearunit(x)/explinearunit(x)
    end
end

function grad_x_logeval!(result, g::Rectifier, x)
    @assert size(result,1) == size(x,1) "Dimension of result and x don't match"
    if g.T=="squared"
        vmap!(xi->dsquare(xi)/square(xi), result, x)
        return result
    elseif g.T=="exponential"
        vmap!(1.0, result, x)
        return result
    elseif g.T=="softplus"
        vmap!(xi->dsoftplus(xi)/softplus(xi), result, x)
        return result
    elseif g.T=="explinearunit"
        vmap!(xi->dexplinearunit(xi)/explinearunit(xi), result, x)
        return result
    end
end

vgrad_x_logeval(g::Rectifier, x) = grad_x_logeval!(zero(x), g, x)

# Compute (g″(x)g(x)-g′(x)^2)/g(x) i.e d^2/dx^2 log(g(x))
function hess_x_logeval(g::Rectifier, x::T) where {T <: Real}
    if g.T=="squared"
        return (d2square(x)*square(x) - dsquare(x)^2)/square(x)^2
    elseif g.T=="exponential"
        return 0.0
    elseif g.T=="softplus"
        return (d2softplus(x)*softplus(x) - dsoftplus(x)^2)/softplus(x)^2
    elseif g.T=="explinearunit"
        return (d2explinearunit(x)*explinearunit(x) - dexplinearunit(x)^2)/explinearunit(x)^2
    end
end

function hess_x_logeval!(result, g::Rectifier, x)
    @assert size(result,1) == size(x,1) "Dimension of result and x don't match"
    if g.T=="squared"
        vmap!(xi->(d2square(xi)*square(xi) - dsquare(xi)^2)/square(xi)^2, result, x)
        return result
    elseif g.T=="exponential"
        vmap!(0.0, result, x)
        return result
    elseif g.T=="softplus"
        vmap!(xi->(d2softplus(xi)*softplus(xi) - dsoftplus(xi)^2)/softplus(xi)^2, result, x)
        return result
    elseif g.T=="explinearunit"
        vmap!(xi->(d2explinearunit(xi)*explinearunit(xi) - dexplinearunit(xi)^2)/explinearunit(xi)^2, result, x)
        return result
    end
end

vhess_x_logeval(g::Rectifier, x) = hess_x_logeval!(zero(x), g, x)

function hess_x(g::Rectifier, x::T) where {T <: Real}
    if g.T=="squared"
        return d2square(x)
    elseif g.T=="exponential"
        return exp(x)
    elseif g.T=="softplus"
        return d2softplus(x)
    elseif g.T=="explinearunit"
        return d2explinearunit(x)
    end
end

function hess_x!(result, g::Rectifier, x)
    @assert size(result,1) == size(x,1) "Dimension of result and x don't match"
    if g.T=="squared"
        vmap!(d2square, result, x)
        return result
    elseif g.T=="exponential"
        vmap!(exp, result, x)
        return result
    elseif g.T=="softplus"
        vmap!(d2softplus, result, x)
        return result
    elseif g.T=="explinearunit"
        vmap!(d2explinearunit, result, x)
        return result
    end
end

vhess_x(g::Rectifier, x) = hess_x!(zero(x), g, x)


text:
################################################################################
#                                                                              #
#                                 OptiMode.jl:                                 #
#                Differentiable Electromagnetic Eigenmode Solver               #
#                                 (an attempt)                                 #
#                                                                              #
################################################################################

module OptiMode

## Imports ##
using LinearAlgebra
using LinearMaps
using ArrayInterface
using StaticArrays
using HybridArrays
# using StructArrays
# using RecursiveArrayTools
using FFTW
using AbstractFFTs
using GeometryPrimitives
using Roots
using ChainRules
using Zygote
using ForwardDiff
using LoopVectorization
using Tullio
using IterativeSolvers
using Rotations
using Colors: Color, RGB, RGBA, @colorant_str
using RuntimeGeneratedFunctions
# using ModelingToolkit, Unitful #, Latexify
using Symbolics
using SymbolicUtils
using UnicodePlots
using Makie
# using KrylovKit
using DFTK: LOBPCG
using KrylovKit: eigsolve
using StaticArrays: Dynamic, SVector
using ChainRulesCore
using ChainRulesCore: rrule, frule, @thunk, @non_differentiable, @not_implemented, NoTangent, ZeroTangent, AbstractZero
# using ChainRules: ChainRulesCore, @non_differentiable,  NoTangent, @thunk, @not_implemented, AbstractZero
using Zygote: Buffer, bufferfrom, @ignore, @adjoint, ignore, dropgrad, forwarddiff, Numeric, literal_getproperty, accum
using GeometryPrimitives: Cylinder
using Makie.GeometryBasics

RuntimeGeneratedFunctions.init(@__MODULE__)
# using Plots: plot, heatmap, plot!, heatmap!
# using MaxwellFDM: kottke_avg_param
## Exports ##
# export plot_ε, test_shapes, ridge_wg, circ_wg, trap_wg, trap_wg2, plot, heatmap, SHM3
export SHM3
export plot_data, uplot, uplot!, xlims, ylims

# Import methods that we will overload for custom types
import Base: size, eltype
import LinearAlgebra: mul!
import ChainRulesCore: rrule

# FFTW settings
FFTW.set_num_threads(1)     # chosen for thread safety when combined with other parallel code, consider increasing

## Abbreviations, aliases, etc. ##
SHM3 = SHermitianCompact{3,Float64,6}   # static Hermitian 3×3 matrix Type, renamed for brevity

## generic plotting methods/aliases
import Colors: JULIA_LOGO_COLORS
logocolors = JULIA_LOGO_COLORS
uplot(x::AbstractVector, y::AbstractVector; kwargs...) = UnicodePlots.lineplot(x, y ; kwargs...)
uplot!(plt::UnicodePlots.Plot, x::AbstractVector, y::AbstractVector; kwargs...) = UnicodePlots.lineplot!(plt, x, y ; kwargs...)
plot_data(x) = x    # nontrivial methods defined for custom types

xlims(x::Real) = x
xlims(x::Complex) = real(x)
xlims(v::AbstractVector{<:Number}) = [extrema(real(v))...]
xlims(vs::AbstractVector) = [extrema( vcat( xlims.(vs) ) )...]
function xlims(a,b)
	xlims_a, xlims_b = xlims(a), xlims(b)
	[ min(xlims_a[1], xlims_b[1]), max(xlims_a[2], xlims_b[2]) ]
end
xlims(plt::UnicodePlots.Plot) = parse.(Float64,getfield.(unique(plt.decorations),:second))

ylims(x::Real) = x
ylims(x::Complex) = real(x)
ylims(v::AbstractVector{<:Number}) = [extrema(real(v))...]
ylims(vs::AbstractVector) = [extrema( vcat( ylims.(vs) ) )...]
function ylims(fns::AbstractVector{<:Function};xlims=(0.5,1.8),nvals=100)
	xrng = range(xlims[1],xlims[2];length=nvals)
	[ extrema(hcat(map.(fns,(xrng,) )...) )... ]
end
function ylims(a,b)
	ylims_a, ylims_b = ylims(a), ylims(b)
	[ min(ylims_a[1], ylims_b[1]), max(ylims_a[2], ylims_b[2]) ]
end
ylims(plt::UnicodePlots.Plot) = parse.(Float64,getfield.(unique(plt.labels_left),:second))

# function ylims(plt::UnicodePlots.Plot, new ; kwargs...)
# 	ylims_plt = parse.(Float64,getfield.(unique(plt.labels_left),:second))
# 	ylims(ylims_plt,ylims(new))
# end

## Add methods to external packages ##
LinearAlgebra.ldiv!(c,A::LinearMaps.LinearMap,b) = mul!(c,A',b)

## Includes ##

# include("plot.jl")
include("linalg.jl")
include("epsilon.jl")
include("materials.jl")
include("grid.jl")
include("geometry.jl")
include("maxwell.jl")
include("constraints.jl")
include("fields.jl")
include("smooth.jl")
include("solve.jl")
include("grads.jl")
include("explicit.jl")
# include("optimize.jl")

## Definitions ##

## More includes ##

end


text:
include("tools_tests.jl")
include("fit_tests.jl")

println("All tests complete")


text:
import Base: size

"""
    nonbacktracking_matrix(g)
Return a non-backtracking matrix `B` and an edgemap storing the oriented
edges' positions in `B`.
Given two arcs ``A_{i j}` and `A_{k l}` in `g`, the
non-backtraking matrix ``B`` is defined as
``B_{A_{i j}, A_{k l}} = δ_{j k} * (1 - δ_{i l})``
"""
function nonbacktracking_matrix(g::AGraph)
    E = edgetype(g)
    edgeidmap = Dict{E,Int}()
    m = 0
    for e in edges(g)
        m += 1
        edgeidmap[e] = m
    end

    if !is_directed(g)
        for e in edges(g)
            m += 1
            edgeidmap[reverse(e)] = m
        end
    end

    B = zeros(Float64, m, m)

    for (e, u) in edgeidmap
        i, j = src(e), dst(e)
        for e in in_edges(g, i)
            k = src(e)
            k == j && continue
            v = edgeidmap[e]
            B[v, u] = 1
        end
    end

    return B, edgeidmap
end


text:
using Pkg
using Dates

PROJECT_DIR = abspath(@__DIR__) |> dirname

# Check whether triggering comment contains '[benchmark]' command
#  and return immediately if not.
# Reference implementation of listening to GH comments:
#  https://github.com/repetitive/actions/blob/master/clojure/entrypoint.sh

# Prepare packages
Pkg.build("Turing")
Pkg.add("BenchmarkTools")

# Run `MCMCBenchmarks`
# Add code

# Report results back to PRs
# Add code


text:
# Invertible conditional HINT network from Kruse et. al (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: February 2020

export NetworkConditionalHINT, NetworkConditionalHINT3D

"""
    CH = NetworkConditionalHINT(n_in, n_hidden, depth; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

    CH = NetworkConditionalHINT3D(n_in, n_hidden, depth; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

 Create a conditional HINT network for data-driven generative modeling based
 on the change of variables formula.

 *Input*:

 - 'n_in': number of input channels

 - `n_hidden`: number of hidden units in residual blocks

 - `depth`: number network layers

 - `k1`, `k2`: kernel size for first and third residual layer (`k1`) and second layer (`k2`)

 - `p1`, `p2`: respective padding sizes for residual block layers

 - `s1`, `s2`: respective strides for residual block layers

 *Output*:

 - `CH`: conditioinal HINT network

 *Usage:*

 - Forward mode: `Zx, Zy, logdet = CH.forward(X, Y)`

 - Inverse mode: `X, Y = CH.inverse(Zx, Zy)`

 - Backward mode: `ΔX, X = CH.backward(ΔZx, ΔZy, Zx, Zy)`

 *Trainable parameters:*

 - None in `CH` itself

 - Trainable parameters in activation normalizations `CH.AN_X[i]` and `CH.AN_Y[i]`,
 and in coupling layers `CH.CL[i]`, where `i` ranges from `1` to `depth`.

 See also: [`ActNorm`](@ref), [`ConditionalLayerHINT!`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct NetworkConditionalHINT <: InvertibleNetwork
    AN_X::AbstractArray{ActNorm, 1}
    AN_Y::AbstractArray{ActNorm, 1}
    CL::AbstractArray{ConditionalLayerHINT, 1}
    logdet::Bool
    is_reversed::Bool
end

@Flux.functor NetworkConditionalHINT

# Constructor
function NetworkConditionalHINT(n_in, n_hidden, depth; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, logdet=true, ndims=2)

    AN_X = Array{ActNorm}(undef, depth)
    AN_Y = Array{ActNorm}(undef, depth)
    CL = Array{ConditionalLayerHINT}(undef, depth)

    # Create layers
    for j=1:depth
        AN_X[j] = ActNorm(n_in; logdet=logdet)
        AN_Y[j] = ActNorm(n_in; logdet=logdet)
        CL[j] = ConditionalLayerHINT(n_in, n_hidden; permute=true, k1=k1, k2=k2, p1=p1, p2=p2,
                                     s1=s1, s2=s2, logdet=logdet, ndims=ndims)
    end

    return NetworkConditionalHINT(AN_X, AN_Y, CL, logdet, false)
end

NetworkConditionalHINT3D(args...;kw...) = NetworkConditionalHINT(args...; kw..., ndims=3)

# Forward pass and compute logdet
function forward(X, Y, CH::NetworkConditionalHINT; logdet=nothing)
    isnothing(logdet) ? logdet = (CH.logdet && ~CH.is_reversed) : logdet = logdet

    depth = length(CH.CL)
    logdet_ = 0f0
    for j=1:depth
        logdet ? (X_, logdet1) = CH.AN_X[j].forward(X) : X_ = CH.AN_X[j].forward(X)
        logdet ? (Y_, logdet2) = CH.AN_Y[j].forward(Y) : Y_ = CH.AN_Y[j].forward(Y)
        logdet ? (X, Y, logdet3) = CH.CL[j].forward(X_, Y_) : (X, Y) = CH.CL[j].forward(X_, Y_)
        logdet && (logdet_ += (logdet1 + logdet2 + logdet3))
    end
    logdet ? (return X, Y, logdet_) : (return X, Y)
end

# Inverse pass and compute gradients
function inverse(Zx, Zy, CH::NetworkConditionalHINT; logdet=nothing)
    isnothing(logdet) ? logdet = (CH.logdet && CH.is_reversed) : logdet = logdet

    depth = length(CH.CL)
    logdet_ = 0f0
    for j=depth:-1:1
        logdet ? (Zx_, Zy_, logdet1) = CH.CL[j].inverse(Zx, Zy; logdet=true) : (Zx_, Zy_) = CH.CL[j].inverse(Zx, Zy; logdet=false)
        logdet ? (Zy, logdet2) = CH.AN_Y[j].inverse(Zy_; logdet=true) : Zy = CH.AN_Y[j].inverse(Zy_; logdet=false)
        logdet ? (Zx, logdet3) = CH.AN_X[j].inverse(Zx_; logdet=true) : Zx = CH.AN_X[j].inverse(Zx_; logdet=false)
        logdet && (logdet_ += (logdet1 + logdet2 + logdet3))
    end
    logdet ? (return Zx, Zy, logdet_) : (return Zx, Zy)
end

# Backward pass and compute gradients
function backward(ΔZx, ΔZy, Zx, Zy, CH::NetworkConditionalHINT; set_grad::Bool=true)
    depth = length(CH.CL)
    if ~set_grad
        Δθ = Array{Parameter, 1}(undef, 0)
        CH.logdet && (∇logdet = Array{Parameter, 1}(undef, 0))
    end
    for j=depth:-1:1
        if set_grad
            ΔZx_, ΔZy_, Zx_, Zy_ = CH.CL[j].backward(ΔZx, ΔZy, Zx, Zy)
            ΔZx, Zx = CH.AN_X[j].backward(ΔZx_, Zx_)
            ΔZy, Zy = CH.AN_Y[j].backward(ΔZy_, Zy_)
        else
            if CH.logdet
                ΔZx_, ΔZy_, Δθcl, Zx_, Zy_, ∇logdetcl = CH.CL[j].backward(ΔZx, ΔZy, Zx, Zy; set_grad=set_grad)
                ΔZx, Δθx, Zx, ∇logdetx = CH.AN_X[j].backward(ΔZx_, Zx_; set_grad=set_grad)
                ΔZy, Δθy, Zy, ∇logdety = CH.AN_Y[j].backward(ΔZy_, Zy_; set_grad=set_grad)
                ∇logdet = cat(∇logdetx, ∇logdety, ∇logdetcl, ∇logdet; dims=1)
            else
                ΔZx_, ΔZy_, Δθcl, Zx_, Zy_ = CH.CL[j].backward(ΔZx, ΔZy, Zx, Zy; set_grad=set_grad)
                ΔZx, Δθx, Zx = CH.AN_X[j].backward(ΔZx_, Zx_; set_grad=set_grad)
                ΔZy, Δθy, Zy = CH.AN_Y[j].backward(ΔZy_, Zy_; set_grad=set_grad)
            end
            Δθ = cat(Δθx, Δθy, Δθcl, Δθ; dims=1)
        end
    end
    if set_grad
        return ΔZx, ΔZy, Zx, Zy
    else
        CH.logdet ? (return ΔZx, ΔZy, Δθ, Zx, Zy, ∇logdet) : (return ΔZx, ΔZy, Δθ, Zx, Zy)
    end
end

# Backward reverse pass and compute gradients
function backward_inv(ΔX, ΔY, X, Y, CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=1:depth
        ΔX_, X_ = backward_inv(ΔX, X, CH.AN_X[j])
        ΔY_, Y_ = backward_inv(ΔY, Y, CH.AN_Y[j])
        ΔX, ΔY, X, Y = backward_inv(ΔX_, ΔY_, X_, Y_, CH.CL[j])
    end
    return ΔX, ΔY, X, Y
end

# Forward pass and compute logdet
function forward_Y(Y, CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=1:depth
        Y_ = CH.AN_Y[j].forward(Y; logdet=false)
        Y = CH.CL[j].forward_Y(Y_)
    end
    return Y
end

# Inverse pass and compute gradients
function inverse_Y(Zy, CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=depth:-1:1
        Zy_ = CH.CL[j].inverse_Y(Zy)
        Zy = CH.AN_Y[j].inverse(Zy_; logdet=false)
    end
    return Zy
end

## Jacobian-related utils

function jacobian(ΔX, ΔY, Δθ, X, Y, CH::NetworkConditionalHINT; logdet=nothing)
    isnothing(logdet) ? logdet = (CH.logdet && ~CH.is_reversed) : logdet = logdet

    depth = length(CH.CL)
    logdet_ = 0f0
    logdet && (GNΔθ = Array{Parameter, 1}(undef, 0))
    n = Int64(length(Δθ)/depth)
    for j=1:depth
        Δθj = Δθ[n*(j-1)+1:n*j]
        if logdet
            ΔX_, X_, logdet1, GNΔθ1 = CH.AN_X[j].jacobian(ΔX, Δθj[1:2], X)
            ΔY_, Y_, logdet2, GNΔθ2 = CH.AN_Y[j].jacobian(ΔY, Δθj[3:4], Y)
            ΔX, ΔY, X, Y, logdet3, GNΔθ3 = CH.CL[j].jacobian(ΔX_, ΔY_, Δθj[5:end], X_, Y_)
            logdet_ += (logdet1 + logdet2 + logdet3)
            GNΔθ = cat(GNΔθ, GNΔθ1, GNΔθ2, GNΔθ3; dims=1)
        else
            ΔX_, X_ = CH.AN_X[j].jacobian(ΔX, Δθj[1:2], X)
            ΔY_, Y_ = CH.AN_Y[j].jacobian(ΔY, Δθj[3:4], Y)
            ΔX, ΔY, X, Y = CH.CL[j].jacobian(ΔX_, ΔY_, Δθj[5:end], X_, Y_)
        end
    end
    logdet ? (return ΔX, ΔY, X, Y, logdet_, GNΔθ) : (return ΔX, ΔY, X, Y)
end

adjointJacobian(ΔZx, ΔZy, Zx, Zy, CH::NetworkConditionalHINT) = backward(ΔZx, ΔZy, Zx, Zy, CH; set_grad=false)

## Other utils

# Clear gradients
function clear_grad!(CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=1:depth
        clear_grad!(CH.AN_X[j])
        clear_grad!(CH.AN_Y[j])
        clear_grad!(CH.CL[j])
    end
end

# Get parameters
function get_params(CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    p = Array{Parameter, 1}(undef, 0)
    for j=1:depth
        p = cat(p, get_params(CH.AN_X[j]); dims=1)
        p = cat(p, get_params(CH.AN_Y[j]); dims=1)
        p = cat(p, get_params(CH.CL[j]); dims=1)
    end
    return p
end

# Set is_reversed flag in full network tree
function tag_as_reversed!(CH::NetworkConditionalHINT, tag::Bool)
    depth = length(CH.CL)
    CH.is_reversed = tag
    for j=1:depth
        tag_as_reversed!(CH.AN_X[j], tag)
        tag_as_reversed!(CH.AN_Y[j], tag)
        tag_as_reversed!(CH.CL[j], tag)
    end
    return CH
end


text:
@info "Test chain of responsibility pattern..."

bank = Bank(100)
paypal = Paypal(300)
bitcoin = Bitcoin(1000)

bank.successor = paypal
paypal.successor = bitcoin

pay(bank, 500)

println(bank.balance)
println(paypal.balance)
println(bitcoin.balance)


text:
"""
    mutable struct Agent

Basic structure which contains the agent based Agent.
"""
mutable struct Agent

    dims::Int

    declaredSymbols::Dict{String,Array{Symbol,1}}
    declaredUpdates::Dict{String,Expr}

    evolve::Function

    function Agent()
        new(0,
            Dict{String,Array{Symbol}}("Local"=>Symbol[],"Identity"=>Symbol[],
                                        "Global"=>Symbol[],"GlobalArray"=>Symbol[],"Medium"=>Symbol[]),
            Dict{String,Expr}(),
            needCompilation)
    end
end

function Base.show(io::IO,abm::Agent)
    print("PARAMETERS\n")
    for i in keys(abm.declaredSymbols)
        if ! isempty(abm.declaredSymbols[i])
            print(i,"\n\t")
            for j in abm.declaredSymbols[i]
                print(" ",j,",")
            end
            print("\n")
        end
    end

    print("\n\nUPDATE RULES\n")
    for i in keys(abm.declaredUpdates)
        if [i for i in abm.declaredUpdates[i].args if typeof(i) != LineNumberNode] != []
            print(i,"\n")
            print(" ",clean(copy(abm.declaredUpdates[i])),"\n\n")
        end
    end
end

"""
    function needCompilation()

Function that ask for compilation before evolve can be used with the current Agent.
"""
function needCompilation(args...;kwargs...)
    error("Agent not yet compiled or modified. Compile it first.")
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using Printf

function write_to_disk(x::Array{Float64,1}, x_domains::Int, y_domains::Int,
                        xcell::Int, ycell::Int, temp1::Float64, filename:: String)
    f = open(filename, "w")

    #add first x-boundary
    for k = 1:ycell*y_domains + 2
        print(f, lstrip(@sprintf("%15.11f", temp1)))
        if (k < ycell*y_domains + 2)
            print(f, " ")
        end
    end
    print(f, "\n")

    #then, add internal cells + y-boundaries
    c = 0
    for k = 1:x_domains
        for m = 1:xcell
            for i = 1:y_domains
                for j = 1:ycell
                    c += 1
                    if (i==1 && j==1)
                        print(f, lstrip(@sprintf("%15.11f", temp1)))
                        print(f, " ")
                    end
                    print(f, lstrip(@sprintf("%15.11f", x[(i-1) * x_domains * xcell * ycell +
                                             (k-1) * xcell * ycell + (j-1) * xcell + m])))
                    print(f, " ")
                    if (i==y_domains && j==ycell)
                        print(f, lstrip(@sprintf("%15.11f", temp1)))
                        print(f, "\n")
                    end
                end
            end
        end
    end

    #add second x-boundary
    for k = 1:ycell*y_domains + 2
        print(f, lstrip(@sprintf("%15.11f", temp1)))
        if (k < ycell*y_domains + 2)
            print(f, " ")
        end
    end

    close(f)
end

max_stars_repo_name: ornl-training/julia-basics

max_issues_repo_name: ornl-training/julia-basics

max_forks_repo_name: ornl-training/julia-basics


text:
function test_poly(s)
    global solver = s

    test_dir = dirname(@__FILE__)

    @testset "Nonconvex Polynomial Models" begin end
end


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 3inch)

plot(dataset("lattice", "singer"), x=:VoicePart, y=:Height,
     Geom.boxplot(suppress_outliers=true, method=[0.1, 0.2, 0.5, 0.8, 0.9]))


text: abstract type TwoDotsModel end


text:
println("*******************************************************")
@testset "blockFGMRES" begin
	@testset "real matrix" begin
		n = 50;
		A  = sprandn(n,n,.1) + SparseMatrixCSC(10.0I, n, n)
		m  = 2;
		D  = Vector(diag(A))
		M2 = x -> D.\x
		rhs = randn(n,m)
		tol = 1e-6;

		# test printing and behaviour for early stopping
		xtt = blockFGMRES(A,rhs ,3,tol=1e-12,maxIter=3,out=2)
		@test xtt[2]==-1

		# test behaviour for zero rhs
		xtt = blockFGMRES(A,0*rhs,3,tol=tol,maxIter=10,out=2)
		@test xtt[2]==-9
		@test all(xtt[1].==0.0)
		@test length(xtt[1])==m*n
		@test eltype(xtt[1])==eltype(rhs)

		println("-----------------------------------")
		x1 = blockFGMRES(A,rhs ,5,tol=tol,maxIter=100,out=1)
		x3 = blockFGMRES(A,rhs,5,tol=tol,maxIter=100,X=randn(size(rhs)),flexible=true)
		x4 = blockFGMRES(A,rhs,5,tol=tol,maxIter=100,M=M2,flexible = true)

		@test norm(A*x1[1]-rhs)/norm(rhs) < tol
		@test norm(A*x3[1]-rhs)/norm(rhs) < tol
		@test norm(A*x4[1]-rhs)/norm(rhs) < tol
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5
	end
	println("*****************************************************************************")
	@testset "complex matrix" begin
		A  = sprandn(100,100,.1) + SparseMatrixCSC(10.0I, 100, 100) + 1im*(sprandn(100,100,.1) + SparseMatrixCSC(10.0I, 100, 100) )
		m  = 3;
		D  = Vector(diag(A))
		M3 = x -> x./D
		rhs = complex(randn(100,m))
		tol = 1e-6;

		# test behaviour for zero rhs
		xtt = blockFGMRES(A,0.0*rhs,5,tol=tol,maxIter=10,out=2)
		@test xtt[2]==-9
		@test all(xtt[1].==0)
		@test length(xtt[1])==100*m
		@test eltype(xtt[1])==eltype(rhs)

		x1 = blockFGMRES(A,rhs ,5,tol=tol,maxIter=100)
		x3 = blockFGMRES(A,rhs,5,tol=tol,maxIter=100,X=randn(size(rhs)),flexible = true)
		x4 = blockFGMRES(A,rhs,5,tol=tol,maxIter=100,M=M3,flexible = true)

		@test norm(A*x1[1]-rhs)/norm(rhs) < tol
		@test norm(A*x3[1]-rhs)/norm(rhs) < tol
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5
	end
end


text:
# Helper type
# n-face of the polytope, i.e., any polytope of lower dimension `N` representing
# its boundary and the polytope itself (for `N` equal to the space dimension `D`)
struct NFace{D}
  anchor::Point{D,Int}
  extrusion::Point{D,Int}
end

"""
    struct ExtrusionPolytope{D} <: Polytope{D}
      extrusion::Point{D,Int}
      # + private fields
    end

Concrete type for polytopes that can be
represented with an "extrusion" tuple. The underlying extrusion is available
in the field `extrusion`. Instances of this type can be obtained with the constructors

- [`Polytope(extrusion::Int...)`](@ref)
- [`ExtrusionPolytope(extrusion::Int...)`](@ref)

"""
struct ExtrusionPolytope{D} <: Polytope{D}
  extrusion::Point{D,Int}
  nfaces::Vector{NFace}
  nfacesdim::Vector{UnitRange{Int}}
  nf_nfs::Vector{Vector{Int}}
  nf_dim::Vector{Vector{UnitRange{Int}}}
end

# Constructors

"""
    Polytope(extrusion::Int...)

Equivalent to `ExtrusionPolytope(extrusion...)`
"""
Polytope(extrusion::Int...) = Polytope(extrusion)

function Polytope(extrusion::NTuple{D,Int}) where D
  ExtrusionPolytope(extrusion)
end

function ExtrusionPolytope(extrusion::NTuple{N,Int}) where N
  ExtrusionPolytope(Point{N,Int}(extrusion))
end

"""
    ExtrusionPolytope(extrusion::Int...)

Generates an `ExtrusionPolytope` from the tuple `extrusion`.
The values in `extrusion` are either equal to the constant
[`HEX_AXIS`](@ref) or the constant [`TET_AXIS`](@ref).

# Examples

Creating a quadrilateral, a triangle, and a wedge

```jldoctest
using InplaceArrays.ReferenceFEs

quad = ExtrusionPolytope(HEX_AXIS,HEX_AXIS)

tri = ExtrusionPolytope(TET_AXIS,TET_AXIS)

wedge = ExtrusionPolytope(TET_AXIS,TET_AXIS,HEX_AXIS)

println(quad == QUAD)
println(tri == TRI)
println(wedge == WEDGE)

# output
true
true
true

```

"""
function ExtrusionPolytope(extrusion::Int...)
  ExtrusionPolytope(extrusion)
end

"""
Constant to be used in order to indicate a hex-like extrusion axis.
"""
const HEX_AXIS = 1

"""
Constant to be used in order to indicate a tet-like extrusion axis.
"""
const TET_AXIS = 2

# Implementation of the interface

get_faces(p::ExtrusionPolytope) = p.nf_nfs

get_dimrange(p::ExtrusionPolytope) = p.nfacesdim

function Polytope{0}(p::ExtrusionPolytope,Dfaceid::Integer)
  @assert Dfaceid in 1:num_vertices(p)
  _vertex()
end

function Polytope{D}(p::ExtrusionPolytope,Dfaceid::Integer)::ExtrusionPolytope{D} where D
  reffaces = _nface_ref_polytopes(p)
  offset = get_offset(p,D)
  id = Dfaceid + offset
  reffaces[id]
end

function Polytope{D}(p::ExtrusionPolytope{D},Dfaceid::Integer) where D
  @assert Dfaceid == 1
  p
end

function (==)(a::ExtrusionPolytope{D},b::ExtrusionPolytope{D}) where D
  a.extrusion == b.extrusion
end

function get_vertex_coordinates(p::ExtrusionPolytope)
  _vertices_coordinates(p)
end

function get_edge_tangents(p::ExtrusionPolytope)
  _edge_tangents(p)
end

function get_facet_normals(p::ExtrusionPolytope)
  n, _ = _face_normals(p)
  n
end

function get_facet_orientations(p::ExtrusionPolytope)
  _, or = _face_normals(p)
  or
end

function get_vertex_permutations(p::ExtrusionPolytope)
  _admissible_permutations(p)
end

# Helpers

function ExtrusionPolytope(extrusion::Point{D,Int}) where D
  zerop = Point{D,Int}(zeros(Int64, D))
  pol_nfs_dim = _polytopenfaces(zerop, extrusion)
  pol_nfs = pol_nfs_dim[1]
  pol_dim = pol_nfs_dim[2]
  nfs_id = Dict(nf => i for (i, nf) in enumerate(pol_nfs))
  nf_nfs_dim = _polytopemesh(pol_nfs, nfs_id)
  nf_nfs = nf_nfs_dim[1]
  nf_dim = nf_nfs_dim[2]
  ExtrusionPolytope{D}(extrusion, pol_nfs, pol_dim, nf_nfs, nf_dim)
end

function ExtrusionPolytope(extrusion::Point{0,Int})
  _vertex()
end

# Generates the array of n-faces of a polytope
function _polytopenfaces(anchor::Point{D,Int}, extrusion::Point{D,Int}) where D
  dnf = _nfdim(extrusion)
  zerop = Point{D,Int}(zeros(Int64, D))
  nf_nfs = []
  nf_nfs = _nfaceboundary!(anchor, zerop, extrusion, true, nf_nfs)
  [sort!(nf_nfs, by = x -> x.anchor[i]) for i = 1:length(extrusion)]
  [sort!(nf_nfs, by = x -> x.extrusion[i]) for i = 1:length(extrusion)]
  [sort!(nf_nfs, by = x -> sum(x.extrusion))]
  numnfs = length(nf_nfs)
  nfsdim = [_nfdim(nf_nfs[i].extrusion) for i = 1:numnfs]
  dimnfs = Array{UnitRange{Int64},1}(undef, dnf + 1)
  dim = 0
  i = 1
  for iface = 1:numnfs
    if (nfsdim[iface] > dim)
      # global dim; # global i
      dim += 1
      dimnfs[dim] = i:iface-1
      i = iface
    end
  end
  dimnfs[dnf+1] = numnfs:numnfs
  return [nf_nfs, dimnfs]
end

_nfdim(a::Point{D,Int}) where D = sum([a[i] > 0 ? 1 : 0 for i = 1:D])

# Generates the list of n-face of a polytope the d-faces for 0 <= d <n on its
# boundary
function _nfaceboundary!(
  anchor::Point{D,Int},
  extrusion::Point{D,Int},
  extend::Point{D,Int},
  isanchor::Bool,
  list) where D

  newext = extend
  list = [list..., NFace{D}(anchor, extrusion)]
  for i = 1:D
    curex = newext[i]
    if (curex > 0) # Perform extension
      func1 = (j -> j == i ? 0 : newext[j])
      newext = Point{D,Int}([func1(i) for i = 1:D])
      func2 = (j -> j == i ? 1 : 0)
      edim = Point{D,Int}([func2(i) for i = 1:D])
      func3 = (j -> j >= i ? anchor[j] : 0)
      tetp = Point{D,Int}([func3(i) for i = 1:D]) + edim
      if (curex == 1) # Quad extension
        list = _nfaceboundary!(anchor + edim, extrusion, newext, false, list)
      elseif (isanchor)
        list = _nfaceboundary!(tetp, extrusion, newext, false, list)
      end
      list = _nfaceboundary!(
        anchor,
        extrusion + edim * curex,
        newext,
        false,
        list)
    end
  end
  return list
end

# Provides for all n-faces of a polytope the d-faces for 0 <= d <n on its
# boundary (e.g., given a face, it provides the set of edges and corners on its
# boundary) using the global n-face numbering of the base polytope
function _polytopemesh(nfaces::Vector{NFace{D}}, nfaceid::Dict) where D
  num_nfs = length(nfaces)
  nfnfs = Vector{Vector{Int64}}(undef, num_nfs)
  nfnfs_dim = Vector{Vector{UnitRange{Int64}}}(undef, num_nfs)
  for (inf, nf) in enumerate(nfaces)
    nfs_dim_nf = _polytopenfaces(nf.anchor, nf.extrusion)
    nf_nfs = nfs_dim_nf[1]
    dimnfs = nfs_dim_nf[2]
    nfnfs[inf] = [get(nfaceid, nf, nf) for nf in nf_nfs]
    nfnfs_dim[inf] = dimnfs
  end
  return [nfnfs, nfnfs_dim]
end

# Returns an array with the reference polytopes for all n-faces (undef for vertices)
function _nface_ref_polytopes(p::ExtrusionPolytope)
  function _eliminate_zeros(a)
    b = Int[]
    for m in a
      if (m != 0)
        push!(b, m)
      end
    end
    return Tuple(b)
  end
  nf_ref_p = Vector{ExtrusionPolytope}(undef, length(p.nfaces))
  ref_nf_ps = ExtrusionPolytope[]
  v = _vertex()
  for (i_nf, nf) in enumerate(p.nfaces)
    r_ext = _eliminate_zeros(nf.extrusion)
    if r_ext != ()
      k = 0
      for (i_p, ref_p) in enumerate(ref_nf_ps)
        if r_ext == ref_p.extrusion
          k = i_p
          nf_ref_p[i_nf] = ref_p
        end
      end
      if k == 0
        ref_p = ExtrusionPolytope(r_ext)
        push!(ref_nf_ps, ref_p)
        k = length(ref_nf_ps) + 1
        nf_ref_p[i_nf] = ref_p
      end
    else
        nf_ref_p[i_nf] = v
    end
  end
  return nf_ref_p
end

function _vertex()
  ext = ()
  nfdim = [[1:1]]
  nfnfs = [[1]]
  nfanc = Point{0,Int}()
  nf = NFace{0}(nfanc,nfanc)
  nfs = [nf]
  return ExtrusionPolytope{0}(ext, nfs, [1:1], nfnfs, nfdim)
end

function _dimfrom_fs_dimto_fs(p::ExtrusionPolytope, dim_from::Int, dim_to::Int)
  @assert dim_to <= dim_from
  dim_from += 1
  dim_to += 1
  dffs_r = p.nf_dim[end][dim_from]
  dffs_dtfs = Vector{Vector{Int}}(undef, dffs_r[end] - dffs_r[1] + 1)
  offs = p.nf_dim[end][dim_to][1] - 1
  for (i_dff, dff) in enumerate(dffs_r)
    dff_nfs = p.nf_nfs[dff]
    dff_dtfs_r = p.nf_dim[dff][dim_to]
    dff_dtfs = dff_nfs[dff_dtfs_r]
    dffs_dtfs[i_dff] = dff_dtfs .- offs
    # @santiagobadia : With or without offset ?
  end
  return dffs_dtfs
end

# It generates the list of coordinates of all vertices in the polytope. It is
# assumed that the polytope has the bounding box [0,1]**dim
function _vertices_coordinates(p::ExtrusionPolytope{D}) where D
  vs = _dimfrom_fs_dimto_fs(p, D, 0)[1]
  vcs = Point{D,Float64}[]
  for i = 1:length(vs)
    cs = convert(Vector{Float64}, [p.nfaces[vs[i]].anchor...])
    push!(vcs, cs)
  end
  return vcs
end

# Return the n-faces vertices coordinates array for a given n-face dimension
function _nfaces_vertices(p,d)
  nc = _num_nfaces(p,d)
  verts = _vertices_coordinates(p)
  faces_vs = _dimfrom_fs_dimto_fs(p,d,0)
  cfvs = collect(LocalToGlobalArray(faces_vs,verts))
end

# Returns number of nfaces of dimension dim
function _num_nfaces(polytope::ExtrusionPolytope, dim::Integer)
  n = length(polytope.nf_dim)
  k = 0
  for nface = 1:n
    d = length(polytope.nf_dim[nface]) - 1
    if d == dim
      k += 1
    end
  end
  k
end

# It generates the outwards normals of the facets of a polytope. It returns two
# arrays, the first one being the outward normal and the second one the orientation.
function _face_normals(p::ExtrusionPolytope{D}) where D
  nf_vs = _dimfrom_fs_dimto_fs(p, D - 1, 0)
  vs = _vertices_coordinates(p)
  f_ns = Point{D,Float64}[]
  f_os = Int[]
  for i_f = 1:length(p.nf_dim[end][end-1])
    n, f_o = _facet_normal(p, nf_vs, vs, i_f)
    push!(f_ns, Point{D,Float64}(n))
    push!(f_os, f_o)
  end
  return f_ns, f_os
end

function _face_normals(p::ExtrusionPolytope{0})
  (VectorValue{0,Float64}[], Int[])
end

function _face_normals(p::ExtrusionPolytope{1})
  (VectorValue{1,Float64}[(-1),(1)], [1,1])
end

# It generates the tangent vectors for polytope edges.
function _edge_tangents(p::ExtrusionPolytope{D}) where D
  ed_vs = _nfaces_vertices(p,1)
  return ts = [(t = vs[2]-vs[1])/norm(t) for vs in ed_vs ]
end

function _edge_tangents(p::ExtrusionPolytope{0})
  VectorValue{0,Float64}[]
end

function _facet_normal(p::ExtrusionPolytope{D}, nf_vs, vs, i_f) where D
  if (length(p.extrusion) > 1)
    v = Float64[]
    for i = 2:length(nf_vs[i_f])
      vi = vs[nf_vs[i_f][i]] - vs[nf_vs[i_f][1]]
      push!(v, vi...)
    end
    n = nullspace(transpose(reshape(v, D, length(nf_vs[i_f]) - 1)))
    n = n .* 1 / sqrt(dot(n, n))
    ext_v = _vertex_not_in_facet(p, i_f, nf_vs)
    v3 = vs[nf_vs[i_f][1]] - vs[ext_v]
    f_or = 1
    if dot(v3, n) < 0.0
      n *= -1
      f_or = -1
    end
  elseif (length(p.extrusion) == 1)
    ext_v = _vertex_not_in_facet(p, i_f, nf_vs)
    n = vs[nf_vs[i_f][1]] - vs[ext_v]
    n = n .* 1 / dot(n, n)
    f_or = 1
  else
    error("O-dim polytopes do not have properly define outward facet normals")
  end
  return n, f_or
end

function _vertex_not_in_facet(p, i_f, nf_vs)
  for i in p.nf_dim[end][1]
    is_in_f = false
    for j in nf_vs[i_f]
      if i == j
        is_in_f = true
        break
      end
    end
    if !is_in_f
      return i
      break
    end
  end
end

# It generates all the admissible permutations of nodes that lead to an
# admissible polytope
function _admissible_permutations(p::ExtrusionPolytope)
  p_dims = length(p.extrusion)
  p_vs = _dimfrom_fs_dimto_fs(p, p_dims, 0)
  vs = p.nfaces[p_vs...]
  num_vs = length(vs)
  ext = p.extrusion
  l = [i for i = 1:num_vs]
  permuted_polytopes = Vector{Int}[]
  for c in Combinatorics.permutations(l)#, p_dims + 1)
    admissible_polytope = true
    c1 = vs[c[1]].anchor
    for j = 2:p_dims+1
      c2 = vs[c[j]].anchor
      if (!_are_nodes_connected(c1, c2, ext))
        admissible_polytope = false
      end
    end
    if (admissible_polytope)
      push!(permuted_polytopes, c)
    end
  end
  return permuted_polytopes
end

# Auxiliary function that determines whether two nodes are connected
function _are_nodes_connected(c1, c2, ext)
  sp_dims = length(c1)
  d = zeros(length(c1))
  for i = 1:length(d)
    d[i] = c2[i] - c1[i]
  end
  dn = sum(d .* d)
  connected = false
  if (dn == 1)
    connected = true
  else
    k = 0
    for j = 1:sp_dims
      if (ext[j] == 2)
        if (c1[j] == 1 || c2[j] == 1)
          k = j
        end
      end
    end
    for l = 1:k
      d[l] = 0
    end
    dn = sum(d .* d)
    if (dn == 0)
      connected = true
    end
  end
  return connected
end

# Some particular cases

"""
    const VERTEX = Polytope()
"""
const VERTEX = Polytope()

"""
    const SEGMENT = Polytope(HEX_AXIS)
"""
const SEGMENT = Polytope(HEX_AXIS)

# TODO use larger names
"""
    const TRI = Polytope(TET_AXIS,TET_AXIS)
"""
const TRI = Polytope(TET_AXIS,TET_AXIS)

"""
    const QUAD = Polytope(HEX_AXIS,HEX_AXIS)
"""
const QUAD = Polytope(HEX_AXIS,HEX_AXIS)

"""
    const TET = Polytope(TET_AXIS,TET_AXIS,TET_AXIS)
"""
const TET = Polytope(TET_AXIS,TET_AXIS,TET_AXIS)

"""
    const HEX = Polytope(HEX_AXIS,HEX_AXIS,HEX_AXIS)
"""
const HEX = Polytope(HEX_AXIS,HEX_AXIS,HEX_AXIS)

"""
    const WEDGE = Polytope(TET_AXIS,TET_AXIS,HEX_AXIS)
"""
const WEDGE = Polytope(TET_AXIS,TET_AXIS,HEX_AXIS)

"""
    const PYRAMID = Polytope(HEX_AXIS,HEX_AXIS,TET_AXIS)
"""
const PYRAMID = Polytope(HEX_AXIS,HEX_AXIS,TET_AXIS)


text:
function to_rna(dna)

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -13.110901270484293 <= q <= 11.380819522292354)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.7773677229963925 + -0.31823111196481735*1/(1 + exp(-(0.13898320039278467 + -0.4573031975765818*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + -0.027363793963717686*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + -0.7118578003946046*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))))))) + -0.05869662098445039*1/(1 + exp(-(0.7177153834386609 + 0.046197114005355644*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + 0.8341427099559482*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + -0.5961814570241306*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))))))) + -0.31920840200501477*1/(1 + exp(-(-0.05621054225465549 + -0.4685096496776473*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + -0.723697578017108*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + 0.643520408356284*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))))))))) + 1/(1 + exp(-(-0.2800668423712529 + 0.5241807567635033*1/(1 + exp(-(0.13898320039278467 + -0.4573031975765818*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + -0.027363793963717686*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + -0.7118578003946046*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))))))) + 0.9218858277065758*1/(1 + exp(-(0.7177153834386609 + 0.046197114005355644*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + 0.8341427099559482*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + -0.5961814570241306*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))))))) + 0.0068491484591501894*1/(1 + exp(-(-0.05621054225465549 + -0.4685096496776473*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + -0.723697578017108*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + 0.643520408356284*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))))))))) + 1/(1 + exp(-(-0.0923150493931284 + -0.18688964829603716*1/(1 + exp(-(0.13898320039278467 + -0.4573031975765818*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + -0.027363793963717686*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + -0.7118578003946046*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))))))) + 0.4699561510424055*1/(1 + exp(-(0.7177153834386609 + 0.046197114005355644*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + 0.8341427099559482*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + -0.5961814570241306*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))))))) + -0.26284635568255643*1/(1 + exp(-(-0.05621054225465549 + -0.4685096496776473*1/(1 + exp(-(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])))) + -0.723697578017108*1/(1 + exp(-(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])))) + 0.643520408356284*1/(1 + exp(-(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
cd(@__DIR__); allNGWPs = false; include("setups/toronto.jl")
gr(dpi = 200)

## pedestrian signal 16 most important VM-NGWP vectors
VM_NGWP = vm_ngwp(𝚽, GP_dual)
fp = load("../datasets/new_toronto.jld", "fp")
G_Sig.f = reshape(fp, (N, 1))
dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
println("================fp-VM-NGWP-top-basis-vectors=================")
for i in 1:16
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    w = VM_NGWP[dr, dc, :]
    j, k, l = NGWP_jkl(GP_dual, dr, dc)
    print("(j, k, l) = ($(j), $(k), $(l))  ")
    if j == jmax
        print("φ_{$(GP_dual.ind[dr]-1)}")
    end
    println()
    gplot(A, X; width = 1)
    scatter_gplot!(X; marker = w, plotOrder = :s2l, ms = 3)
    plt = plot!(cbar = false, clims = (-0.075,0.075))
    savefig(plt,
        "../figs/Toronto_fp_DAG_VM_NGWP_ibv$(lpad(i,2,"0"))_j$(j)_k$(k)_l$(l).png")
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2020, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Basics (Modia3D/Basics/_module.jl)
#

# Epsilon and sign
const neps = sqrt( eps() )
sign_eps(value::Float64; seps::Float64 = 100*neps)::Float64 = value > seps ? 1.0 : (value < -seps ? -1.0 : 0.0)

function normalizeVector(n::SVector{3,Float64})::SVector{3,Float64}
  nabs = norm(n)
  if nabs <= neps
    println("neps ", neps)
    println("nabs ", nabs)
    @assert(nabs > neps) # && norm(vec) > eps()
    # return nothing
  else
    return n/nabs
  end
end

# MVector / MMatrix
zeroMVector()   = MVector{3,Float64}(0.0, 0.0, 0.0)
onesMVector()   = MVector{3,Float64}(1.0, 1.0, 1.0)
nullMRotation() = MMatrix{3,3,Float64,9}(1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0)

# Standard constants
const degToRad      = pi/180.0
const radToDeg      = 180.0/pi
const ZeroMVector   = MVector{3,Float64}(0.0,0.0,0.0)
const NoNameDefined = "_NoNameDefined"

# Standard types
const PositionMVector = MVector{3,Float64}
const RotationMMatrix = MMatrix{3,3,Float64,9}

"""    mutable struct BoundingBox - Smallest box that contains a visual element"""
mutable struct BoundingBox
   x_min::Float64
   x_max::Float64
   y_min::Float64
   y_max::Float64
   z_min::Float64
   z_max::Float64
   BoundingBox() = new(0.0,0.0,0.0,0.0,0.0,0.0)
   BoundingBox(x_min,x_max,y_min,y_max,z_min,z_max) = new(x_min,x_max,y_min,y_max,z_min,z_max)
end

# Time functions
linearMovement(delta_x, tStart, tEnd, time) = delta_x*(time-tStart)/(tEnd-tStart)

# Trailing part of type name
function trailingPartOfTypeAsString(obj)::String
   name = string( typeof(obj) )

   # Determine trailing solid (after last ".")
   i = first(something(findlast(".", name), 0:-1))
   return i > 0 && i < length(name) ? name[i+1:end] : name
end

function trailingPartOfName(name::AbstractString)::String
   # Determine trailing part of name (after last ".")
   i = first(something(findlast(".", name), 0:-1))
   return i > 0 && i < length(name) ? name[i+1:end] : name
end

"""
    readDictOfStructsFromJSON(fileName, StructType)

Read a JSON file from `fileName` and return a `Dict{String, StructType}` dictionary.
`StructType` must be a mutable struct type with a constructor `StructType()`.
"""
function readDictOfStructsFromJSON(fileName, StructType)
    dict1 = JSON.parsefile(fileName)
    palette = Dict{String, StructType}()
    for (key1,value1) in dict1
        obj = StructType()
        for (key2,value2) in value1
            setfield!(obj, Symbol(key2), value2)
        end
        palette[key1] = obj
    end
    return palette
end

"""
    listKeys(dict)

List the keys of dictionary `dict` in sorted order.
"""
function listKeys(dict)
    for key in sort(collect(keys(dict)))
        if typeof(key) == String
            println("    \"$key\"")
        else
            println("   $key")
        end
    end
end

deleteItem(container, item) = deleteat!(container, findall(x->x==item, container))

function remove!(container, items)
    for item in items
        deleteat!(container, findall(x->x==item, container))
    end
end


text:
# Before running this, please make sure to activate and instantiate the
# environment with [this `Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-powergen/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-powergen/Manifest.toml).
# For instance, copy these files to a folder 'EX-powergen', `cd` to it and
#
# ```julia
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# **Main author**: [Geoffroy Dolphin](https://github.com/gd1989)
#
# ## Initial data processing
# In this tutorial we are fitting solar and wind power generation output for Germany using weather data.
# We focus on the use of a simple linear estimator. Let's load the required packages to get started.

using MLJ
using UrlDownload
using PyPlot
import DataFrames: DataFrame, describe, names, select!
using Statistics

LinearRegressor = @load LinearRegressor pkg=MLJLinearModels

# The power generation and weather data come from two separate datasets.
# We downloaded both datasets from [Open Power networks](https://open-power-system-data.org/).
# The power generation data is available [here](https://data.open-power-system-data.org/time_series/) and the weather data is available [here](https://data.open-power-system-data.org/weather_data/).
# Note that the first row in both datasets contains headers; hence we don't need to provide any.

data_repo = "https://raw.githubusercontent.com/tlienart/DataScienceTutorialsData.jl/master/data"

url_power   = data_repo * "/power_syst/DE_power_hourly.csv"
url_weather = data_repo * "/power_syst/DE_weather_data.csv"

power   = DataFrame(urldownload(url_power))
weather = DataFrame(urldownload(url_weather));

# We've loaded both datasets but only have a rough idea of what's in each of them.
# Let's get a quick overview of the power dataset:

describe(power, :mean, :nmissing)

# and the weather dataset (we only show the first 20 rows as there's 150 features):

first(describe(weather, :mean, :nmissing), 20)

# Let's first consider the columns of the power DataFrame.
# For the purpose of this tutorial we are only interested in actual wind and solar generation.
# So we select a subset of the power dataframe containing only `:utc_timestamp, :DE_solar_generation_actual` and `:DE_wind_generation_actual`:

select!(power, [
    :utc_timestamp,
    :DE_solar_generation_actual,
    :DE_wind_generation_actual]);

# Inspection of the column names in the weather DataFrame shows that we have weather data at the regional level (38 NUTS-2 statistical regions), which consists of four different weather variables: `_windspeed_10m, _temperature, _radiation_direct_horizontal, _radiation_diffuse_horizontal`.
# To match the structure of the power data, we need to aggregate at the country level.
# We do this by calculating the country-level average across all 38 regions, for each weather variable.
# This means we need to calculate an average across columns of the DataFrame that refer to the same weather variable.
# To do this we define a simple function to filter columns based on weather variable name and create new dataframes.

colnames = names(weather)

filter_by_name(name, cols) =
    filter(cn -> occursin(name, String(cn)), cols)

wind   = weather[:, filter_by_name("windspeed", colnames)]
temp   = weather[:, filter_by_name("temperature", colnames)]
raddir = weather[:, filter_by_name("radiation_direct", colnames)]
raddif = weather[:, filter_by_name("radiation_diffuse", colnames)];

# Next, we create a new column in each DataFrame to store the country-level mean and calculate the mean.
# We use a nested loop.
# The lower level loop iterates over all rows of a given DataFrame, the higher-level loop iterates over all dataframes in the df Array.

dfs = [wind, temp, raddir, raddif]
col_mean = [:windspeed_mean, :temp_mean, :raddir_mean, :raddif_mean];

# the zip function associates elements of two objects in the same position with one another:

n_rows = size(first(dfs), 1)
for (df, name) in zip(dfs, col_mean)
    df[!, name] = zeros(n_rows)
    for (i, row) in enumerate(eachrow(df))
      df[i, name] = mean(row)
    end
end;

# Now that we have all variables we need to conduct our analysis, let's gather them in a single DataFrame...

data = DataFrame(
    Timestamp     = weather.utc_timestamp,
    Solar_gen     = power.DE_solar_generation_actual,
    Wind_gen      = power.DE_wind_generation_actual,
    Windspeed     = wind.windspeed_mean,
    Temperature   = temp.temp_mean,
    Radiation_dir = raddir.raddir_mean,
    Radiation_dif = raddif.raddif_mean);

# ...and have a look at their summary statistics

describe(data, :mean, :median, :nmissing)

# Note that the `describe()` function provides you with information about missing values for each of the columns.
# Fortunately, there are none.
#
# ### Adjusting the scientific types
#
# Let's check the default scientific type that's currently associated with the data features:

schema(data)

# It is important that the scientific type of the variables corresponds to one of the types allowed for use with the models you are planning to use.
# (For more guidance on this, see the [Scientific Type](https://alan-turing-institute.github.io/DataScienceTutorials.jl/data/scitype/) tutorial.
# The scientific type of both `Wind_gen` and `Solar_gen` is currently `Count`. Let's coerce them to `Continuous`.

coerce!(data, :Wind_gen => Continuous)
coerce!(data, :Solar_gen => Continuous)

schema(data)

# We're now ready to go!

# ## Exploratory Data Analysis

# To get a better understanding of our targets, let's plot their respective distributions.

figure(figsize=(8, 6))
hist(data.Solar_gen, color="blue", edgecolor="white", bins=100,
     density=true, alpha=0.5)
xlabel("Solar power generation (MWh)", fontsize=14)
ylabel("Frequency", fontsize=14)
xticks(fontsize=12)
yticks([0, 1e-3, 2e-3], fontsize=12)

# \figalt{Histogram of the solar power generated}{hist_solar.svg}

# As one might expect, the sun doesn't always shine (and certainly not at night), hence there is a very high proportion of observations whose value is equal or close to 0.
# The distribution of wind power generation looks markedly different

figure(figsize=(8, 6))
hist(data.Wind_gen, color="blue", edgecolor = "white", bins=50,
     density=true, alpha=0.5)
xlabel("Wind power generation (MWh)", fontsize=14)
ylabel("Frequency", fontsize=14)

# \figalt{Histogram of the wind power generated}{hist_wind.svg}

# Finally, before fitting the estimator, we might want to gauge what to expect from them by looking at scatter plots.
# Let's look at solar power first.

fig = figure(figsize=(15, 15))

subplot(221)
scatter(data.Solar_gen, data.Radiation_dir)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Solar radiation - directional", fontsize=14)

subplot(222)
scatter(data.Solar_gen, data.Radiation_dif)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Solar radiation - diffuse", fontsize=14)

subplot(223)
scatter(data.Solar_gen, data.Windspeed)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Wind speed (m/s)", fontsize=14)

subplot(224)
scatter(data.Solar_gen, data.Temperature)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Temperature (C)", fontsize=14)

# @@img-wide \figalt{Solar power scatter plots}{solar_scatter.png} @@

# Then at wind generation

fig = figure(figsize=(15, 15))

subplot(221)
scatter(data.Wind_gen, data.Radiation_dir)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Solar radiation - directional", fontsize=14)

subplot(222)
scatter(data.Wind_gen, data.Radiation_dif)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Solar radiation - diffuse", fontsize=14)

subplot(223)
scatter(data.Wind_gen, data.Windspeed)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Wind speed (m/s)", fontsize=14)

subplot(224)
scatter(data.Wind_gen, data.Temperature)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Temperature (C)", fontsize=14)

# @@img-wide \figalt{Wind power scatter plots}{wind_scatter.png} @@

# As expected, solar power generation shows a strong relationship to solar irradiance while wind power generation denotes a strong relationship to wind speed.
#
# ## Models
#
# Let's fit a linear regression to our data.
# We focus on fitting the wind power generation but the same procedure could be applied for the solar power generation (a good exercise!).

y_wind = data.Wind_gen
X = data[:, [:Windspeed, :Temperature, :Radiation_dir, :Radiation_dif]];

# Next, we partition the data in training and test set; we choose the usual 70-30 split:

train, test = partition(collect(eachindex(y_wind)), 0.7, shuffle=true, rng=5);

# then we instantiate a model and fit it:

linReg = LinearRegressor()
m_linReg = machine(linReg, X, y_wind)
fit!(m_linReg, rows=train);

# ### Model evaluation
#
# We've now fitted the model for wind power generation (`Wind_gen`).
# Let's use it to predict values over the test set and investigate the performance:

y_hat = MLJ.predict(m_linReg, rows=test);

# We can start by visualising the observed and predicted valzes of wind power generation.

figure(figsize=(8, 6))
plot(y_hat, color="blue", label="Predicted")
plot(y_wind[test], color="red", label="Observed")
xlabel("Time", fontsize=14)
ylabel("Power generation", fontsize=14)
xticks([])
yticks(fontsize=12)
xlim(0, 100)
legend(fontsize=14)

# \figalt{Observed vs Predicted}{obs_v_pred.svg}

# Let's look at the RMSE on the training and test sets.

rms(y_wind[train], MLJ.predict(m_linReg, rows=train))

# on the test set...

rms(y_wind[test], y_hat)

# Finally, let's plot the residuals to see if there is any obvious remaining structure in them.

res = y_hat .- y_wind[test];

# Let's look at the stem plot of the residuals to check if there's any structure we might not have picked up:

figure(figsize=(12, 6))
stem(res)
xlim(0, length(res))

# @@img-wide \figalt{Residuals}{residuals.png} @@

# Nothing really stands out, the distribution also looks ok:

figure(figsize=(8, 6))
hist(res, color="blue", edgecolor="white", bins=50,
     density=true, alpha=0.5)

# \figalt{Histogram of the residuals}{hist_residuals.svg}

# We leave it at that for now, I hope you found this tutorial interesting.

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using LinearAlgebra, MixedModels, SparseArrays, Test
using MixedModels: mulαβ!

@testset "mulαβ!" begin
    for (m, p, n, q, k) in (
                            (10, 0.7, 5, 0.3, 15),
                            (100, 0.01, 100, 0.01, 20),
                            (100, 0.1, 100, 0.2, 100),
                           )
        a = sprand(m, n, p)
        b = sprand(n, k, q)
        as = sparse(a')
        bs = sparse(b')
        c = zeros(m, k)
        ab = a * b
        arbt = Array(b')
        aab = Array(a) * Array(b)
        @test aab ≈ mulαβ!(c, a, bs', true, true)
        @test aab ≈ mulαβ!(c, a, bs')
        @test aab ≈ mulαβ!(c, a, arbt')
        @test aab ≈ mulαβ!(c, a, arbt')
        @test aab ≈ mulαβ!(fill!(c, 0.0), a, arbt', true, true)
        @test maximum(abs, mulαβ!(c, a, arbt', -1.0, true)) ≤ sqrt(eps())
        @test maximum(abs.(ab - aab)) < 100*eps()
        @test a*bs' == ab
        @test as'*b == ab
        @test as'*bs' == ab
        f = Diagonal(rand(n))
        @test Array(a*f) == Array(a)*f
        @test Array(f*b) == f*Array(b)
    end
end


text:
using Parquet
using Test

function test_col_cursor(file::String, parcompat::String=joinpath(dirname(@__FILE__), "parquet-compatibility"))
    p = ParFile(joinpath(parcompat, file))
    println("loaded ", file)

    nr = nrows(p)
    cnames = colnames(p)
    for cname in cnames
        rr = 1:nr
        println("\tcolumn ", cname, " rows:", rr)
        println("\tvalue, defn, repn, next idx")
        t1 = time()
        cc = ColCursor(p, rr, cname)
        num_read = 0
        for (v,i) in enumerate(cc)
            val,defn,repn = v
            num_read += 1
        end
        println("\t\t", isnull(val) ? nothing : get(val), ", ", defn, ", ", repn, ", ", i)
        println("\t\tread $num_read values in $(time()-t1) time")
    end
end

function test_juliabuilder_row_cursor(file::String, typename::Symbol, parcompat::String=joinpath(dirname(@__FILE__), "parquet-compatibility"))
    p = ParFile(joinpath(parcompat, file))
    println("loaded ", file)

    t1 = time()
    nr = nrows(p)
    cnames = colnames(p)
    schema(JuliaConverter(Main), p, typename)
    jb = JuliaBuilder(p, getfield(Main, typename))
    rc = RecCursor(p, 1:nr, colnames(p), jb)
    rec = nothing
    for i in rc
        rec = i
    end
    println("\t\tlast record: $rec")
    println("\t\tread $nr records in $(time()-t1) time")
end

function test_col_cursor_all_files()
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            test_col_cursor("parquet-testdata/impala/1.1.1-$encformat/$fname.impala.parquet")
        end
    end
end

function test_juliabuilder_row_cursor_all_files()
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            test_juliabuilder_row_cursor("parquet-testdata/impala/1.1.1-$encformat/$fname.impala.parquet", Symbol(encformat * fname))
        end
    end
end

#test_col_cursor_all_files()
test_juliabuilder_row_cursor_all_files()


text:
function printlist(l)
    tl = typeof(l)
    if tl == JLPList
        println("$(typeof(l.val) <: JLPNil ? "" : l.val.val)(")
        _printlist(l, 1)
        println(")")
    else
        println(l)
    end
    return
end

function _printlist(l,index = 0)
    tl = typeof(l)
    if tl == JLPNil
        index == 0 && return
        print(repeat(" ", index*2 - 2))
        return
    elseif tl == JLPList && typeof(l.car) == JLPList
        print(repeat(" ", index*2))
        println("$(typeof(l.car.val) <: JLPNil ? "" : l.car.val.val)(")
        _printlist(l.car, index+1)
        println(")")
    else
        print(repeat(" ", index*2))
        typeof(l) == JLPList ? _printlist(l.car) : println(l)
    end
    tl == JLPList && _printlist(l.cdr, index)
end

function car(l :: JLPList)
    return l.car
end

function cdr(l :: JLPList)
    return l.cdr
end


text:
println(abs(0))
println(abs(-2.7))
println(abs(2))

max_stars_repo_name: lucianobajr/Hacktobers

max_issues_repo_name: lucianobajr/Hacktobers

max_forks_repo_name: lucianobajr/Hacktobers


text:
const TIME_SERIES_CACHE_SIZE_BYTES = 1024 * 1024

abstract type TimeSeriesCache end

function Base.iterate(cache::TimeSeriesCache, state=nothing)
    if state === nothing
        reset!(cache)
    end

    ta = get_next_time_series_array!(cache)
    if ta === nothing
        reset!(cache)
        return nothing
    end

    # Unlike standard Julia iterators, state is maintained internally.
    return ta, 1
end

Base.length(c::TimeSeriesCache) = _get_num_iterations(c)

"""
Return the TimeSeries.TimeArray starting at timestamp.
Reads from storage if the data is not already in cache.

Timestamps must be read sequentially. Random access may be added in the future.

# Arguments

  - `cache::StaticTimeSeriesCache`: cached instance
  - `timestamp::Dates.DateTime`: starting timestamp for the time series array
"""
function get_time_series_array!(cache::TimeSeriesCache, timestamp::Dates.DateTime)
    next_time = get_next_time(cache)
    if next_time === nothing
        throw(InvalidValue("timestamp = $timestamp is beyond the range of the time series"))
    end

    if timestamp != next_time
        throw(InvalidValue("random access is not currently supported. next = $next_time"))
    end

    if next_time > _get_last_cached_time(cache)
        @debug "get_next_time_series_array! update cache" _group = LOG_GROUP_TIME_SERIES next_time
        _update!(cache)
    else
        @debug "get_next_time_series_array! data is in cache" _group = LOG_GROUP_TIME_SERIES next_time
    end

    len = _get_length(cache)
    ta = get_time_series_array(
        _get_component(cache),
        _get_time_series(cache),
        next_time;
        len=len,
        ignore_scaling_factors=_get_ignore_scaling_factors(cache),
    )
    _increment_next_time!(cache, len)
    _decrement_iterations_remaining!(cache)
    return ta
end

"""
Return the next TimeSeries.TimeArray.

Returns `nothing` when all data has been read. Call [`reset!`](@ref) to restart.
Call [`get_next_time`](@ref) to check the start time.

Reads from storage if the data is not already in cache.

# Arguments

  - `cache::StaticTimeSeriesCache`: cached instance
"""
function get_next_time_series_array!(cache::TimeSeriesCache)
    next_time = get_next_time(cache)
    if next_time === nothing
        return
    end

    return get_time_series_array!(cache, next_time)
end

"""
Return the timestamp for the next read with [`get_next_time_series_array!`](@ref).

Return `nothing` if all data has been read.
"""
function get_next_time(cache::TimeSeriesCache)
    if _get_iterations_remaining(cache) == 0
        return
    end

    return cache.common.next_time[]
end

"""
Reset parameters in order to start reading data from the beginning with
[`get_next_time_series_array!`](@ref)
"""
reset!(cache::TimeSeriesCache) = _reset!(cache.common)

_get_component(c::TimeSeriesCache) = _get_component(c.common)
_get_last_cached_time(c::TimeSeriesCache) = c.common.last_cached_time[]
_get_length_available(c::TimeSeriesCache) = c.common.length_available[]
_set_length_available!(c::TimeSeriesCache, len) = c.common.length_available[] = len
_get_length_remaining(c::TimeSeriesCache) = c.common.length_remaining[]
_decrement_length_remaining!(c::TimeSeriesCache, num) = c.common.length_remaining[] -= num
_get_name(c::TimeSeriesCache) = c.common.name
_get_num_iterations(c::TimeSeriesCache) = c.common.num_iterations
_get_ignore_scaling_factors(c::TimeSeriesCache) = c.common.ignore_scaling_factors
_get_type(c::TimeSeriesCache) = typeof(c.common.ts[])
_get_time_series(c::TimeSeriesCache) = c.common.ts[]
_set_time_series!(c::TimeSeriesCache, ts) = c.common.ts[] = ts
_get_iterations_remaining(c::TimeSeriesCache) = c.common.iterations_remaining[]
_decrement_iterations_remaining!(c::TimeSeriesCache) = c.common.iterations_remaining[] -= 1
_get_resolution(cache::TimeSeriesCache) = get_resolution(_get_time_series(cache))

struct TimeSeriesCacheCommon{T <: TimeSeriesData, U <: InfrastructureSystemsComponent}
    ts::Base.RefValue{T}
    component::U
    name::String
    orig_next_time::Dates.DateTime
    next_time::Base.RefValue{Dates.DateTime}
    last_cached_time::Base.RefValue{Dates.DateTime}
    "Total length"
    len::Int
    "Cached data available to read"
    length_available::Base.RefValue{Int}
    "Length remaining to be read on disk"
    length_remaining::Base.RefValue{Int}
    "Total iterations to traverse all data"
    num_iterations::Int
    iterations_remaining::Base.RefValue{Int}
    ignore_scaling_factors::Bool

    function TimeSeriesCacheCommon(;
        ts,
        component,
        name,
        next_time,
        len,
        num_iterations,
        ignore_scaling_factors,
    )
        return new{typeof(ts), typeof(component)}(
            Ref(ts),
            component,
            name,
            next_time,
            Ref(next_time),
            Ref(next_time - Dates.Minute(1)),
            len,
            Ref(0),
            Ref(len),
            num_iterations,
            Ref(num_iterations),
            ignore_scaling_factors,
        )
    end
end

_get_component(c::TimeSeriesCacheCommon) = c.component

function _reset!(common::TimeSeriesCacheCommon)
    common.next_time[] = common.orig_next_time
    common.last_cached_time[] = common.orig_next_time - Dates.Minute(1)
    common.length_available[] = common.len
    common.length_remaining[] = common.len
    common.iterations_remaining[] = common.num_iterations
    return
end

struct ForecastCache{T <: TimeSeriesData, U <: InfrastructureSystemsComponent} <:
       TimeSeriesCache
    common::TimeSeriesCacheCommon{T, U}
    in_memory_count::Int
    horizon::Int
end

"""
Construct ForecastCache to automatically control caching of forecast data.
Maintains some count of forecast windows in memory based on `cache_size_bytes`.

Call Base.iterate or [`get_next_time_series_array!`](@ref) to retrieve data. Each iteration
will return a TimeSeries.TimeArray covering one forecast window of length ``horizon``.

# Arguments

  - `::Type{T}`: subtype of Forecast
  - `component::InfrastructureSystemsComponent`: component
  - `name::AbstractString`: forecast name
  - `start_time::Union{Nothing, Dates.DateTime} = nothing`: forecast start time
  - `horizon::Union{Nothing, Int} = nothing`: forecast horizon
  - `cache_size_bytes = TIME_SERIES_CACHE_SIZE_BYTES`: maximum size of data to keep in memory
  - `ignore_scaling_factors = false`: controls whether to ignore `scaling_factor_multiplier`
    in the time series instance
"""
function ForecastCache(
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::AbstractString;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    horizon::Union{Nothing, Int}=nothing,
    cache_size_bytes=TIME_SERIES_CACHE_SIZE_BYTES,
    ignore_scaling_factors=false,
) where {T <: Forecast}
    metadata_type = time_series_data_to_metadata(T)
    ts_metadata = get_time_series(metadata_type, component, name)
    initial_timestamp = get_initial_timestamp(ts_metadata)
    if start_time === nothing
        start_time = initial_timestamp
    end
    if horizon === nothing
        horizon = get_horizon(ts_metadata)
    end

    # Get one instance to assess data size.
    ts = get_time_series(
        T,
        component,
        name;
        start_time=start_time,
        len=get_horizon(ts_metadata),
    )
    vals = get_time_series_values(component, ts, start_time, len=get_horizon(ts_metadata))
    row_size = _get_row_size(vals)

    count = get_count(ts_metadata)
    if start_time != initial_timestamp
        count -=
            Dates.Millisecond(start_time - initial_timestamp) ÷
            Dates.Millisecond(get_interval(ts_metadata))
    end

    window_size = row_size * horizon
    in_memory_count = minimum((cache_size_bytes ÷ window_size, count))
    @debug "ForecastCache" _group = LOG_GROUP_TIME_SERIES row_size window_size in_memory_count

    return ForecastCache(
        TimeSeriesCacheCommon(
            ts=ts,
            component=component,
            name=name,
            next_time=start_time,
            len=count,
            num_iterations=count,
            ignore_scaling_factors=ignore_scaling_factors,
        ),
        in_memory_count,
        horizon,
    )
end

function _get_count(cache::ForecastCache)
    return minimum((cache.in_memory_count, _get_length_remaining(cache)))
end

_get_length(c::ForecastCache) = c.horizon

function _update!(cache::ForecastCache)
    if _get_length_remaining(cache) == 0
        throw(ArgumentError("Exceeded time series range"))
    end

    count = _get_count(cache)
    next_time = get_next_time(cache)
    len = _get_length(cache)
    ts = get_time_series(
        _get_type(cache),
        _get_component(cache),
        _get_name(cache);
        start_time=next_time,
        len=len,
        count=count,
    )
    _set_length_available!(cache, len)
    _set_time_series!(cache, ts)
    _set_last_cached_time!(cache, next_time)
    _decrement_length_remaining!(cache, count)
    return
end

function _increment_next_time!(cache::ForecastCache, len)
    cache.common.next_time[] += get_interval(_get_time_series(cache))
    return
end

function _set_last_cached_time!(cache::ForecastCache, next_time)
    interval = get_interval(_get_time_series(cache))
    cache.common.last_cached_time[] = next_time + (cache.in_memory_count - 1) * interval
    return
end

struct StaticTimeSeriesCache{T <: TimeSeriesData, U <: InfrastructureSystemsComponent} <:
       TimeSeriesCache
    common::TimeSeriesCacheCommon{T, U}
    in_memory_rows::Int
end

"""
Construct StaticTimeSeriesCache to automatically control caching of time series data.
Maintains rows of data in memory based on `cache_size_bytes`.

Call Base.iterate or [`get_time_series_array`](@ref) to retrieve data. Each iteration will
return a TimeSeries.TimeArray of size 1.

# Arguments

  - `::Type{T}`: subtype of StaticTimeSeries
  - `component::InfrastructureSystemsComponent`: component
  - `name::AbstractString`: time series name
  - `cache_size_bytes = TIME_SERIES_CACHE_SIZE_BYTES`: maximum size of data to keep in memory
  - `ignore_scaling_factors = false`: controls whether to ignore scaling_factor_multiplier
    in the time series instance
"""
function StaticTimeSeriesCache(
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::AbstractString;
    cache_size_bytes=TIME_SERIES_CACHE_SIZE_BYTES,
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    ignore_scaling_factors=false,
) where {T <: StaticTimeSeries}
    metadata_type = time_series_data_to_metadata(T)
    ts_metadata = get_time_series(metadata_type, component, name)
    initial_timestamp = get_initial_timestamp(ts_metadata)
    if start_time === nothing
        start_time = initial_timestamp
    end

    total_length = length(ts_metadata)
    if start_time != initial_timestamp
        total_length -= (start_time - initial_timestamp) ÷ get_resolution(ts_metadata)
    end

    # Get an instance to assess data size.
    ts = get_time_series(T, component, name; start_time=start_time, len=1)
    vals = get_time_series_values(component, ts, start_time, len=1)
    row_size = _get_row_size(vals)

    if row_size > cache_size_bytes
        @warn "Increasing cache size to $row_size in order to accommodate" row_size
        cache_size_bytes = row_size
    end
    in_memory_rows = minimum((cache_size_bytes ÷ row_size, total_length))

    @debug "StaticTimeSeriesCache" _group = LOG_GROUP_TIME_SERIES total_length in_memory_rows
    return StaticTimeSeriesCache(
        TimeSeriesCacheCommon(
            ts=ts,
            component=component,
            name=name,
            next_time=start_time,
            len=total_length,
            num_iterations=total_length,
            ignore_scaling_factors=ignore_scaling_factors,
        ),
        in_memory_rows,
    )
end

_get_count(c::StaticTimeSeriesCache) = nothing
_get_length(cache::StaticTimeSeriesCache) = 1

function _update!(cache::StaticTimeSeriesCache)
    if _get_length_remaining(cache) == 0
        throw(ArgumentError("Exceeded time series range"))
    end

    len = minimum((cache.in_memory_rows, _get_length_remaining(cache)))
    next_time = get_next_time(cache)
    ts = get_time_series(
        _get_type(cache),
        _get_component(cache),
        _get_name(cache);
        start_time=next_time,
        len=len,
    )
    _set_length_available!(cache, len)
    _set_time_series!(cache, ts)
    _set_last_cached_time!(cache, next_time, len)
    _decrement_length_remaining!(cache, len)
    return
end

function _set_last_cached_time!(cache::StaticTimeSeriesCache, next_time, len)
    cache.common.last_cached_time[] = next_time + (len - 1) * _get_resolution(cache)
    return
end

function _increment_next_time!(cache::StaticTimeSeriesCache, len)
    cache.common.next_time[] += len * get_resolution(_get_time_series(cache))
    return
end

function _get_row_size(vals)
    dims = ndims(vals)
    if dims == 1
        row_size = sizeof(vals[1])
    elseif dims == 2
        row_size = sizeof(vals[1, :])
    elseif dims == 3
        row_size = sizeof(vals[1, :, :])
    else
        error("dims=$dims is not supported")
    end

    return row_size
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct PhaseShiftingTransformer <: ACBranch
        name::String
        available::Bool
        activepower_flow::Float64
        reactivepower_flow::Float64
        arc::Arc
        r::Float64
        x::Float64
        primaryshunt::Float64
        tap::Float64
        α::Float64
        rate::Union{Nothing, Float64}
        services::Vector{Service}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `activepower_flow::Float64`
- `reactivepower_flow::Float64`
- `arc::Arc`
- `r::Float64`: System per-unit value, validation range: (0, 4), action if invalid: error
- `x::Float64`: System per-unit value, validation range: (-2, 4), action if invalid: error
- `primaryshunt::Float64`, validation range: (0, 2), action if invalid: error
- `tap::Float64`, validation range: (0, 2), action if invalid: error
- `α::Float64`, validation range: (-1.571, 1.571), action if invalid: warn
- `rate::Union{Nothing, Float64}`, validation range: (0, nothing), action if invalid: error
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct PhaseShiftingTransformer <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    "System per-unit value"
    r::Float64
    "System per-unit value"
    x::Float64
    primaryshunt::Float64
    tap::Float64
    α::Float64
    rate::Union{Nothing, Float64}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function PhaseShiftingTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    PhaseShiftingTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, services, ext, forecasts, InfrastructureSystemsInternal(), )
end

function PhaseShiftingTransformer(; name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    PhaseShiftingTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, services, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function PhaseShiftingTransformer(::Nothing)
    PhaseShiftingTransformer(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        primaryshunt=0.0,
        tap=1.0,
        α=0.0,
        rate=0.0,
        services=Device[],
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::PhaseShiftingTransformer) = value.name
"""Get PhaseShiftingTransformer available."""
get_available(value::PhaseShiftingTransformer) = value.available
"""Get PhaseShiftingTransformer activepower_flow."""
get_activepower_flow(value::PhaseShiftingTransformer) = value.activepower_flow
"""Get PhaseShiftingTransformer reactivepower_flow."""
get_reactivepower_flow(value::PhaseShiftingTransformer) = value.reactivepower_flow
"""Get PhaseShiftingTransformer arc."""
get_arc(value::PhaseShiftingTransformer) = value.arc
"""Get PhaseShiftingTransformer r."""
get_r(value::PhaseShiftingTransformer) = value.r
"""Get PhaseShiftingTransformer x."""
get_x(value::PhaseShiftingTransformer) = value.x
"""Get PhaseShiftingTransformer primaryshunt."""
get_primaryshunt(value::PhaseShiftingTransformer) = value.primaryshunt
"""Get PhaseShiftingTransformer tap."""
get_tap(value::PhaseShiftingTransformer) = value.tap
"""Get PhaseShiftingTransformer α."""
get_α(value::PhaseShiftingTransformer) = value.α
"""Get PhaseShiftingTransformer rate."""
get_rate(value::PhaseShiftingTransformer) = value.rate
"""Get PhaseShiftingTransformer services."""
get_services(value::PhaseShiftingTransformer) = value.services
"""Get PhaseShiftingTransformer ext."""
get_ext(value::PhaseShiftingTransformer) = value.ext

InfrastructureSystems.get_forecasts(value::PhaseShiftingTransformer) = value.forecasts
"""Get PhaseShiftingTransformer internal."""
get_internal(value::PhaseShiftingTransformer) = value.internal

InfrastructureSystems.set_name!(value::PhaseShiftingTransformer, val::String) = value.name = val
"""Set PhaseShiftingTransformer available."""
set_available!(value::PhaseShiftingTransformer, val::Bool) = value.available = val
"""Set PhaseShiftingTransformer activepower_flow."""
set_activepower_flow!(value::PhaseShiftingTransformer, val::Float64) = value.activepower_flow = val
"""Set PhaseShiftingTransformer reactivepower_flow."""
set_reactivepower_flow!(value::PhaseShiftingTransformer, val::Float64) = value.reactivepower_flow = val
"""Set PhaseShiftingTransformer arc."""
set_arc!(value::PhaseShiftingTransformer, val::Arc) = value.arc = val
"""Set PhaseShiftingTransformer r."""
set_r!(value::PhaseShiftingTransformer, val::Float64) = value.r = val
"""Set PhaseShiftingTransformer x."""
set_x!(value::PhaseShiftingTransformer, val::Float64) = value.x = val
"""Set PhaseShiftingTransformer primaryshunt."""
set_primaryshunt!(value::PhaseShiftingTransformer, val::Float64) = value.primaryshunt = val
"""Set PhaseShiftingTransformer tap."""
set_tap!(value::PhaseShiftingTransformer, val::Float64) = value.tap = val
"""Set PhaseShiftingTransformer α."""
set_α!(value::PhaseShiftingTransformer, val::Float64) = value.α = val
"""Set PhaseShiftingTransformer rate."""
set_rate!(value::PhaseShiftingTransformer, val::Union{Nothing, Float64}) = value.rate = val
"""Set PhaseShiftingTransformer services."""
set_services!(value::PhaseShiftingTransformer, val::Vector{Service}) = value.services = val
"""Set PhaseShiftingTransformer ext."""
set_ext!(value::PhaseShiftingTransformer, val::Dict{String, Any}) = value.ext = val

InfrastructureSystems.set_forecasts!(value::PhaseShiftingTransformer, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set PhaseShiftingTransformer internal."""
set_internal!(value::PhaseShiftingTransformer, val::InfrastructureSystemsInternal) = value.internal = val


text:
#
#--------------------------------------#
export lapl
#--------------------------------------#
"""
 for v,u in H¹₀(Ω)

 (v,-∇² u) = (vx,ux) + (vy,uy)\n
          := a(v,u)\n
           = v' * A * u\n
           = (Q*R'*v)'*A_l*(Q*R'*u)\n
           = v'*R*Q'*A_l*Q*R'*u\n

 implemented as

 R'R * QQ' * A_l * u_loc
 where A_l is

 [Dr]'*[rx sx]'*[B 0]*[rx sx]*[Dr]
 [Ds]  [ry sy]  [0 B] [ry sy] [Ds]

 = [Dr]' * [G11 G12]' * [Dr]
   [Ds]    [G12 G22]    [Ds]

"""
function lapl(u::Array
             ,msh::Mesh)

    @unpack Dr,Ds,G11,G12,G22 = msh

    Au = laplace(u,Dr,Ds,G11,G12,G22)
#   Au = gatherScatter(Au,msh)
#   Au = mask(Au,M)

    return Au
end
#--------------------------------------#
function lapl(u::Array
             ,ν::Array
             ,msh::Mesh)

    # -∇⋅(ν∇u)

    return ν .* lapl(u,msh)
end
#--------------------------------------#
function lapl(u::Array
             ,msh1::Mesh
             ,msh2::Mesh)
    # Dealiased implementation
    return lapl(u,msh1)
end
#--------------------------------------#
function lapl(u,M,Jr,Js,QQtx,QQty,Dr,Ds
             ,G11,G12,G22,mult)

Au = u;

Au=laplace(Au,Jr,Js,Dr,Ds,G11,G12,G22);

#Au = Zygote.hook(d->hmp(d),Au);
Au = Zygote.hook(d->d .* mult,Au);

Au = gatherScatter(Au,QQtx,QQty);
Au = mask(Au,M);

return Au
end
#--------------------------------------#
function laplace(u,Dr,Ds,G11,G12,G22)

ur = ABu([],Dr,u)
us = ABu(Ds,[],u)

wr = @. G11*ur + G12*us
ws = @. G12*ur + G22*us

Au = ABu([],Dr',wr) + ABu(Ds',[],ws)

return Au
end
#--------------------------------------#
function laplace(u,Jr,Js,Dr,Ds
                ,G11,G12,G22)

# dealiased

ur = ABu([],Dr,u)
us = ABu(Ds,[],u)

Jur = ABu(Js,Jr,ur)
Jus = ABu(Js,Jr,us)

vr = @. G11*Jur + G12*Jus
vs = @. G12*Jur + G22*Jus

wr = ABu(Js',Jr',vr)
ws = ABu(Js',Jr',vs)

Au = ABu([],Dr',wr) + ABu(Ds',[],ws)

return Au
end
#--------------------------------------#


text:
function talusplot(σ2::Vector, k=1:length(σ2)-1, args...; kwargs...)
    talus = -diff(log2.(σ2))/2 # log₂(σₖ)-log₂(σₖ₊₁)

    df = DataFrame(k=k, talus=talus[k])

    plot(layer(df, x=:k, y=:talus, Geom.line()),
         layer(df, x=:k, y=:talus, Geom.point()),
         Guide.xlabel("k (Principal Component)"), Guide.ylabel("log₂(σₖ)-log₂(σₖ₊₁)"),
         args...; kwargs...)
end

function talusplot(X::Matrix, k=[], args...; kwargs...)
    K = size(X,2)<=size(X,1) ? Symmetric(X'X) : Symmetric(X*X') # Work with the smaller of the kernel matrices

    # compute σᵢ² and talus plot
    σ2 = eigvals(K)
    sort!(σ2,rev=true) # largest to smallest
    σ2 = σ2[1:findlast(σ2.>1e-9)] # get rid of eigenvalues at the end that might be 0 (or below, due to numerical problems)

    talusplot(σ2, isempty(k) ? (1:length(σ2)-1) : k, args...; kwargs...)
end
talusplot(X::Matrix, k::Integer, args...; kwargs...) = talusplot(X, 1:k, args...; kwargs...)


text:
# Statistics for the DQMC simulation runtime. This has nothin to do with
# measurements/observables.
include("statistics.jl")

# Some tapes, checkerboard stuff (which has been neglected a lot)
include("abstract.jl")

# contains `DQMCParameters` which holds parameters relevant to DQMC and the
# general simulation.
include("parameters.jl")

# There are a bunch of functions that require the DQMC type, but logically fit
# with one of its constituents. For example `initialize_stack` requires
# a ::DQMC to figure out how large a bunch of matrices need to be, but
# logically fits in `stack.jl`.
mutable struct DQMC{
        M <: Model, CB <: Checkerboard, ConfType <: Any, RT <: AbstractRecorder,
        Stack <: AbstractDQMCStack, UTStack <: AbstractDQMCStack,
        US <: AbstractUpdateScheduler
    } <: MonteCarloFlavor

    model::M
    conf::ConfType
    temp_conf::ConfType
    last_sweep::Int

    stack::Stack
    ut_stack::UTStack
    scheduler::US
    parameters::DQMCParameters
    analysis::DQMCAnalysis

    recorder::RT
    thermalization_measurements::Dict{Symbol, AbstractMeasurement}
    measurements::Dict{Symbol, AbstractMeasurement}
    lattice_iterator_cache::LatticeIteratorCache

    function DQMC{M, CB, ConfType, RT, Stack, UTStack, US}(args...) where {
            M <: Model, CB <: Checkerboard, ConfType <: Any,
            RT <: AbstractRecorder,
            Stack <: AbstractDQMCStack, UTStack <: AbstractDQMCStack,
            US <: AbstractUpdateScheduler
        }

        @assert isconcretetype(M)
        @assert isconcretetype(ConfType)
        @assert isconcretetype(Stack)
        @assert isconcretetype(UTStack)
        @assert isconcretetype(RT)
        @assert isconcretetype(US)

        new{M, CB, ConfType, RT, Stack, UTStack, US}(args..., LatticeIteratorCache())
    end
end

# Simplified constructor
function DQMC(
        CB, model::M, conf::ConfType, temp_conf::ConfType, last_sweep,
        stack::Stack, ut_stack::UTStack, scheduler::US,
        parameters, analysis,
        recorder::RT,
        thermalization_measurements, measurements
    ) where {M, ConfType, RT, Stack, UTStack, US}

    DQMC{M, CB, ConfType, RT, Stack, UTStack, US}(
        model, conf, temp_conf, last_sweep, stack, ut_stack,
        scheduler, parameters, analysis, recorder,
        thermalization_measurements, measurements
    )
end

# copy constructor
function DQMC(
        mc::DQMC{x, CBT};
        CB = CBT, model::M = mc.model, conf::ConfType = mc.conf,
        temp_conf::ConfType = mc.conf, last_sweep = mc.last_sweep,
        stack::Stack = mc.stack, ut_stack::UTStack = mc.ut_stack,
        scheduler::US = mc.scheduler, parameters = mc.parameters,
        analysis = mc.analysis, recorder::RT = mc.recorder,
        thermalization_measurements = mc.thermalization_measurements,
        measurements = mc.measurements
    ) where {x, CBT, M, ConfType, RT, Stack, UTStack, US}

    DQMC{M, CB, ConfType, RT, Stack, UTStack, US}(
        model, conf, temp_conf, last_sweep, stack, ut_stack,
        scheduler, parameters, analysis, recorder,
        thermalization_measurements, measurements
    )
end

# Contains `DQMCStack`, `propagate` and related functions. The stack keeps track
# of all the matrices that go into calculating the path integral. `propagate`
# moves imaginary time along, between `l = 1` and `l = slices`.
# This also contains some functions for calculating greens functions.
include("stack.jl")

# Contains the `UnequalTimeStack`, related greens calculation methods and
# iterators (i.e. CombinedGreens). This is used for unequal time measurements.
include("unequal_time_stack.jl")

# Contains functions for computations using the matrices representing a time
# slice. Also contains special methods for Checkerboard decomposed matrices.
include("slice_matrices.jl")

# Contains global update schedulers, the updates themselves as well as some
# functions necessary to calculate and perform those updates.
# include("global.jl")
include("updates/scheduler.jl")
include("updates/local_updates.jl")
include("updates/global_updates.jl")
include("updates/parallel_updates.jl")

# Contains the "main/top-level" simulation code. Specifically construction and
# initialization of the DQMC type, convenience functions, `run!` and `replay!`,
# and local updates.
include("DQMC.jl")

# Contains mandatory and optional method definitions expected to exist for a
# given model.
include("DQMC_mandatory.jl")
include("DQMC_optional.jl")

# saving and loading of DQMC
include("FileIO.jl")

# Contains some greens methods
# some others are in stack.jl and unequal_time_stack.jl
include("greens.jl")

# Contains code related to make measurements. Specifically:
# The overall structure
include("measurements/generic.jl")
# Quick constructers and measurement kernels (i.e. applied Wicks theorem)
include("measurements/measurements.jl")
# Superfluid stiffness stuff
include("measurements/superfluid_stiffness.jl")
# Contains some post processing tools
include("measurements/extensions.jl")
# structs and conversions from the old system
include("measurements/deprecated.jl")