text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
reset_default_graph()
include("nnutil.jl")

testtype = "NeuralNetwork1D"
include("NNTrussPull_Domain.jl")

prop = Dict("name"=> testtype, "rho"=> 0.1, "E"=> 200.0, "B"=> 10.0,
            "sigmaY"=>0.300, "K"=>1/9*200, "A0"=> 1.0, "eta"=> 10.0, "nn"=>post_nn)
elements = []
for i = 1:nx
    elnodes = [i, i+1]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))
end
# domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
@load "Data/domain.jld2" domain
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
assembleMassMatrix!(globdat, domain)
# # # need to update state in domain from globdat
# updateStates!(domain, globdat)

# # # # for debug
# X, Y = prepare_strain_stress_data1D(domain)
# x = constant(X)
# y = squeeze(ae(x, [20,20,20,20,1], "nn"))
# sess = Session(); init(sess)
# close("all")
# ADCME.load(sess, "Data/learned_nn.mat")
# out = run(sess, y)
# plot(X[:,2], out,"+", label="NN")
# plot(X[:,2], Y, ".", label="Exact")
# legend()

T = 0.5
NT = 20
Δt = T/NT

nntype = "ae_scaled"
E = prop["E"]
H0 = zeros(1,1)
H0[1,1] = E

state_history, fext_history = read_data("$(@__DIR__)/Data/1.dat")
loss = DynamicMatLawLoss(domain, globdat, state_history, fext_history, nn,Δt)

variable_scope(nntype) do
    global opt = AdamOptimizer().minimize(loss)
end
sess = Session(); init(sess)
# ADCME.load(sess, "Data/learned_nn.mat")
@show run(sess, loss)

for i = 1:2000
    l, _ = run(sess, [loss, opt])
    @show i,l
end

X, Y = prepare_strain_stress_data1D(domain)
y = squeeze(nn(constant(X[:,1]),constant(X[:,2]),constant(X[:,3])))
close("all")
out = run(sess, y)
plot(X[:,2], out,"+", label="NN")
plot(X[:,2], Y, ".", label="Exact")
legend()


text:
export soundspeed, absorption, waterdensity, reflectioncoef, surfaceloss, doppler, bubbleresonance

"""
$(SIGNATURES)
Compute sound speed in water in m/s, given:
- water `temperature` in °C
- `salinity` in ppm
- `depth` in meters
- void fraction (`voidfrac`) in bubbly water
- sound speed in gas (`cgas`), if `voidfrac` > 0
- ratio of density of water to gas (`reldensity`), if `voidfrac` > 0

Implementation based on Mackenzie (1981), Wood (1964) and Buckingham (1997).
"""
function soundspeed(temperature=27.0, salinity=35.0, depth=10.0; voidfrac=0.0, cgas=340.0, reldensity=1000.0)
  c = 1448.96 + 4.591*temperature - 5.304e-2*temperature^2 + 2.374e-4*temperature^3
  c += 1.340*(salinity-35) + 1.630e-2*depth + 1.675e-7*depth^2
  c += -1.025e-2*temperature*(salinity-35) - 7.139e-13*temperature*depth^3
  if voidfrac > 0.0
    m = √reldensity
    c = 1.0/(1.0/c*√((voidfrac*(c/cgas)^2*m+(1.0-voidfrac)/m)*(voidfrac/m+(1-voidfrac)*m)))
  end
  return c
end

"""
$(SIGNATURES)
Compute volume acoustic absorption coefficient in water, given:
- `frequency` in Hz
- `distance` in meters
- `salinity` in ppm
- water `temperature` in °C
- `depth` in meters
- `pH` of water

Implementation based on the Francois-Garrison model.
"""
function absorption(frequency, distance=1000.0, salinity=35.0, temperature=27.0, depth=10.0, pH=8.1)
  f = frequency/1000.0
  d = distance/1000.0
  c = 1412.0 + 3.21*temperature + 1.19*salinity + 0.0167*depth
  A1 = 8.86/c * 10.0^(0.78*pH-5.0)
  P1 = 1.0
  f1 = 2.8*√(salinity/35.0) * 10.0^(4.0-1245.0/(temperature+273.0))
  A2 = 21.44*salinity/c*(1.0+0.025*temperature)
  P2 = 1.0 - 1.37e-4*depth + 6.2e-9*depth*depth
  f2 = 8.17 * 10^(8.0-1990.0/(temperature+273.0)) / (1.0+0.0018*(salinity-35.0))
  P3 = 1.0 - 3.83e-5*depth + 4.9e-10*depth*depth
  if temperature < 20.0
    A3 = 4.937e-4 - 2.59e-5*temperature + 9.11e-7*temperature^2 - 1.5e-8*temperature^3
  else
    A3 = 3.964e-4 - 1.146e-5*temperature + 1.45e-7*temperature^2 - 6.5e-10*temperature^3
  end
  a = A1*P1*f1*f*f/(f1*f1+f*f) + A2*P2*f2*f*f/(f2*f2+f*f) + A3*P3*f*f
  db2amp(-a*d)
end

"""
$(SIGNATURES)
Compute density of water (kg/m^3), given `temperature` in °C and `salinity` in ppm.

Implementation based on Fofonoff (1985 - IES 80).
"""
function waterdensity(temperature=27, salinity=35)
  t = temperature
  A = 1.001685e-04 + t * (-1.120083e-06 + t * 6.536332e-09)
  A = 999.842594 + t * (6.793952e-02 + t * (-9.095290e-03 + t * A))
  B = 7.6438e-05 + t * (-8.2467e-07 + t * 5.3875e-09)
  B = 0.824493 + t * (-4.0899e-03 + t * B)
  C = -5.72466e-03 + t * (1.0227e-04 - t * 1.6546e-06)
  D = 4.8314e-04
  A + salinity * (B + C*√(salinity) + D*salinity)
end

"""
$(SIGNATURES)
Compute complex reflection coefficient at a fluid-fluid boundary, given:
- angle of incidence `θ` (angle to the surface normal)
- relative density of the reflecting medium to incidence medium `ρᵣ`
- relative sound speed of the reflecting medium to incidence medium `cᵣ`
- dimensionless absorption coefficient `δ`

Implementation based on Brekhovskikh & Lysanov. Dimensionless absorption
coefficient based on APL-UW Technical Report 9407.
"""
function reflectioncoef(θ, ρᵣ, cᵣ, δ=0.0)
  n = Complex(1.0, δ) / cᵣ
  t1 = ρᵣ * cos(θ)
  t2 = n*n - sin(θ)^2
  t3 = √abs(t2) * cis(angle(t2)/2)   # ForwardDiff friendly complex √
  (t1 - t3) / (t1 + t3)
end

"""
$(SIGNATURES)
Compute surface reflection coefficient, given:
- `windspeed` in m/s
- `frequency` in Hz
- angle of incidence `θ` (angle to the surface normal)

Implementation based on the APL-UW Technical Report 9407 II-21.
"""
function surfaceloss(windspeed, frequency, θ)
  β = π/2 - θ
  f = frequency/1000.0
  if windspeed >= 6.0
    a = 1.26e-3/sin(β) * windspeed^1.57 * f^0.85
  else
    a = 1.26e-3/sin(β) * 6^1.57 * f^0.85 * exp(1.2*(windspeed-6.0))
  end
  db2amp(-a)
end

"""
$(SIGNATURES)
Compute Doppler frequency, given relative speed between transmitter and
receiver in m/s. `c` is the nominal sound speed in water.
"""
doppler(speed, frequency, c=soundspeed()) = (1.0+speed/c)*frequency

"""
$(SIGNATURES)
Compute resonance frequency of a freely oscillating has bubble in water, given:
- bubble `radius` in meters
- `depth` of bubble in water in meters
- gas ratio of specific heats 'γ', default: 1.4 (for air)
- atmospheric pressure 'p0', default: 1.013e5
- density of water 'ρ' in kg/m³, default: 1022.476

This ignores surface-tension, thermal, viscous and acoustic damping effects, and the pressure-volume relationship is taken to be adiabatic.
Implementation based on Medwin & Clay (1998).
"""
function bubbleresonance(radius, depth=0.0, γ=1.4, p0=1.013e5, ρ=1022.476)
  g = 9.80665 #acceleration due to gravity
  pₐ = p0 + ρ*g*depth
  1 / (2π * radius) * √(3γ * pₐ/ρ)
end

ext: jl


text:
@testset "Edmonds Karp" begin
    # Construct DiGraph
    flow_graph = Graphs.DiGraph(8)

    # Load custom dataset
    flow_edges = [
        (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
        (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
        (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
    ]

    for fg in (flow_graph, Graphs.DiGraph{UInt8}(flow_graph), Graphs.DiGraph{Int16}(flow_graph))
        capacity_matrix = zeros(Int, 8, 8)
        for e in flow_edges
            u, v, f = e
            Graphs.add_edge!(fg, u, v)
            capacity_matrix[u, v] = f
        end
        residual_graph = @inferred(GraphsFlows.residual(fg))

        # Test with default distances
        @test @inferred(GraphsFlows.edmonds_karp_impl(residual_graph, 1, 8, GraphsFlows.DefaultCapacity(residual_graph)))[1] == 3

        # Test with capacity matrix
        @test @inferred(GraphsFlows.edmonds_karp_impl(residual_graph, 1, 8, capacity_matrix))[1] == 28

        # Test the types of the values returned by fetch_path
        function test_find_path_types(residual_graph, s, t, flow_matrix, capacity_matrix)
            v, P, S, flag = GraphsFlows.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
            @test typeof(P) == Vector{Int}
            @test typeof(S) == Vector{Int}
            @test typeof(flag) == Int
            @test typeof(v) == eltype(residual_graph)
        end

        # Test the value of the flags returned.
        function test_find_path_disconnected(residual_graph, s, t, flow_matrix, capacity_matrix)
            h = copy(residual_graph)
            for dst in collect(Graphs.neighbors(residual_graph, s))
                Graphs.rem_edge!(residual_graph, s, dst)
            end
            v, P, S, flag = GraphsFlows.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
            @test flag == 1
            for dst in collect(Graphs.neighbors(h, t))
                Graphs.rem_edge!(h, t, dst)
            end
            v, P, S, flag = GraphsFlows.fetch_path(h, s, t, flow_matrix, capacity_matrix)
            @test flag == 0
            for i in collect(Graphs.inneighbors(h, t))
                Graphs.rem_edge!(h, i, t)
            end
            v, P, S, flag = GraphsFlows.fetch_path(h, s, t, flow_matrix, capacity_matrix)
            @test flag == 2
        end
        flow_matrix = zeros(Int, Graphs.nv(residual_graph), Graphs.nv(residual_graph))
        test_find_path_types(residual_graph, 1, 8, flow_matrix, capacity_matrix)
        test_find_path_disconnected(residual_graph, 1, 8, flow_matrix, capacity_matrix)
    end
end


text:
using Test

module Dependency
  import Base: convert
  export abst, pabst, convert, convertsave

  abstract type abst end
  abstract type pabst end

  convert(::Type{P}, ::T) where {P <: pabst, T <: abst} =
          getfield(T.name.module, Symbol("Packed$(T.name.name)"))
  convertsave(t) = convert(pabst, t)
end

module Extend
  using Main.Dependency
  import Main.Dependency: convert

  export T1, PackedT1, convertsave

  mutable struct T1 <: abst  end
  mutable struct PackedT1 <: pabst  end
end

using Main.Extend

@testset "Ensure converter types can be run from extending namespaces..." begin

@test convertsave(T1()) == Extend.PackedT1

end

#


text:
# This file is copied from Julia itself

module TestAdjointTranspose

using Test, LinearAlgebra, SparseArrays
using ModifiedArrays

# <--- removed inner constructor tests

# <--- removed outer constructor tests

@testset "Adjoint and Transpose basic AbstractArray functionality" begin
    # vectors and matrices with real scalar eltype, and their adjoints/transposes
    intvec, intmat = [1, 2], [1 2 3; 4 5 6]
    tintvec, tintmat = [1 2], [1 4; 2 5; 3 6]
    @testset "length methods" begin
        @test length(modified_adjoint(intvec)) == length(intvec)
        @test length(modified_adjoint(intmat)) == length(intmat)
        @test length(modified_transpose(intvec)) == length(intvec)
        @test length(modified_transpose(intmat)) == length(intmat)
    end
    @testset "size methods" begin
        @test size(modified_adjoint(intvec)) == (1, length(intvec))
        @test size(modified_adjoint(intmat)) == reverse(size(intmat))
        @test size(modified_transpose(intvec)) == (1, length(intvec))
        @test size(modified_transpose(intmat)) == reverse(size(intmat))
    end
    @testset "indices methods" begin
        @test axes(modified_adjoint(intvec)) == (Base.OneTo(1), Base.OneTo(length(intvec)))
        @test axes(modified_adjoint(intmat)) == reverse(axes(intmat))
        @test axes(modified_transpose(intvec)) == (Base.OneTo(1), Base.OneTo(length(intvec)))
        @test axes(modified_transpose(intmat)) == reverse(axes(intmat))
    end
    @testset "IndexStyle methods" begin
        @test IndexStyle(modified_adjoint(intvec)) == IndexLinear()
        @test IndexStyle(modified_adjoint(intmat)) == IndexCartesian()
        @test IndexStyle(modified_transpose(intvec)) == IndexLinear()
        @test IndexStyle(modified_transpose(intmat)) == IndexCartesian()
    end
    # vectors and matrices with complex scalar eltype, and their adjoints/transposes
    complexintvec, complexintmat = [1im, 2im], [1im 2im 3im; 4im 5im 6im]
    tcomplexintvec, tcomplexintmat = [1im 2im], [1im 4im; 2im 5im; 3im 6im]
    acomplexintvec, acomplexintmat = conj.(tcomplexintvec), conj.(tcomplexintmat)
    # vectors and matrices with real-vector and real-matrix eltype, and their adjoints/transposes
    intvecvec = [[1, 2], [3, 4]]
    tintvecvec = [[[1 2]] [[3 4]]]
    intmatmat = [[[1 2]] [[3  4]] [[ 5  6]];
                 [[7 8]] [[9 10]] [[11 12]]]
    tintmatmat = [[hcat([1, 2])] [hcat([7, 8])];
                  [hcat([3, 4])] [hcat([9, 10])];
                  [hcat([5, 6])] [hcat([11, 12])]]
    # vectors and matrices with complex-vector and complex-matrix eltype, and their adjoints/transposes
    complexintvecvec, complexintmatmat = im .* (intvecvec, intmatmat)
    tcomplexintvecvec, tcomplexintmatmat = im .* (tintvecvec, tintmatmat)
    acomplexintvecvec, acomplexintmatmat = conj.(tcomplexintvecvec), conj.(tcomplexintmatmat)
    @testset "getindex methods, elementary" begin
        # implicitly test elementary definitions, for arrays with concrete real scalar eltype
        @test modified_adjoint(intvec) == tintvec
        @test modified_adjoint(intmat) == tintmat
        @test modified_transpose(intvec) == tintvec
        @test modified_transpose(intmat) == tintmat
        # implicitly test elementary definitions, for arrays with concrete complex scalar eltype
        @test modified_adjoint(complexintvec) == acomplexintvec
        @test modified_adjoint(complexintmat) == acomplexintmat
        @test modified_transpose(complexintvec) == tcomplexintvec
        @test modified_transpose(complexintmat) == tcomplexintmat
        # implicitly test elementary definitions, for arrays with concrete real-array eltype
        @test modified_adjoint(intvecvec) == tintvecvec
        @test modified_adjoint(intmatmat) == tintmatmat
        @test modified_transpose(intvecvec) == tintvecvec
        @test modified_transpose(intmatmat) == tintmatmat
        # implicitly test elementary definitions, for arrays with concrete complex-array type
        @test modified_adjoint(complexintvecvec) == acomplexintvecvec
        @test modified_adjoint(complexintmatmat) == acomplexintmatmat
        @test modified_transpose(complexintvecvec) == tcomplexintvecvec
        @test modified_transpose(complexintmatmat) == tcomplexintmatmat
    end
    @testset "getindex(::AdjOrTransVec, ::Colon, ::AbstractArray{Int}) methods that preserve wrapper type" begin
        # for arrays with concrete scalar eltype
        @test modified_adjoint(intvec)[:, [1, 2]] == modified_adjoint(intvec)
        @test modified_transpose(intvec)[:, [1, 2]] == modified_transpose(intvec)
        @test modified_adjoint(complexintvec)[:, [1, 2]] == modified_adjoint(complexintvec)
        @test modified_transpose(complexintvec)[:, [1, 2]] == modified_transpose(complexintvec)
        # for arrays with concrete array eltype
        @test modified_adjoint(intvecvec)[:, [1, 2]] == modified_adjoint(intvecvec)
        @test modified_transpose(intvecvec)[:, [1, 2]] == modified_transpose(intvecvec)
        @test modified_adjoint(complexintvecvec)[:, [1, 2]] == modified_adjoint(complexintvecvec)
        @test modified_transpose(complexintvecvec)[:, [1, 2]] == modified_transpose(complexintvecvec)
    end
    @testset "getindex(::AdjOrTransVec, ::Colon, ::Colon) methods that preserve wrapper type" begin
        # for arrays with concrete scalar eltype
        @test modified_adjoint(intvec)[:, :] == modified_adjoint(intvec)
        @test modified_transpose(intvec)[:, :] == modified_transpose(intvec)
        @test modified_adjoint(complexintvec)[:, :] == modified_adjoint(complexintvec)
        @test modified_transpose(complexintvec)[:, :] == modified_transpose(complexintvec)
        # for arrays with concrete array elype
        @test modified_adjoint(intvecvec)[:, :] == modified_adjoint(intvecvec)
        @test modified_transpose(intvecvec)[:, :] == modified_transpose(intvecvec)
        @test modified_adjoint(complexintvecvec)[:, :] == modified_adjoint(complexintvecvec)
        @test modified_transpose(complexintvecvec)[:, :] == modified_transpose(complexintvecvec)
    end
    @testset "getindex(::AdjOrTransVec, ::Colon, ::Int) should preserve wrapper type on result entries" begin
        # for arrays with concrete scalar eltype
        @test modified_adjoint(intvec)[:, 2] == intvec[2:2]
        @test modified_transpose(intvec)[:, 2] == intvec[2:2]
        @test modified_adjoint(complexintvec)[:, 2] == conj.(complexintvec[2:2])
        @test modified_transpose(complexintvec)[:, 2] == complexintvec[2:2]
        # for arrays with concrete array eltype
        # @test modified_adjoint(intvecvec)[:, 2] == Adjoint.(intvecvec[2:2])
        # @test modified_transpose(intvecvec)[:, 2] == Transpose.(intvecvec[2:2])
        # @test modified_adjoint(complexintvecvec)[:, 2] == Adjoint.(complexintvecvec[2:2])
        # @test modified_transpose(complexintvecvec)[:, 2] == Transpose.(complexintvecvec[2:2])
    end
    @testset "setindex! methods" begin
        # for vectors with real scalar eltype
        @test (wv = modified_adjoint(copy(intvec));
                wv === setindex!(wv, 3, 2) &&
                 wv == setindex!(copy(tintvec), 3, 1, 2)    )
        @test (wv = modified_transpose(copy(intvec));
                wv === setindex!(wv, 4, 2) &&
                 wv == setindex!(copy(tintvec), 4, 1, 2)    )
        # for matrices with real scalar eltype
        @test (wA = modified_adjoint(copy(intmat));
                wA === setindex!(wA, 7, 3, 1) &&
                 wA == setindex!(copy(tintmat), 7, 3, 1)    )
        @test (wA = modified_transpose(copy(intmat));
                wA === setindex!(wA, 7, 3, 1) &&
                 wA == setindex!(copy(tintmat), 7, 3, 1)    )
        # for vectors with complex scalar eltype
        @test (wz = modified_adjoint(copy(complexintvec));
                wz === setindex!(wz, 3im, 2) &&
                 wz == setindex!(copy(acomplexintvec), 3im, 1, 2)   )
        @test (wz = modified_transpose(copy(complexintvec));
                wz === setindex!(wz, 4im, 2) &&
                 wz == setindex!(copy(tcomplexintvec), 4im, 1, 2)   )
        # for  matrices with complex scalar eltype
        @test (wZ = modified_adjoint(copy(complexintmat));
                wZ === setindex!(wZ, 7im, 3, 1) &&
                 wZ == setindex!(copy(acomplexintmat), 7im, 3, 1)   )
        @test (wZ = modified_transpose(copy(complexintmat));
                wZ === setindex!(wZ, 7im, 3, 1) &&
                 wZ == setindex!(copy(tcomplexintmat), 7im, 3, 1)   )
        # for vectors with concrete real-vector eltype
        @test (wv = modified_adjoint(copy(intvecvec));
                wv === setindex!(wv, modified_adjoint([5, 6]), 2) &&
                 wv == setindex!(copy(tintvecvec), [5 6], 2))
        @test (wv = modified_transpose(copy(intvecvec));
                wv === setindex!(wv, modified_transpose([5, 6]), 2) &&
                 wv == setindex!(copy(tintvecvec), [5 6], 2))
        # for matrices with concrete real-matrix eltype
        @test (wA = modified_adjoint(copy(intmatmat));
                wA === setindex!(wA, modified_adjoint([13 14]), 3, 1) &&
                 wA == setindex!(copy(tintmatmat), hcat([13, 14]), 3, 1))
        @test (wA = modified_transpose(copy(intmatmat));
                wA === setindex!(wA, modified_transpose([13 14]), 3, 1) &&
                 wA == setindex!(copy(tintmatmat), hcat([13, 14]), 3, 1))
        # for vectors with concrete complex-vector eltype
        @test (wz = modified_adjoint(copy(complexintvecvec));
                wz === setindex!(wz, modified_adjoint([5im, 6im]), 2) &&
                 wz == setindex!(copy(acomplexintvecvec), [-5im -6im], 2))
        @test (wz = modified_transpose(copy(complexintvecvec));
                wz === setindex!(wz, modified_transpose([5im, 6im]), 2) &&
                 wz == setindex!(copy(tcomplexintvecvec), [5im 6im], 2))
        # for matrices with concrete complex-matrix eltype
        @test (wZ = modified_adjoint(copy(complexintmatmat));
                wZ === setindex!(wZ, modified_adjoint([13im 14im]), 3, 1) &&
                 wZ == setindex!(copy(acomplexintmatmat), hcat([-13im, -14im]), 3, 1))
        @test (wZ = modified_transpose(copy(complexintmatmat));
                wZ === setindex!(wZ, modified_transpose([13im 14im]), 3, 1) &&
                 wZ == setindex!(copy(tcomplexintmatmat), hcat([13im, 14im]), 3, 1))
    end
end

# @testset "Adjoint and Transpose convert methods that convert underlying storage" begin
#     intvec, intmat = [1, 2], [1 2 3; 4 5 6]
#     @test convert(Adjoint{Float64,Vector{Float64}}, modified_adjoint(intvec))::Adjoint{Float64,Vector{Float64}} == modified_adjoint(intvec)
#     @test convert(Adjoint{Float64,Matrix{Float64}}, modified_adjoint(intmat))::Adjoint{Float64,Matrix{Float64}} == modified_adjoint(intmat)
#     @test convert(Transpose{Float64,Vector{Float64}}, modified_transpose(intvec))::Transpose{Float64,Vector{Float64}} == modified_transpose(intvec)
#     @test convert(Transpose{Float64,Matrix{Float64}}, modified_transpose(intmat))::Transpose{Float64,Matrix{Float64}} == modified_transpose(intmat)
# end

# @testset "Adjoint and Transpose similar methods" begin
#     intvec, intmat = [1, 2], [1 2 3; 4 5 6]
#     # similar with no additional specifications, vector (rewrapping) semantics
#     @test size(similar(modified_adjoint(intvec))::Adjoint{Int,Vector{Int}}) == size(modified_adjoint(intvec))
#     @test size(similar(modified_transpose(intvec))::Transpose{Int,Vector{Int}}) == size(modified_transpose(intvec))
#     # similar with no additional specifications, matrix (no-rewrapping) semantics
#     @test size(similar(modified_adjoint(intmat))::Matrix{Int}) == size(modified_adjoint(intmat))
#     @test size(similar(modified_transpose(intmat))::Matrix{Int}) == size(modified_transpose(intmat))
#     # similar with element type specification, vector (rewrapping) semantics
#     @test size(similar(modified_adjoint(intvec), Float64)::Adjoint{Float64,Vector{Float64}}) == size(modified_adjoint(intvec))
#     @test size(similar(modified_transpose(intvec), Float64)::Transpose{Float64,Vector{Float64}}) == size(modified_transpose(intvec))
#     # similar with element type specification, matrix (no-rewrapping) semantics
#     @test size(similar(modified_adjoint(intmat), Float64)::Matrix{Float64}) == size(modified_adjoint(intmat))
#     @test size(similar(modified_transpose(intmat), Float64)::Matrix{Float64}) == size(modified_transpose(intmat))
#     # similar with element type and arbitrary dims specifications
#     shape = (2, 2, 2)
#     @test size(similar(modified_adjoint(intvec), Float64, shape)::Array{Float64,3}) == shape
#     @test size(similar(modified_adjoint(intmat), Float64, shape)::Array{Float64,3}) == shape
#     @test size(similar(modified_transpose(intvec), Float64, shape)::Array{Float64,3}) == shape
#     @test size(similar(modified_transpose(intmat), Float64, shape)::Array{Float64,3}) == shape
# end

@testset "Adjoint and Transpose parent methods" begin
    intvec, intmat = [1, 2], [1 2 3; 4 5 6]
    @test parent(modified_adjoint(intvec)) === intvec
    @test parent(modified_adjoint(intmat)) === intmat
    @test parent(modified_transpose(intvec)) === intvec
    @test parent(modified_transpose(intmat)) === intmat
end

@testset "Adjoint and Transpose vector vec methods" begin
    intvec = [1, 2]
    @test vec(modified_adjoint(intvec)) == intvec
    @test vec(modified_transpose(intvec)) === intvec
    cvec = [1 + 1im]
    @test vec(cvec')[1] == cvec[1]'
end

# @testset "horizontal concatenation of Adjoint/Transpose-wrapped vectors and Numbers" begin
#     # horizontal concatenation of Adjoint/Transpose-wrapped vectors and Numbers
#     # should preserve the Adjoint/Transpose-wrapper to preserve semantics downstream
#     vec, tvec, avec = [1im, 2im], [1im 2im], [-1im -2im]
#     vecvec = [[1im, 2im], [3im, 4im]]
#     tvecvec = [[[1im 2im]] [[3im 4im]]]
#     avecvec = [[[-1im -2im]] [[-3im -4im]]]
#     # for arrays with concrete scalar eltype
#     @test hcat(modified_adjoint(vec), modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == hcat(avec, avec)
#     @test hcat(modified_adjoint(vec), 1, modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == hcat(avec, 1, avec)
#     @test hcat(modified_transpose(vec), modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == hcat(tvec, tvec)
#     @test hcat(modified_transpose(vec), 1, modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == hcat(tvec, 1, tvec)
#     # for arrays with concrete array eltype
#     @test hcat(modified_adjoint(vecvec), modified_adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == hcat(avecvec, avecvec)
#     @test hcat(modified_transpose(vecvec), modified_transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == hcat(tvecvec, tvecvec)
# end

# @testset "map/broadcast over Adjoint/Transpose-wrapped vectors and Numbers" begin
#     # map and broadcast over Adjoint/Transpose-wrapped vectors and Numbers
#     # should preserve the Adjoint/Transpose-wrapper to preserve semantics downstream
#     vec, tvec, avec = [1im, 2im], [1im 2im], [-1im -2im]
#     vecvec = [[1im, 2im], [3im, 4im]]
#     tvecvec = [[[1im 2im]] [[3im 4im]]]
#     avecvec = [[[-1im -2im]] [[-3im -4im]]]
#     # unary map over wrapped vectors with concrete scalar eltype
#     @test map(-, modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == -avec
#     @test map(-, modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == -tvec
#     # unary map over wrapped vectors with concrete array eltype
#     @test map(-, modified_adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -avecvec
#     @test map(-, modified_transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -tvecvec
#     # binary map over wrapped vectors with concrete scalar eltype
#     @test map(+, modified_adjoint(vec), modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec
#     @test map(+, modified_transpose(vec), modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec
#     # binary map over wrapped vectors with concrete array eltype
#     @test map(+, modified_adjoint(vecvec), modified_adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == avecvec + avecvec
#     @test map(+, modified_transpose(vecvec), modified_transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == tvecvec + tvecvec
#     # unary broadcast over wrapped vectors with concrete scalar eltype
#     @test broadcast(-, modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == -avec
#     @test broadcast(-, modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == -tvec
#     # unary broadcast over wrapped vectors with concrete array eltype
#     @test broadcast(-, modified_adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -avecvec
#     @test broadcast(-, modified_transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -tvecvec
#     # binary broadcast over wrapped vectors with concrete scalar eltype
#     @test broadcast(+, modified_adjoint(vec), modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec
#     @test broadcast(+, modified_transpose(vec), modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec
#     # binary broadcast over wrapped vectors with concrete array eltype
#     @test broadcast(+, modified_adjoint(vecvec), modified_adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == avecvec + avecvec
#     @test broadcast(+, modified_transpose(vecvec), modified_transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == tvecvec + tvecvec
#     # trinary broadcast over wrapped vectors with concrete scalar eltype and numbers
#     @test broadcast(+, modified_adjoint(vec), 1, modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec .+ 1
#     @test broadcast(+, modified_transpose(vec), 1, modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec .+ 1
#     @test broadcast(+, modified_adjoint(vec), 1im, modified_adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec .+ 1im
#     @test broadcast(+, modified_transpose(vec), 1im, modified_transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec .+ 1im
#     # ascertain inference friendliness, ref. https://github.com/JuliaLang/julia/pull/25083#issuecomment-353031641
#     sparsevec = SparseVector([1.0, 2.0, 3.0])
#     @test map(-, modified_adjoint(sparsevec), modified_adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}}
#     @test map(-, modified_transpose(sparsevec), modified_transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}}
#     @test broadcast(-, modified_adjoint(sparsevec), modified_adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}}
#     @test broadcast(-, modified_transpose(sparsevec), modified_transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}}
#     @test broadcast(+, modified_adjoint(sparsevec), 1.0, modified_adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}}
#     @test broadcast(+, modified_transpose(sparsevec), 1.0, modified_transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}}
# end

@testset "Adjoint/Transpose-wrapped vector multiplication" begin
    realvec, realmat = [1, 2, 3], [1 2 3; 4 5 6; 7 8 9]
    complexvec, complexmat = [1im, 2, -3im], [1im 2 3; 4 5 -6im; 7im 8 9]
    # Adjoint/Transpose-vector * vector
    @test modified_adjoint(realvec) * realvec == dot(realvec, realvec)
    @test modified_transpose(realvec) * realvec == dot(realvec, realvec)
    @test modified_adjoint(complexvec) * complexvec == dot(complexvec, complexvec)
    @test modified_transpose(complexvec) * complexvec == dot(conj(complexvec), complexvec)
    # vector * Adjoint/Transpose-vector
    @test realvec * modified_adjoint(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3)))
    @test realvec * modified_transpose(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3)))
    @test complexvec * modified_adjoint(complexvec) == broadcast(*, complexvec, reshape(conj(complexvec), (1, 3)))
    @test complexvec * modified_transpose(complexvec) == broadcast(*, complexvec, reshape(complexvec, (1, 3)))
    # Adjoint/Transpose-vector * matrix
    ## Removed type assertions
    @test (modified_adjoint(realvec) * realmat) ==
        reshape(copy(modified_adjoint(realmat)) * realvec, (1, 3))
    @test (modified_transpose(realvec) * realmat) ==
        reshape(copy(modified_transpose(realmat)) * realvec, (1, 3))
    @test (modified_adjoint(complexvec) * complexmat) ==
        reshape(conj(copy(modified_adjoint(complexmat)) * complexvec), (1, 3))
    @test (modified_transpose(complexvec) * complexmat) ==
        reshape(copy(modified_transpose(complexmat)) * complexvec, (1, 3))
    # Adjoint/Transpose-vector * Adjoint/Transpose-matrix
    @test (modified_adjoint(realvec) * modified_adjoint(realmat)) ==
        reshape(realmat * realvec, (1, 3))
    @test (modified_transpose(realvec) * modified_transpose(realmat)) ==
        reshape(realmat * realvec, (1, 3))
    @test (modified_adjoint(complexvec) * modified_adjoint(complexmat)) ==
        reshape(conj(complexmat * complexvec), (1, 3))
    @test (modified_transpose(complexvec) * modified_transpose(complexmat)) ==
        reshape(complexmat * complexvec, (1, 3))
end

@testset "Adjoint/Transpose-wrapped vector pseudoinversion" begin
    realvec, complexvec = [1, 2, 3, 4], [1im, 2, 3im, 4]
    rowrealvec, rowcomplexvec = reshape(realvec, (1, 4)), reshape(complexvec, (1, 4))
    # pinv(Adjoint/Transpose-vector) should match matrix equivalents
    # TODO tighten type asserts once pinv yields Transpose/Adjoint
    @test pinv(modified_adjoint(realvec))::Vector{Float64} ≈ pinv(rowrealvec)
    @test pinv(modified_transpose(realvec))::Vector{Float64} ≈ pinv(rowrealvec)
    @test pinv(modified_adjoint(complexvec))::Vector{Complex{Float64}} ≈ pinv(conj(rowcomplexvec))
    @test pinv(modified_transpose(complexvec))::Vector{Complex{Float64}} ≈ pinv(rowcomplexvec)
end

@testset "Adjoint/Transpose-wrapped vector left-division" begin
    realvec, complexvec = [1., 2., 3., 4.,], [1.0im, 2., 3.0im, 4.]
    rowrealvec, rowcomplexvec = reshape(realvec, (1, 4)), reshape(complexvec, (1, 4))
    # \(Adjoint/Transpose-vector, Adjoint/Transpose-vector) should mat matrix equivalents
    @test modified_adjoint(realvec)\modified_adjoint(realvec) ≈ rowrealvec\rowrealvec
    @test modified_transpose(realvec)\modified_transpose(realvec) ≈ rowrealvec\rowrealvec
    @test modified_adjoint(complexvec)\modified_adjoint(complexvec) ≈ conj(rowcomplexvec)\conj(rowcomplexvec)
    @test modified_transpose(complexvec)\modified_transpose(complexvec) ≈ rowcomplexvec\rowcomplexvec
end

@testset "Adjoint/Transpose-wrapped vector right-division" begin
    realvec, realmat = [1, 2, 3], [1 0 0; 0 2 0; 0 0 3]
    complexvec, complexmat = [1im, 2, -3im], [2im 0 0; 0 3 0; 0 0 -5im]
    rowrealvec, rowcomplexvec = reshape(realvec, (1, 3)), reshape(complexvec, (1, 3))
    # /(Adjoint/Transpose-vector, matrix)
    @test (modified_adjoint(realvec) / realmat) ≈ rowrealvec / realmat
    @test (modified_adjoint(complexvec) / complexmat) ≈ conj(rowcomplexvec) / complexmat
    @test (modified_transpose(realvec) / realmat) ≈ rowrealvec / realmat
    @test (modified_transpose(complexvec) / complexmat) ≈ rowcomplexvec / complexmat
    # /(Adjoint/Transpose-vector, Adjoint matrix)
    @test (modified_adjoint(realvec) / modified_adjoint(realmat)) ≈ rowrealvec / copy(modified_adjoint(realmat))
    @test (modified_adjoint(complexvec) / modified_adjoint(complexmat)) ≈ conj(rowcomplexvec) / copy(modified_adjoint(complexmat))
    @test (modified_transpose(realvec) / modified_adjoint(realmat)) ≈ rowrealvec / copy(modified_adjoint(realmat))
    @test (modified_transpose(complexvec) / modified_adjoint(complexmat)) ≈ rowcomplexvec / copy(modified_adjoint(complexmat))
    # /(Adjoint/Transpose-vector, Transpose matrix)
    @test (modified_adjoint(realvec) / modified_transpose(realmat)) ≈ rowrealvec / copy(modified_transpose(realmat))
    @test (modified_adjoint(complexvec) / modified_transpose(complexmat)) ≈ conj(rowcomplexvec) / copy(modified_transpose(complexmat))
    @test (modified_transpose(realvec) / modified_transpose(realmat)) ≈ rowrealvec / copy(modified_transpose(realmat))
    @test (modified_transpose(complexvec) / modified_transpose(complexmat)) ≈ rowcomplexvec / copy(modified_transpose(complexmat))
end

@testset "norm and opnorm of Adjoint/Transpose-wrapped vectors" begin
    # definitions are in base/linalg/generic.jl
    realvec, complexvec = [3, -4], [3im, -4im]
    # one norm result should be sum(abs.(realvec)) == 7
    # two norm result should be sqrt(sum(abs.(realvec))) == 5
    # inf norm result should be maximum(abs.(realvec)) == 4
    for v in (realvec, complexvec)
        @test norm(modified_adjoint(v)) ≈ 5
        @test norm(modified_adjoint(v), 1) ≈ 7
        @test norm(modified_adjoint(v), Inf) ≈ 4
        @test norm(modified_transpose(v)) ≈ 5
        @test norm(modified_transpose(v), 1) ≈ 7
        @test norm(modified_transpose(v), Inf) ≈ 4
    end
    # one opnorm result should be maximum(abs.(realvec)) == 4
    # two opnorm result should be sqrt(sum(abs.(realvec))) == 5
    # inf opnorm result should be sum(abs.(realvec)) == 7
    for v in (realvec, complexvec)
        @test opnorm(modified_adjoint(v)) ≈ 5
        @test opnorm(modified_adjoint(v), 1) ≈ 4
        @test opnorm(modified_adjoint(v), Inf) ≈ 7
        @test opnorm(modified_transpose(v)) ≈ 5
        @test opnorm(modified_transpose(v), 1) ≈ 4
        @test opnorm(modified_transpose(v), Inf) ≈ 7
    end
end

@testset "adjoint and transpose of Numbers" begin
    @test adjoint(1) == 1
    @test adjoint(1.0) == 1.0
    @test adjoint(1im) == -1im
    @test adjoint(1.0im) == -1.0im
    @test transpose(1) == 1
    @test transpose(1.0) == 1.0
    @test transpose(1im) == 1im
    @test transpose(1.0im) == 1.0im
end

@testset "adjoint!(a, b) return a" begin
    a = fill(1.0+im, 5)
    b = fill(1.0+im, 1, 5)
    @test adjoint!(a, b) === a
    @test adjoint!(b, a) === b
end

@testset "aliasing with adjoint and transpose" begin
    A = collect(reshape(1:25, 5, 5)) .+ rand.().*im
    B = copy(A)
    B .= B'
    @test B == A'
    B = copy(A)
    B .= transpose(B)
    @test B == transpose(A)
    B = copy(A)
    B .= B .* B'
    @test B == A .* A'
end

@testset "test show methods for $t of Factorizations" for t in (Adjoint, Transpose)
    A = randn(4, 4)
    F = lu(A)
    Fop = t(F)
    @test "LinearAlgebra."*sprint(show, Fop) ==
                  "$t of "*sprint(show, parent(Fop))
    @test "LinearAlgebra."*sprint((io, t) -> show(io, MIME"text/plain"(), t), Fop) ==
                  "$t of "*sprint((io, t) -> show(io, MIME"text/plain"(), t), parent(Fop))
end

const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl"))
using .Main.OffsetArrays

@testset "offset axes" begin
    s = Base.Slice(-3:3)'
    @test axes(s) === (Base.OneTo(1), Base.IdentityUnitRange(-3:3))
    @test collect(LinearIndices(s)) == reshape(1:7, 1, 7)
    @test collect(CartesianIndices(s)) == reshape([CartesianIndex(1,i) for i = -3:3], 1, 7)
    @test s[1] == -3
    @test s[7] ==  3
    @test s[4] ==  0
    @test_throws BoundsError s[0]
    @test_throws BoundsError s[8]
    @test s[1,-3] == -3
    @test s[1, 3] ==  3
    @test s[1, 0] ==  0
    @test_throws BoundsError s[1,-4]
    @test_throws BoundsError s[1, 4]
end

end # module TestAdjointTranspose


text:
#!/usr/bin/env julia

function write_toml(f::Function, names::String...)
    path = joinpath(names...) * ".toml"
    mkpath(dirname(path))
    open(path, "w") do io
        f(io)
    end
end

toml_key(str::String) = ismatch(r"[^\w-]", str) ? repr(str) : str
toml_key(strs::String...) = join(map(toml_key, [strs...]), '.')

prefix = joinpath(homedir(), ".julia", "registries", "Uncurated")

write_toml(prefix, "registry") do io
    repo = "https://github.com/JuliaRegistries/Uncurated.git"
    uuid = string(uuid5(uuid_registry, repo))
    println(io, "name = ", repr("Uncurated"))
    println(io, "uuid = ", repr(uuid))
    println(io, "repo = ", repr(repo))
    println(io, "\ndescription = \"\"\"")
    print(io, """
        Official uncurated Julia package registry where people can
        register any package they want without too much debate about
        naming and without enforced standards on documentation or
        testing. We nevertheless encourage documentation, testing and
        some amount of consideration when choosing package names.
        """)
    println(io, "\"\"\"")
    println(io, "\n[packages]")
    for (pkg, p) in sort!(collect(pkgs), by=(p->p.uuid.value)∘last)
        bucket = string(uppercase(first(pkg)))
        path = joinpath(bucket, pkg)
        println(io, p.uuid, " = { name = ", repr(pkg), ", path = ", repr(path), " }")
    end
end

buckets = Dict()
for (pkg, p) in pkgs
    bucket = string(uppercase(first(pkg)))
    push!(get!(buckets, bucket, []), (pkg, p))
end

include("utils.jl")

for (bucket, b_pkgs) in buckets
    for (pkg, p) in b_pkgs
        url = p.url
        startswith(url, "git://github.com") && (url = "https"*url[4:end])
        write_toml(prefix, bucket, pkg, "package") do io
            println(io, "name = ", repr(pkg))
            println(io, "uuid = ", repr(string(p.uuid)))
            println(io, "repo = ", repr(url))
        end
        write_toml(prefix, bucket, pkg, "versions") do io
            for (i, (ver, v)) in enumerate(sort!(collect(p.versions), by=first))
                i > 1 && println(io)
                println(io, "[", toml_key(string(ver)), "]")
                println(io, "hash-sha1 = ", repr(v.sha1))
            end
        end
        verv = filter(vv->!isempty(vv[end].requires), collect(p.versions))
        !isempty(verv) && write_toml(prefix, bucket, pkg, "requirements") do io
            for (i, (ver, v)) in enumerate(sort!(verv, by=first))
                i > 1 && println(io)
                println(io, "[", toml_key(string(ver)), "]")
                for (req, r) in sort!(collect(v.requires), by=first, lt=packagelt)
                    vers = compress_versions(r.versions, collect(keys(pkgs[req].versions)))
                    println(io, toml_key(req), " = ", versions_repr(vers))
                end
            end
        end
    end
end


text:
## R code 14.12
m14.1 <- ulam(
    alist(
        wait ~ normal( mu , sigma ),
        mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
        c(a_cafe,b_cafe)[cafe] ~ multi_normal( c(a,b) , Rho , sigma_cafe ),
        a ~ normal(5,2),
        b ~ normal(-1,0.5),
        sigma_cafe ~ exponential(1),
        sigma ~ exponential(1),
        Rho ~ lkj_corr(2)
    ) , data=d , chains=4 , cores=4 )

## R code 14.18
library(rethinking)
data(chimpanzees)
d <- chimpanzees
d$block_id <- d$block
d$treatment <- 1L + d$prosoc_left + 2L*d$condition

dat <- list(
    L = d$pulled_left,
    tid = d$treatment,
    actor = d$actor,
    block_id = as.integer(d$block_id) )

m14.2 <- ulam(
    alist(
        L ~ binomial(1,p),
        logit(p) <- g[tid] + alpha[actor,tid] + beta[block_id,tid],

        # adaptive priors
        vector[4]:alpha[actor] ~ multi_normal(0,Rho_actor,sigma_actor),
        vector[4]:beta[block_id] ~ multi_normal(0,Rho_block,sigma_block),

        # fixed priors
        g[tid] ~ dnorm(0,1),
        sigma_actor ~ dexp(1),
        Rho_actor ~ dlkjcorr(4),
        sigma_block ~ dexp(1),
        Rho_block ~ dlkjcorr(4)
    ) , data=dat , chains=4 , cores=4 )

## R code 14.19
m14.3 <- ulam(
    alist(
        L ~ binomial(1,p),
        logit(p) <- g[tid] + alpha[actor,tid] + beta[block_id,tid],

        # adaptive priors - non-centered
        transpars> matrix[actor,4]:alpha <-
                compose_noncentered( sigma_actor , L_Rho_actor , z_actor ),
        transpars> matrix[block_id,4]:beta <-
                compose_noncentered( sigma_block , L_Rho_block , z_block ),
        matrix[4,actor]:z_actor ~ normal( 0 , 1 ),
        matrix[4,block_id]:z_block ~ normal( 0 , 1 ),

        # fixed priors
        g[tid] ~ normal(0,1),
        vector[4]:sigma_actor ~ dexp(1),
        cholesky_factor_corr[4]:L_Rho_actor ~ lkj_corr_cholesky( 2 ),
        vector[4]:sigma_block ~ dexp(1),
        cholesky_factor_corr[4]:L_Rho_block ~ lkj_corr_cholesky( 2 ),

        # compute ordinary correlation matrixes from Cholesky factors
        gq> matrix[4,4]:Rho_actor <<- multiply_lower_tri_self_transpose(L_Rho_actor),
        gq> matrix[4,4]:Rho_block <<- multiply_lower_tri_self_transpose(L_Rho_block)
    ) , data=dat , chains=4 , cores=4 , log_lik=TRUE )

## R code 14.24
m14.4 <- ulam(
    alist(
        W ~ dnorm( mu , sigma ),
        mu <- aW + bEW*E,
        aW ~ dnorm( 0 , 0.2 ),
        bEW ~ dnorm( 0 , 0.5 ),
        sigma ~ dexp( 1 )
    ) , data=dat_sim , chains=4 , cores=4 )
precis( m14.4 )

## R code 14.25
m14.5 <- ulam(
    alist(
        c(W,E) ~ multi_normal( c(muW,muE) , Rho , Sigma ),
        muW <- aW + bEW*E,
        muE <- aE + bQE*Q,
        c(aW,aE) ~ normal( 0 , 0.2 ),
        c(bEW,bQE) ~ normal( 0 , 0.5 ),
        Rho ~ lkj_corr( 2 ),
        Sigma ~ exponential( 1 )
    ), data=dat_sim , chains=4 , cores=4 )

## R code 14.30
kl_data <- list(
    N = nrow(kl_dyads),
    N_households = max(kl_dyads$hidB),
    did = kl_dyads$did,
    hidA = kl_dyads$hidA,
    hidB = kl_dyads$hidB,
    giftsAB = kl_dyads$giftsAB,
    giftsBA = kl_dyads$giftsBA
)

m14.4 <- ulam(
    alist(
        giftsAB ~ poisson( lambdaAB ),
        giftsBA ~ poisson( lambdaBA ),
        log(lambdaAB) <- a + gr[hidA,1] + gr[hidB,2] + d[did,1] ,
        log(lambdaBA) <- a + gr[hidB,1] + gr[hidA,2] + d[did,2] ,
        a ~ normal(0,1),

       ## gr matrix of varying effects
        vector[2]:gr[N_households] ~ multi_normal(0,Rho_gr,sigma_gr),
        Rho_gr ~ lkj_corr(4),
        sigma_gr ~ exponential(1),

       ## dyad effects
        transpars> matrix[N,2]:d <-
                compose_noncentered( rep_vector(sigma_d,2) , L_Rho_d , z ),
        matrix[2,N]:z ~ normal( 0 , 1 ),
        cholesky_factor_corr[2]:L_Rho_d ~ lkj_corr_cholesky( 8 ),
        sigma_d ~ exponential(1),

       ## compute correlation matrix for dyads
        gq> matrix[2,2]:Rho_d <<- Chol_to_Corr( L_Rho_d )
    ), data=kl_data , chains=4 , cores=4 , iter=2000 )

## R code 14.38
data(Kline2) # load the ordinary data, now with coordinates
d <- Kline2
d$society <- 1:10 # index observations

dat_list <- list(
    T = d$total_tools,
    P = d$population,
    society = d$society,
    Dmat=islandsDistMatrix )

m14.7 <- ulam(
    alist(
        T ~ dpois(lambda),
        lambda <- (a*P^b/g)*exp(k[society]),
        vector[10]:k ~ multi_normal( 0 , SIGMA ),
        matrix[10,10]:SIGMA <- cov_GPL2( Dmat , etasq , rhosq , 0.01 ),
        c(a,b,g) ~ dexp( 1 ),
        etasq ~ dexp( 2 ),
        rhosq ~ dexp( 0.5 )
    ), data=dat_list , chains=4 , cores=4 , iter=2000 )

## R code 14.45
m14.7nc <- ulam(
    alist(
        T ~ dpois(lambda),
        lambda <- (a*P^b/g)*exp(k[society]),

        # non-centered Gaussian Process prior
        transpars> vector[10]: k <<- L_SIGMA * z,
        vector[10]: z ~ normal( 0 , 1 ),
        transpars> matrix[10,10]: L_SIGMA <<- cholesky_decompose( SIGMA ),
        transpars> matrix[10,10]: SIGMA <- cov_GPL2( Dmat , etasq , rhosq , 0.01 ),

        c(a,b,g) ~ dexp( 1 ),
        etasq ~ dexp( 2 ),
        rhosq ~ dexp( 0.5 )
    ), data=dat_list , chains=4 , cores=4 , iter=2000 )

## R code 14.48
dat_list <- list(
    N_spp = nrow(dstan),
    M = standardize(log(dstan$body)),
    B = standardize(log(dstan$brain)),
    G = standardize(log(dstan$group_size)),
    Imat = diag( nrow(dstan) )
)

m14.8 <- ulam(
    alist(
        B ~ multi_normal( mu , SIGMA ),
        mu <- a + bM*M + bG*G,
        matrix[N_spp,N_spp]: SIGMA <- Imat * sigma_sq,
        a ~ normal( 0 , 1 ),
        c(bM,bG) ~ normal( 0 , 0.5 ),
        sigma_sq ~ exponential( 1 )
    ), data=dat_list , chains=4 , cores=4 )

# Brownian motion model
m14.9 <- ulam(
    alist(
        B ~ multi_normal( mu , SIGMA ),
        mu <- a + bM*M + bG*G,
        matrix[N_spp,N_spp]: SIGMA <- R * sigma_sq,
        a ~ normal( 0 , 1 ),
        c(bM,bG) ~ normal( 0 , 0.5 ),
        sigma_sq ~ exponential( 1 )
    ), data=dat_list , chains=4 , cores=4 )
precis( m14.9 )

## R code 14.51
# add scaled and reordered distance matrix
dat_list$Dmat <- Dmat[ spp_obs , spp_obs ] / max(Dmat)

m14.10 <- ulam(
    alist(
        B ~ multi_normal( mu , SIGMA ),
        mu <- a + bM*M + bG*G,
        matrix[N_spp,N_spp]: SIGMA <- cov_GPL1( Dmat , etasq , rhosq , 0.01 ),
        a ~ normal(0,1),
        c(bM,bG) ~ normal(0,0.5),
        etasq ~ half_normal(1,0.25),
        rhosq ~ half_normal(3,0.25)
    ), data=dat_list , chains=4 , cores=4 )


text:
local_path(args...) = normpath(@__DIR__, "..", args...)
filename(i::Int) = string("ref", i, i in Plots._animation_examples ? ".gif" : ".png")

Plots.theme(:default)

function generate_reference_images(be::Symbol, overwrite::Bool = true)
    for i in eachindex(Plots._examples)
        i in Plots._backend_skips[be] || generate_reference_image(be, i, overwrite)
    end
end

function generate_reference_image(be::Symbol, i::Int, overwrite::Bool = true)
    Plots.backend(be)
    map(ex -> Base.eval(Main, ex), Plots._examples[i].exprs)

    dir = local_path("PlotDocs", string(be))
    isdir(dir) || mkpath(dir)

    fn = joinpath(dir, filename(i))
    overwrite || !isfile(fn) || return nothing

    if i in _animation_inds
        anim = @eval Main anim
        Plots.gif(anim, fn, fps = 15)
    else
        Plots.savefig(fn)
    end

    return nothing
end

_animation_inds = (2, 30)

################################################################################

# Alternative implementation of the script PlotDocs/generate_doc_images.jl
# within the module PlotReferanceImages.

# Usage:
# ```julia
# using PlotReferenceImages, StatPlots, RDatasets, ProgressMeter, DataFrames, Distributions, StatsBase, Iterators
# generate_doc_images()
# ```

# Currently this implementation has issues with the plotlyjs backend.

#
# _doc_images = Dict(
#     "lorenz_attractor" => [:(begin
#         using Plots
#         gr()
#         # define the Lorenz attractor
#         mutable struct Lorenz
#         dt; σ; ρ; β; x; y; z
#         end
#
#         function step!(l::Lorenz)
#         dx = l.σ*(l.y - l.x)       ; l.x += l.dt * dx
#         dy = l.x*(l.ρ - l.z) - l.y ; l.y += l.dt * dy
#         dz = l.x*l.y - l.β*l.z     ; l.z += l.dt * dz
#         end
#
#         attractor = Lorenz((dt = 0.02, σ = 10., ρ = 28., β = 8//3, x = 1., y = 1., z = 1.)...)
#
#
#         # initialize a 3D plot with 1 empty series
#         plt = plot3d(1, xlim=(-25,25), ylim=(-25,25), zlim=(0,50),
#                     title = "Lorenz Attractor", marker = 2)
#
#         # build an animated gif by pushing new points to the plot, saving every 10th frame
#         anim = @gif for i=1:1500
#         step!(attractor)
#         push!(plt, attractor.x, attractor.y, attractor.z)
#         end every 10
#     end)],
#     "waves" => [:(begin
#     using Plots, ProgressMeter
#         pyplot(leg=false, ticks=nothing) #change to the pyplot backend and define some defaults
#         x = y = range(-5, stop = 5, length = 40)
#         zs = zeros(0,40)
#         n = 100
#
#         # create a progress bar for tracking the animation generation
#         prog = Progress(n,1)
#
#         anim = @gif for i in range(0, stop = 2π, length = n)
#         f(x,y) = sin(x + 10sin(i)) + cos(y)
#
#         # create a plot with 3 subplots and a custom layout
#         l = @layout [a{0.7w} b; c{0.2h}]
#         p = plot(x, y, f, st = [:surface, :contourf], layout=l)
#
#         # induce a slight oscillating camera angle sweep, in degrees (azimuth, altitude)
#         plot!(p[1], camera=(15*cos(i), 40))
#
#         # add a tracking line
#         fixed_x = zeros(40)
#         z = map(f, fixed_x, y)
#         plot!(p[1], fixed_x, y, z, line = (:black, 5, 0.2))
#         vline!(p[2], [0], line = (:black, 5))
#
#         # add to and show the tracked values over time
#         global zs = vcat(zs, z')
#         plot!(p[3], zs, alpha = 0.2, palette = cgrad(:blues).colors)
#
#         # increment the progress bar
#         next!(prog)
#         end
#     end)],
#     "iris" => [:(begin
#         gr()
#         # load a dataset
#         using RDatasets
#         iris = dataset("datasets", "iris");
#
#         # load the StatPlots recipes (for DataFrames) available via:
#         # Pkg.add("StatPlots")
#         using StatPlots
#         gr()
#
#         # Scatter plot with some custom settings
#         @df iris scatter(:SepalLength, :SepalWidth, group=:Species,
#             title = "My awesome plot",
#             xlabel = "Length", ylabel = "Width",
#             m=(0.5, [:cross :hex :star7], 12),
#             bg=RGB(.2,.2,.2))
#     end)],
#     "lines_1" => [:(begin
#         using Plots; gr()
#         x = 1:10; y = rand(10); # These are the plotting data
#         plot(x,y)
#     end)],
#     "lines_2" => [:(begin
#         x = 1:10; y = rand(10,2) # 2 columns means two lines
#         plot(x,y)
#     end)],
#     "lines_3" => [:(begin
#         z = rand(10)
#         plot!(x,z)
#     end)],
#     "lines_4" => [:(begin
#         x = 1:10; y = rand(10,2) # 2 columns means two lines
#         p = plot(x,y)
#         z = rand(10)
#         plot!(p,x,z)
#     end)],
#     "attr_1" => [:(begin
#         x = 1:10; y = rand(10,2) # 2 columns means two lines
#         plot(x,y,title="Two Lines",label=["Line 1" "Line 2"],lw=3)
#     end)],
#     "attr_2" => [:(begin
#         xlabel!("My x label")
#     end)],
#     "backends_1" => [:(begin
#         x = 1:10; y = rand(10,2) # 2 columns means two lines
#         plotlyjs() # Set the backend to Plotly
#         using ORCA
#         plot(x,y,title="This is Plotted using Plotly")
#     end)],
#     "backends_2" => [:(begin
#         gr() # Set the backend to GR
#         plot(x,y,title="This is Plotted using GR") # This plots using GR
#     end)],
#     "scatter_1" => [:(begin
#         gr() # We will continue onward using the GR backend
#         plot(x,y,seriestype=:scatter,title="My Scatter Plot")
#     end)],
#     "scatter_2" => [:(begin
#         gr()
#         scatter(x,y,title="My Scatter Plot")
#     end)],
#     "subplots_1" => [:(begin
#         y = rand(10,4)
#         plot(x,y,layout=(4,1))
#     end)],
#     "subplots_2" => [:(begin
#         p1 = plot(x,y) # Make a line plot
#         p2 = scatter(x,y) # Make a scatter plot
#         p3 = plot(x,y,xlabel="This one is labelled",lw=3,title="Subtitle")
#         p4 = histogram(x,y) # Four histograms each with 10 points? Why not!
#         plot(p1,p2,p3,p4,layout=(2,2),legend=false)
#     end)],
#     "user_recipes_1" => [:(begin
#         # Pkg.add("StatPlots")
#         using StatPlots # Required for the DataFrame user recipe
#         # Now let's create the DataFrame
#         using DataFrames
#         df = DataFrame(a = 1:10, b = 10*rand(10), c = 10 * rand(10))
#         # Plot the DataFrame by declaring the points by the column names
#         @df df plot(:a, [:b :c]) # x = :a, y = [:b :c]. Notice this is two columns!
#     end)],
#     "user_recipes_2" => [:(begin
#         @df df scatter(:a, :b, title="My DataFrame Scatter Plot!") # x = :a, y = :b
#     end)],
#     "type_recipes" => [:(begin
#         using Distributions
#         plot(Normal(3,5),lw=3)
#     end)],
#     "plot_recipes" => [:(begin
#         #Pkg.add("RDatasets")
#         using RDatasets
#         iris = dataset("datasets","iris")
#     end)],
#     "series_recipes_1" => [:(begin
#         y = rand(100,4) # Four series of 100 points each
#         violin(["Series 1" "Series 2" "Series 3" "Series 4"],y,leg=false)
#     end)],
#     "series_recipes_2" => [:(begin
#         boxplot!(["Series 1" "Series 2" "Series 3" "Series 4"],y,leg=false)
#     end)],
#     "columns_are_series" => [:(begin
#         using Plots; gr()
#
#         # 10 data points in 4 series
#         xs = 0 : 2π/10 : 2π
#         data = [sin.(xs) cos.(xs) 2sin.(xs) 2cos.(xs)]
#
#         # We put labels in a row vector: applies to each series
#         labels = ["Apples" "Oranges" "Hats" "Shoes"]
#
#         # Marker shapes in a column vector: applies to data points
#         markershapes = [:circle, :star5]
#
#         # Marker colors in a matrix: applies to series and data points
#         markercolors = [:green :orange :black :purple
#                     :red   :yellow :brown :white]
#
#         plot(xs, data, label = labels, shape = markershapes, color = markercolors,
#          markersize = 10)
#     end)],
#     "groups" => [:(begin
#         using Plots; plotlyjs()
#
#         function rectangle_from_coords(xb,yb,xt,yt)
#         [
#             xb yb
#             xt yb
#             xt yt
#             xb yt
#             xb yb
#             NaN NaN
#         ]
#         end
#
#         some_rects=[
#         rectangle_from_coords(1 ,1 ,5 ,5 )
#         rectangle_from_coords(10,10,15,15)
#         ]
#         other_rects=[
#         rectangle_from_coords(1 ,10,5 ,15)
#         rectangle_from_coords(10,1 ,15,5 )
#         ]
#
#         plot(some_rects[:,1], some_rects[:,2],label="some group")
#         plot!(other_rects[:,1], other_rects[:,2],label="other group")
#     end)],
#     "dataframes" => [:(begin
#         gr() # maybe remove again
#         using StatPlots, RDatasets
#         iris = dataset("datasets", "iris")
#         @df iris scatter(:SepalLength, :SepalWidth, group=:Species,
#             m=(0.5, [:+ :h :star7], 12), bg=RGB(.2,.2,.2))
#     end)],
#     "rand_anim" => [:(begin
#         using Plots, Iterators
#
#         itr = repeatedly(()->rand(10), 20)
#         anim = animate(itr, ylims=(0,1), c=:red, fps=5)
#     end)],
#     "my_count_1" => [:(begin
#         using Plots, StatsBase
#         gr(size=(300,300), leg=false)
#
#         @userplot MyCount
#         @recipe function f(mc::MyCount)
#             # get the array from the args field
#             arr = mc.args[1]
#
#             T = typeof(arr)
#             if T.parameters[1] == Float64
#                 seriestype := :histogram
#                 arr
#             else
#                 seriestype := :bar
#                 cm = countmap(arr)
#                 x = sort!(collect(keys(cm)))
#                 y = [cm[xi] for xi=x]
#                 x, y
#             end
#         end
#         mycount(rand(500))
#     end)],
#     "my_count_2" => [:(begin
#         mycount(rand(["A","B","C"],100))
#     end)],
#     "my_count_3" => [:(begin
#         mycount(rand(1:500, 500))
#     end)],
#     "pipeline_1" => [:(begin
#         using Plots; pyplot()
#         plot(x, y, line = (0.5, [4 1 0], [:path :scatter :density]),
#         marker=(10, 0.5, [:none :+ :none]), fill=0.5,
#         orientation = :h, title = "My title")
#     end)],
#     "pipeline_2" => [:(begin
#         gr()
#         mutable struct MyVecWrapper
#             v::Vector{Float64}
#         end
#         mv = MyVecWrapper(rand(100))
#
#         @recipe function f(mv::MyVecWrapper)
#             markershape --> :circle
#             markersize  --> 30
#             mv.v
#         end
#
#         plot(
#             plot(mv.v),
#             plot(mv)
#         )
#     end)],
#     "pipeline_3" => [:(begin
#         scatter(rand(100), group = rand(1:3, 100), marker = (10,0.3,[:s :o :x]))
#     end)],
#     "layout_1" => [:(begin
#         # create a 2x2 grid, and map each of the 4 series to one of the subplots
#         plot(rand(100,4), layout = 4)
#     end)],
#     "layout_2" => [:(begin
#         # create a 4x1 grid, and map each of the 4 series to one of the subplots
#         plot(rand(100,4), layout = (4,1))
#     end)],
#     "layout_3" => [:(begin
#         plot(rand(100,4), layout = grid(4,1,heights=[0.1,0.4,0.4,0.1]))
#     end)],
#     "layout_4" => [:(begin
#         l = @layout [  a{0.3w} [grid(3,3)
#                                          b{0.2h} ]]
#         plot(
#             rand(10,11),
#             layout = l, legend = false, seriestype = [:bar :scatter :path],
#             title = ["($i)" for j = 1:1, i=1:11], titleloc = :right, titlefont = font(8)
#         )
#     end)],
#     "layout_5" => [:(begin
#         # boxplot is defined in StatPlots
#         using StatPlots
#         gr(leg=false, bg=:lightgrey)
#
#         # Create a filled contour and boxplot side by side.
#         plot(contourf(randn(10,20)), boxplot(rand(1:4,1000),randn(1000)))
#
#         # Add a histogram inset on the heatmap.
#         # We set the (optional) position relative to bottom-right of the 1st subplot.
#         # The call is `bbox(x, y, width, height, origin...)`, where numbers are treated as "percent of parent"
#         histogram!(randn(1000), inset = (1, bbox(0.05,0.05,0.5,0.25,:bottom,:right)), ticks=nothing, subplot=3, bg_inside=nothing)
#
#         # Add sticks floating in the window (inset relative to the window, as opposed to being relative to a subplot)
#         sticks!(randn(100), inset = bbox(0,-0.2,200Plots.px,100Plots.px,:center), ticks=nothing, subplot=4)
#     end)],
# )
#
# _animation_imgs = ["lorenz_attractor", "waves", "rand_anim"]
#
# _doc_img_files = Dict(
#     "lorenz_attractor" => "index",
#     "waves" => "index",
#     "iris" => "index",
#     "lines_1" => "tutorial",
#     "lines_2" => "tutorial",
#     "lines_3" => "tutorial",
#     "lines_4" => "tutorial",
#     "attr_1" => "tutorial",
#     "attr_2" => "tutorial",
#     "backends_1" => "tutorial",
#     "backends_2" => "tutorial",
#     "scatter_1" => "tutorial",
#     "scatter_2" => "tutorial",
#     "subplots_1" => "tutorial",
#     "subplots_2" => "tutorial",
#     "user_recipes_1" => "tutorial",
#     "user_recipes_2" => "tutorial",
#     "type_recipes" => "tutorial",
#     "plot_recipes" => "tutorial",
#     "series_recipes_1" => "tutorial",
#     "series_recipes_2" => "tutorial",
#     "columns_are_series" => "input",
#     "groups" => "input",
#     "dataframes" => "input",
#     "rand_anim" => "animations",
#     "my_count_1" => "contributing",
#     "my_count_2" => "contributing",
#     "my_count_3" => "contributing",
#     "pipeline_1" => "pipeline",
#     "pipeline_2" => "pipeline",
#     "pipeline_3" => "pipeline",
#     "layout_1" => "layout",
#     "layout_2" => "layout",
#     "layout_3" => "layout",
#     "layout_4" => "layout",
#     "layout_5" => "layout",
# )
#
# _doc_img_names = [
#     "lorenz_attractor",
#     "waves",
#     "iris",
#     "lines_1",
#     "lines_2",
#     "lines_3",
#     "lines_4",
#     "attr_1",
#     "attr_2",
#     "backends_1",
#     "backends_2",
#     "scatter_1",
#     "scatter_2",
#     "subplots_1",
#     "subplots_2",
#     "user_recipes_1",
#     "user_recipes_2",
#     "type_recipes",
#     "plot_recipes",
#     "series_recipes_1",
#     "series_recipes_2",
#     "columns_are_series",
#     "groups",
#     "dataframes",
#     # "rand_anim", # requires Iterators to work
#     "my_count_1",
#     "my_count_2",
#     "my_count_3",
#     "pipeline_1",
#     "pipeline_2",
#     "pipeline_3",
#     "layout_1",
#     "layout_2",
#     "layout_3",
#     "layout_4",
#     "layout_5",
# ]
#
# function generate_doc_images()
#     for id in _doc_img_names
#         generate_doc_image(id)
#     end
# end
#
# function generate_doc_image(id::AbstractString)
#     map(ex -> Base.eval(Main, ex), _doc_images[id])
#
#     dir = local_path("PlotDocs", _doc_img_files[id])
#     isdir(dir) || mkpath(dir)
#
#     if id in _animation_imgs
#         fn = joinpath(dir, string(id, ".gif"))
#         anim = @eval Main anim
#         if typeof(anim) <: Plots.AnimatedGif
#             mv(anim.filename, fn, force = true)
#         else
#             Plots.gif(anim, fn, fps = 15)
#         end
#     else
#         fn = joinpath(dir, string(id, ".png"))
#         Plots.savefig(fn)
#     end
# end


text:
module ExtendedVariantDictionary

export IntegerVariants, FloatVariants, BooleanVariants

const IntegerVariants = [
	"AddSeekers",
	"BadelineBossCount",
	"BadelineBossNodeCount",
	"ChaserCount",
	"DashCount",
	"JellyfishEverywhere",
	"OshiroCount",
	"ReverseOshiroCount",
	"Stamina"
]

const FloatVariants = [
	"AirFriction",
	"AnxietyEffect",
	"BackgroundBlurLevel",
	"BackgroundBrightness",
	"BadelineLag",
	"BlurLevel",
	"BoostMultiplier",
	"CoyoteTime",
	"DashLength",
	"DashSpeed",
	"DelayBetweenBadelines",
	"ExplodeLaunchSpeed",
	"FallSpeed",
	"ForegroundEffectOpacity",
	"Friction",
	"GameSpeed",
	"GlitchEffect",
	"Gravity",
	"HiccupStrength",
	"HyperdashSpeed",
	"JumpHeight",
	"RegularHiccups",
	"RisingLavaSpeed",
	"RoomLighting",
	"RoomBloom",
	"ScreenShakeIntensity",
	"SnowballDelay",
	"SpeedX",
	"SuperdashSteeringSpeed",
	"SwimmingSpeed",
	"WallBouncingSpeed",
	"ZoomLevel"
]

const BooleanVariants = [
	"AffectExistingChasers",
	"AllStrawberriesAreGoldens",
	"AllowLeavingTheoBehind",
	"AllowThrowingTheoOffscreen",
	"AlwaysInvisible",
	"BadelineBossesEverywhere",
	"BadelineChasersEverywhere",
	"BounceEverywhere",
	"ChangePatternsOfExistingBosses",
	"DashTrailAllTheTime",
	"DisableClimbJumping",
	"DisableClimbingUpOrDown",
	"DisableMadelineSpotlight",
	"DisableNeutralJumping",
	"DisableOshiroSlowdown",
	"DisableRefillsOnScreenTransition",
	"DisableSeekerSlowdown",
	"DisableSuperBoosts",
	"DisableWallJumping",
	"DisplayDashCount",
	"DontRefillStaminaOnGround",
	"EveryJumpIsUltra",
	"EverythingIsUnderwater",
	"FirstBadelineSpawnRandom",
	"ForceDuckOnGround",
	"FriendlyBadelineFollower",
	"HeldDash",
	"InvertDashes",
	"InvertGrab",
	"InvertHorizontalControls",
	"InvertVerticalControls",
	"LegacyDashSpeedBehavior",
	"NoFreezeFrames",
	"OshiroEverywhere",
	"PreserveExtraDashesUnderwater",
	"RefillJumpsOnDashRefill",
	"RestoreDashesOnRespawn",
	"RisingLavaEverywhere",
	"SnowballsEverywhere",
	"TheoCrystalsEverywhere",
	"UpsideDown"
]

const SpecialHandlingVariants = [
	"BadelineAttackPattern", # technically integer, but not easy to use
	"DashDirection", # bool[][]
	"ColorGrading", # string
	"DisplaySpeedometer", # enum
	"DontRefillDashOnGround", # enum
	"JumpCount", # integer but with special "infinite" value
	"MadelineBackpackMode", # enum
	"WindEverywhere" # enum
]

const UnavailableVariants = [
	"MadelineIsSilhouette", # use Madeline Silhouette Trigger instead
	"MadelineHasPonytail" # use Madeline Ponytail Trigger instead
]

end


text:
function get_incompatible_devices(devices_template::Dict)
    incompatible_device_types = Set{DataType}()
    for model in values(devices_template)
        formulation = get_formulation(model)
        if formulation == FixedOutput
            if !isempty(get_services(model))
                @info "$(formulation) for $(get_component_type(model)) is not compatible with the provision of reserve services"
            end
            push!(incompatible_device_types, get_component_type(model))
        end
    end
    return incompatible_device_types
end

function map_contributing_devices_by_type!(
    service_model::ServiceModel,
    contributing_devices,
)
    types = unique(typeof.(contributing_devices))
    for S in types
        if isabstracttype(S)
            error("The device type $S can't be abstract")
        end
        _devices = [d for d in contributing_devices if isa(d, S)]
        add_contributing_devices_map!(service_model, S, _devices)
    end
    return
end

function populate_aggregated_service_model!(template, sys::PSY.System)
    services_template = get_service_models(template)
    for (key, service_model) in services_template
        attributes = get_attributes(service_model)
        if get(attributes, "aggregated_service_model", false)
            delete!(services_template, key)
            D = get_component_type(service_model)
            B = get_formulation(service_model)
            for service in get_available_components(D, sys)
                new_key = (PSY.get_name(service), Symbol(D))
                if !haskey(services_template, new_key)
                    set_service_model!(template, ServiceModel(D, B, PSY.get_name(service)))
                end
            end
        end
    end
    return
end

function populate_contributing_devices!(template, sys::PSY.System)
    service_models = get_service_models(template)
    isempty(service_models) && return

    device_models = get_device_models(template)
    incompatible_device_types = get_incompatible_devices(device_models)
    services_mapping = PSY.get_contributing_device_mapping(sys)
    for (service_key, service_model) in service_models
        S = get_component_type(service_model)
        service = PSY.get_component(S, sys, get_service_name(service_model))
        if service === nothing
            @info "The data doesn't include services of type $(S) and name $(get_service_name(service_model)), consider changing the service models" _group =
                LOG_GROUP_SERVICE_CONSTUCTORS
            continue
        end
        contributing_devices_ =
            services_mapping[(type=S, name=PSY.get_name(service))].contributing_devices
        contributing_devices = [
            d for d in contributing_devices_ if
            typeof(d) ∉ incompatible_device_types && PSY.get_available(d)
        ]
        if isempty(contributing_devices)
            @warn "The contributing devices for service $(PSY.get_name(service)) is empty, consider removing the service from the system" _group =
                LOG_GROUP_SERVICE_CONSTUCTORS
            continue
        end
        map_contributing_devices_by_type!(service_model, contributing_devices)
    end
    return
end

function construct_services!(
    container::OptimizationContainer,
    sys::PSY.System,
    stage::ArgumentConstructStage,
    services_template::ServicesModelContainer,
    devices_template::DevicesModelContainer,
)
    isempty(services_template) && return
    incompatible_device_types = get_incompatible_devices(devices_template)

    groupservice = nothing

    for (key, service_model) in services_template
        if get_formulation(service_model) === GroupReserve  # group service needs to be constructed last
            groupservice = key
            continue
        end
        isempty(get_contributing_devices(service_model)) && continue
        get_contributing_devices(service_model)
        construct_service!(
            container,
            sys,
            stage,
            service_model,
            devices_template,
            incompatible_device_types,
        )
    end
    groupservice === nothing || construct_service!(
        container,
        sys,
        stage,
        services_template[groupservice],
        devices_template,
        incompatible_device_types,
    )
    return
end

function construct_services!(
    container::OptimizationContainer,
    sys::PSY.System,
    stage::ModelConstructStage,
    services_template::ServicesModelContainer,
    devices_template::DevicesModelContainer,
)
    isempty(services_template) && return
    incompatible_device_types = get_incompatible_devices(devices_template)

    groupservice = nothing
    for (key, service_model) in services_template
        if get_formulation(service_model) === GroupReserve  # group service needs to be constructed last
            groupservice = key
            continue
        end
        isempty(get_contributing_devices(service_model)) && continue
        construct_service!(
            container,
            sys,
            stage,
            service_model,
            devices_template,
            incompatible_device_types,
        )
    end
    groupservice === nothing || construct_service!(
        container,
        sys,
        stage,
        services_template[groupservice],
        devices_template,
        incompatible_device_types,
    )
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::ServiceModel{SR, RangeReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.Reserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    add_parameters!(container, RequirementTimeSeriesParameter, service, model)
    contributing_devices = get_contributing_devices(model)

    add_variables!(
        container,
        ActivePowerReserveVariable,
        service,
        contributing_devices,
        RangeReserve(),
    )
    add_to_expression!(container, ActivePowerReserveVariable, model, devices_template)
    add_feedforward_arguments!(container, model, service)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::ServiceModel{SR, RangeReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.Reserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)

    add_constraints!(container, RequirementConstraint, service, contributing_devices, model)

    objective_function!(container, service, model)

    add_feedforward_constraints!(container, model, service)

    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::ServiceModel{SR, RangeReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.StaticReserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)

    add_variables!(
        container,
        ActivePowerReserveVariable,
        service,
        contributing_devices,
        RangeReserve(),
    )
    add_to_expression!(container, ActivePowerReserveVariable, model, devices_template)
    add_feedforward_arguments!(container, model, service)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::ServiceModel{SR, RangeReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.StaticReserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)

    add_constraints!(container, RequirementConstraint, service, contributing_devices, model)

    objective_function!(container, service, model)

    add_feedforward_constraints!(container, model, service)

    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::ServiceModel{SR, StepwiseCostReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.Reserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)
    add_variable!(container, ServiceRequirementVariable(), [service], StepwiseCostReserve())
    add_variables!(
        container,
        ActivePowerReserveVariable,
        service,
        contributing_devices,
        StepwiseCostReserve(),
    )
    add_to_expression!(container, ActivePowerReserveVariable, model, devices_template)
    add_expressions!(container, ProductionCostExpression, [service], model)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::ServiceModel{SR, StepwiseCostReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.Reserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)

    add_constraints!(container, RequirementConstraint, service, contributing_devices, model)

    objective_function!(container, service, model)

    add_feedforward_constraints!(container, model, service)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::ServiceModel{S, T},
    devices_template::Dict{Symbol, DeviceModel},
    ::Set{<:DataType},
) where {S <: PSY.AGC, T <: AbstractAGCFormulation}
    services = get_available_components(S, sys)
    agc_areas = PSY.get_area.(services)
    areas = PSY.get_components(PSY.Area, sys)
    if !isempty(setdiff(areas, agc_areas))
        throw(
            IS.ConflictingInputsError(
                "All area most have an AGC service assigned in order to model the System's Frequency regulation",
            ),
        )
    end

    add_variables!(container, SteadyStateFrequencyDeviation)
    add_variables!(container, AreaMismatchVariable, areas, T())
    add_variables!(container, SmoothACE, areas, T())
    add_variables!(container, LiftVariable, areas, T())
    add_variables!(container, ActivePowerVariable, areas, T())
    add_variables!(container, DeltaActivePowerUpVariable, areas, T())
    add_variables!(container, DeltaActivePowerDownVariable, areas, T())
    add_variables!(container, AdditionalDeltaActivePowerUpVariable, areas, T())
    add_variables!(container, AdditionalDeltaActivePowerDownVariable, areas, T())

    add_initial_condition!(container, services, T(), AreaControlError())

    add_to_expression!(
        container,
        EmergencyUp,
        AdditionalDeltaActivePowerUpVariable,
        areas,
        model,
    )

    add_to_expression!(
        container,
        EmergencyDown,
        AdditionalDeltaActivePowerDownVariable,
        areas,
        model,
    )

    add_to_expression!(container, RawACE, SteadyStateFrequencyDeviation, services, model)

    add_feedforward_arguments!(container, model, services)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::ServiceModel{S, T},
    devices_template::Dict{Symbol, DeviceModel},
    ::Set{<:DataType},
) where {S <: PSY.AGC, T <: AbstractAGCFormulation}
    areas = PSY.get_components(PSY.Area, sys)
    services = get_available_components(S, sys)

    add_constraints!(container, AbsoluteValueConstraint, LiftVariable, areas, model)
    add_constraints!(
        container,
        FrequencyResponseConstraint,
        SteadyStateFrequencyDeviation,
        areas,
        model,
        sys,
    )
    add_constraints!(
        container,
        SACEPIDAreaConstraint,
        SteadyStateFrequencyDeviation,
        areas,
        model,
        sys,
    )
    add_constraints!(container, BalanceAuxConstraint, SmoothACE, areas, model, sys)

    add_feedforward_constraints!(container, model, services)

    objective_function!(container, areas, model)
    return
end

"""
    Constructs a service for StaticReserveGroup.
"""
function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::ServiceModel{SR, GroupReserve},
    ::Dict{Symbol, DeviceModel},
    ::Set{<:DataType},
) where {SR <: PSY.StaticReserveGroup}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_services = PSY.get_contributing_services(service)
    # check if variables exist
    check_activeservice_variables(container, contributing_services)

    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::ServiceModel{SR, GroupReserve},
    ::Dict{Symbol, DeviceModel},
    ::Set{<:DataType},
) where {SR <: PSY.StaticReserveGroup}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_services = PSY.get_contributing_services(service)

    add_constraints!(
        container,
        RequirementConstraint,
        service,
        contributing_services,
        model,
    )

    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::ServiceModel{SR, RampReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.Reserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)
    add_parameters!(container, RequirementTimeSeriesParameter, service, model)

    add_variables!(
        container,
        ActivePowerReserveVariable,
        service,
        contributing_devices,
        RampReserve(),
    )
    add_to_expression!(container, ActivePowerReserveVariable, model, devices_template)
    add_feedforward_arguments!(container, model, service)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::ServiceModel{SR, RampReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.Reserve}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)

    add_constraints!(container, RequirementConstraint, service, contributing_devices, model)
    add_constraints!(container, RampConstraint, service, contributing_devices, model)

    objective_function!(container, service, model)

    add_feedforward_constraints!(container, model, service)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::ServiceModel{SR, NonSpinningReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.ReserveNonSpinning}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)
    add_parameters!(container, RequirementTimeSeriesParameter, service, model)

    add_variables!(
        container,
        ActivePowerReserveVariable,
        service,
        contributing_devices,
        NonSpinningReserve(),
    )
    add_feedforward_arguments!(container, model, service)
    return
end

function construct_service!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::ServiceModel{SR, NonSpinningReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Set{<:DataType},
) where {SR <: PSY.ReserveNonSpinning}
    name = get_service_name(model)
    service = PSY.get_component(SR, sys, name)
    contributing_devices = get_contributing_devices(model)

    add_constraints!(container, RequirementConstraint, service, contributing_devices, model)
    add_constraints!(
        container,
        ReservePowerConstraint,
        service,
        contributing_devices,
        model,
    )

    objective_function!(container, service, model)

    add_feedforward_constraints!(container, model, service)
    return
end


text:
using StatsBase
using LinearAlgebra

"""
	Ray

Represents a ray, characterised by a starting point (p), and a direction (d).
From these, a unit vector, u can be computed.

# Fields
- `p::Vector{Float64}`  : Starting point
- `d::Vector{Float64}`  : Direction vector
- `u::Vector{Float64}`  : Unit vector

"""
struct Ray
    p::Vector{Float64}
    d::Vector{Float64}
    u::Vector{Float64}

    function Ray(p , d)
        new(p, d, unit_vector(d,p))
    end
end

"""
    unit_vector(p0::Vector{Real}, p1::Vector{Real})

Return a unit vector in the direction defined by points p0 and p1
"""
unit_vector(p0::Vector{Real}, p1::Vector{Real}) = (p1 - p0) ./ mag(p1 - p0)

"""
	propagate_ray(r::Ray)

Return a function that propagates the ray "r" from "p", along "e".

# Fields
- `r::Ray`  : Ray to be propagated
"""
function propagate_ray(r::Ray)
    function xray(t::Real)
        return r.p + t * r.d
    end
    return xray
end

"""
    Represents a Cylinder
"""
struct Cylinder
    r   :: Float64
    zmin:: Float64
    zmax:: Float64
end

"""Normal to the cylinder barrel in point p

Uses equation of cylynder: F(x,y,z) = x^2 + y^2 - r^2 = 0
then n = Grad(F)_P /Norm(Grad(F))_P
n = (2x, 2y, 0)/sqrt(4x^2 + 4y^2) = (x,y,0)/r (P)
"""
normal_to_barrel(c::Cylinder, p:: Vector{Float64}) = [p[1], p[2], 0.0] ./c.r

"""
    cylinder_equation(c::Cylinder, p:: Vector{Real})

equation of cylynder: F(x,y,z) = x^2 + y^2 - r^2 = 0
"""
cylinder_equation(c::Cylinder, p:: Vector{Float64}) = p[1]^2 + p[2]^2 - c.r^2

length(c::Cylinder) = abs.(c.zmax - c.zmin)

perimeter(c::Cylinder) = 2 * π * c.r

area_barrel(c::Cylinder) = 2 * π * c.r * length(c)

area_endcap(c::Cylinder) = π * c.r^2

area(c::Cylinder) = area_barrel(c) + 2 * area_endcap(c)

volume(c::Cylinder)= area_endcap(c) * length(c)

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
using PolyChaos
using LinearAlgebra
using Distributions

## This is an example on the usage of PolyChaos for non-intrusive PCE using
## the projection and the regression approach and a comarison of them in the univariate setting.
## f = e^x, x ~ N(1, 0.1)

μ, σ = 1, 0.1
maxdegree = 7

# model equation
function model(x)
    exp.(μ + σ * x)
end

# Analytic refernce
k_values = 0:maxdegree
computeCoefficient_analytic(k::Int) = exp(μ + 0.5*σ^2) * σ^k / factorial(k)
y_ana = computeCoefficient_analytic.(k_values)

# Setup and compute PCE coefficients of x
op = GaussOrthoPoly(maxdegree, Nrec=2*maxdegree)
t2 = Tensor(2, op)

# ------- Projection -------
println("\t == Projection approach ==")

# Perform projection via numerical integration
function computeProjection(k::Int, model::Function)
    g(t) = model(t) * evaluate(k, t, op)
    γ = t2.get([k,k])
    integrate(g, op.quad) / γ
end

y_proj = computeProjection.(0:maxdegree, model)
println.(y_proj)
println()

# Comparison to analytical solution
println("Comapre coefficients analytic <-> projection:")
println(norm(y_ana - y_proj, Inf), "\n")

# ------- Regression -------
include("../src/regression.jl")

println("\t == Regression approach ==")

# FUTURE: Perform truncation somehow

# Draw n samples, where N > P has to hold
nSamples = maxdegree * 20
X = sampleMeasure(nSamples, op)

# Evaluate model -> vector Y
Y = model.(X)

# Build matrix Φ with ϕ(x(i))
# Φ = Array{Float64}(undef, nSamples, maxdegree+1)
Φ = [ evaluate(j, X[i], op) for i = 1:nSamples, j = 0:maxdegree]

# Ordinary least squares regression
y_reg = leastSquares(Φ, Y)
println.(y_reg)
println()

# Validation of PCE model
println("Comapre coefficients analytic <-> regression:")
println(norm(y_ana - y_reg, Inf))

# genError = empError(Y, Φ, y_reg)
# println("Determination coefficient R² (normalized empicial error): ", 1 - genError)

ϵLoo = looError(Y, Φ, y_reg)
println("Determination coefficient Q² (leave-one-out error): ", 1- ϵLoo)
println()

# ------- Monte-Carlo -------
# Evaluate model function on same set of previously drawn samples X
y_mc = model.(X)

# ------- Comparison of moments -------
println("Comparison of moments to analytic solution")

# Analytic moments for y
mean_ana = exp(μ + σ^2/2)
std_ana = sqrt(exp(2*μ + σ^2) * (exp(σ^2) - 1))
skew_ana = sqrt(exp(σ^2) - 1) * (exp(σ^2) + 2)

# MC moments
mean_mc, std_mc, skew_mc = mean(y_mc), std(y_mc), skewness(y_mc)
error_mean_mc = abs(mean_ana - mean_mc)
error_std_mc = abs(std_ana - std_mc)
println("\t\t\t error MC, mean: \t $(error_mean_mc)")
println("\t\t\t error MC, std: \t $(error_std_mc)")

# Projection moments
mean_proj = mean(y_proj, op)
std_proj = std(y_proj, op)
error_mean_proj = abs(mean_ana - mean_proj)
error_std_proj = abs(std_ana - std_proj)
println("\t\t\t error proj, mean: \t $(error_mean_proj)")
println("\t\t\t error proj, std: \t $(error_std_proj)")

# Regression moments
mean_reg = mean(y_reg, op)
std_reg = std(y_reg, op)
error_mean_reg = abs(mean_ana - mean_reg)
error_std_reg = abs(std_ana - std_reg)
println("\t\t\t error reg, mean: \t $(error_mean_reg)")
println("\t\t\t error reg, std: \t $(error_std_reg)")

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos


text:
# back propagation
export backprop

function backprop(net::Vector{Layer},inputBase::AFArray{Float32,3},outputBase::AFArray{Float32,3};nIter=100)
  err=Float64[]
  dataLayerIndex=find(layer->typeof(layer)==JuLearn.dataLayer,net)
  lossLayerIndex=find(layer->typeof(layer)==JuLearn.lossLayer,net)
  nData=size(inputBase)[end-1]
  for iter=1:nIter
    idx=(iter-1)*batchsize+1
    range=map(x->mod1(x,nData),idx:idx+batchsize-1)
    range=AFArray(range)
    for (i,index) in enumerate(dataLayerIndex)
      net[index].rₒ=inputBase[:,range,i]
      net[index].rᵢ=outputBase[:,range,i]
    end
    forward!(net)
    backward!(net)
    update!(net,iter)
    shareParams!(net)
    push!(err,sum(i->net[i].loss,lossLayerIndex))
    println("iteration=",iter)
  end
   err
end

function backprop(net::Vector{Layer},inputBase::AFArray{Float32,2},outputBase::AFArray{Float32,2};nIter=100)
  err=Float64[]
  dataLayerIndex=find(layer->typeof(layer)==JuLearn.dataLayer,net)
  lossLayerIndex=find(layer->typeof(layer)==JuLearn.lossLayer,net)
  nData=size(inputBase)[end]
  for iter=1:nIter
    idx=(iter-1)*batchsize+1
    range=map(x->mod1(x,nData),idx:idx+batchsize-1)
    range=AFArray(range)
    net[dataLayerIndex[1]].rₒ=inputBase[:,range]
    net[dataLayerIndex[1]].rᵢ=outputBase[:,range]

    forward!(net)
    backward!(net)
    update!(net,iter)
    shareParams!(net)
    push!(err,sum(i->net[i].loss,lossLayerIndex))
    println("iteration=",iter)
  end
   err
end


text:
export DeterministicPolicy, RandomPolicy, EpsilonGreedyPolicy, LinearPolicy, FunctionalPolicy

include("deterministic_policy.jl")
include("random_policy.jl")
include("epsilon_greedy_policy.jl")
include("linear_policy.jl")
include("functional_policy.jl")


text:
type LinearSolution
    names::Dict{Symbol, Array{Symbol, 1}}
    g_0::Array{Float64,1}
    g_y::Array{Float64,2}
    g_e::Array{Float64,2}
    sigma::Array{Float64,2}
end

function residuals(model, y)
    f = model.functions.f_static
    y = [model.calibration[p] for p in model.symbols[:endogenous]]
    e = [model.calibration[p] for p in model.symbols[:exogenous]]
    p = [model.calibration[p] for p in model.symbols[:parameters]]
    # check residuals are zero
    res = (f(Der{0}, [y; e],p))
    return res
end

function residuals(model)
    y = [model.calibration[p] for p in model.symbols[:endogenous]]
    return residuals(model, y)
end

function solve(model)

    symbols = model.symbols

    n_v = length(model.symbols[:endogenous])
    n_e = length(model.symbols[:exogenous])

    endogenous = symbols[:endogenous]
    exogenous = symbols[:exogenous]
    parameters = symbols[:parameters]

    it = Dolang.IncidenceTable(model.equations)

    fut_syms = it.by_date[1]
    cur_syms = it.by_date[0]
    past_syms = it.by_date[-1]

    fut_inds = [i for (i,v) in enumerate(endogenous) if v in it.by_date[1]]
    cur_inds = [i for (i,v) in enumerate(endogenous) if v in it.by_date[0]]
    past_inds = [i for (i,v) in enumerate(endogenous) if v in it.by_date[-1]]
    shocks_inds = [i for (i,v) in enumerate(exogenous) if v in it.by_date[0]]
    # parms_inds = [i for (i,v) in enumerate(parameters) if v in it.by_date[0]]
    parms_inds = 1:length(model.symbols[:parameters])

    y_ss = model.calibration_grouped[:endogenous]
    e_ss = model.calibration_grouped[:exogenous]
    p_ss = model.calibration_grouped[:parameters]

    v = cat(1, y_ss[fut_inds], y_ss[cur_inds], y_ss[past_inds], e_ss[shocks_inds])
    p = p_ss[parms_inds]

    f = model.functions.f_dynamic

    # check residuals are zero
    # res = (f(Der{0}, v,p))
    # assert(maximum(abs(res))<1e-8)

    jac = f(Der{1}, v,p)

    n_v = length(y_ss)
    n_e = length(e_ss)

    F = zeros(n_v, n_v)
    G = zeros(n_v, n_v)
    H = zeros(n_v, n_v)

    K = zeros(n_v, n_e)

    n1 = length(fut_inds)
    n2 = length(cur_inds)
    n3 = length(past_inds)

    F[:,fut_inds] = jac[:,1:n1]
    G[:,cur_inds] = jac[:,n1+1:n1+n2]
    H[:,past_inds] = jac[:,n1+n2+1:n1+n2+n3]
    K[:,shocks_inds] = jac[:,n1+n2+n3+1:end]

    g_y = solve_second_order(F,G,H)
    g_e = -(F*g_y+G)\K

    names = Dict(
        :endogenous => model.symbols[:endogenous],
        :exogenous => model.symbols[:exogenous],
    )

    g_0 = model.calibration_grouped[:endogenous]
    sigma = model.exogenous.Σ.mat

    return LinearSolution(names, g_0, g_y, g_e, sigma)

end

function solve_second_order(F,G,H)

    n_v = size(F,1)
    N = size(F,1)

    D = [ zeros(N,N) F ;
                eye(N) zeros(N,N)]

    E = [ -H           -G ;
          zeros(N, N) eye(N)] ;

    sf = schurfact(D, E)

    if false in .~((sf.alpha .== 0) .*  (sf.beta .== 0))
        error("Indeterminate eigenvalues.")
    end

    ev = sf.beta ./ abs.(sf.alpha)
    ev[abs.(sf.alpha).==0.0] = Inf
    sort!(ev)
    delta = ev[n_v+1] - ev[n_v]
    if delta == 0.0
        error("MOD solution indeterminate.")
    end
    cutoff = ev[n_v] + (ev[n_v+1] - ev[n_v])/2

    select = (cutoff*abs.(sf.alpha)) .> sf.beta

    ordschur!(sf, select)

    T = sf.S
    S = sf.T

    eigvals = diag(sf.S)./diag(sf.T)

    Z = sf.Z'
    Q = sf.Q

    Z_21 = Z[n_v+1:end, 1:n_v]
    Z_22 = Z[n_v+1:end, n_v+1:end]

    g_1 = - Z_22 \ Z_21
    g_1 = real(g_1)

    return g_1
end

function simulate(sol::LinearSolution, sigma::Matrix{Float64}; T=20)
    n_v = size(sol.g_0,1)
    n_e = size(sol.sigma, 1)
    y_vec = zeros(n_v, (T+1))
    e_vec = zeros(n_e, (T+1))
    rnrm = MultivariateNormal(zeros(n_e), sigma)
    A = sol.g_y
    B = sol.g_e
    y = sol.g_0
    y_0 = sol.g_0
    y_vec[:,1] = y
    e_vec[:,2:end] = rand(rnrm,T)
    for t=2:T+1
        y = y_0 + A*(y-y_0) + B*e_vec[:,t]
        y_vec[:,t] = y
    end
    res = [y_vec; e_vec]

    tvec = collect(1:size(res,2)) - 1
    res = [tvec';res]

    colnames = cat(1, [:t], sol.names[:endogenous], sol.names[:exogenous])
    columns = Dict(colnames[i] => res[i,:] for i=1:size(res,1))

    return DataFrame(columns)

end

function simulate(sol; T=20)
    return simulate(sol, sol.sigma)
end

function simulate(model::Model; T=20)
    sol = solve(model)
    return simulate(sol, T=T)
end

function response(dr::LinearSolution, e0::Vector{Float64}; T=40)
    ysim = zeros((T+1,length(dr.g_0)))
    esim = zeros((T+1,size(dr.sigma,1)))
    y0 = dr.g_0
    ysim[1,:] = y0
    esim[2,:] = e0
    for t=1:T
        yp = ysim[t]
        ysim[t+1,:] = y0 + dr.g_y*(ysim[t,:]-y0) + dr.g_e*esim[t+1,:]
    end
    res = cat(2,ysim, esim)
    headers = [dr.names[:endogenous]; dr.names[:exogenous]]
    arg = Dict([headers[i]=>res[:,i] for i=1:size(res,2)])
    return DataFrames.DataFrame(arg)
end

function response(dr::LinearSolution, var::Symbol, value::Float64; T=40)
    index = findfirst(dr.names[:exogenous], var)
    e0 = zeros(size(dr.sigma,1))
    e0[index] = value
    return response(dr, e0; T=T)
end

function response(dr::LinearSolution, var::Symbol; T=40)
    index = findfirst(dr.names[:exogenous], var)
    e0 = zeros(size(dr.sigma,1))
    e0[index] = sqrt( dr.sigma[index,index] )
    return response(dr, e0; T=T)
end


text:
using DataStructures: OrderedDict

export
    DisplayStack,
    render

struct DisplayStack
    elements::OrderedDict{Symbol, DisplayObj}
end

function DisplayStack()
    return DisplayStack(OrderedDict{Symbol, DisplayObj}())
end

function Base.push!(stack::DisplayStack, e::Pair)
    push!(stack.elements, e)
end

function Base.push!(stack::DisplayStack, e::DisplayObj)
    push!(stack.elements, gensym("component")=>e)
end

function Base.getindex(stack::DisplayStack, key::Symbol)
    return stack.elements[key]
end

function render(stream::IO, stack::DisplayStack)
    T.buffered(stream) do buffer
        foreach(e->render(buffer, e.second), stack.elements)
        T.cmove_line_last(buffer)
    end
end

render(stack::DisplayStack) = render(T.out_stream, stack)


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9]);F = [2704*x1^2*x3^2+12050*x1^2*x6^2+3040*x1^2*x6*x9+2657*x1^2*x9^2+5096*x1*x2*x3*x8+9408*x1*x2*x4*x9-10750*x1*x2*x6*x9-1792*x1*x2*x9^2+6968*x1*x3*x4^2+4180*x1*x3*x4*x6-19506*x1*x3*x6^2-3040*x1*x3*x6*x9+12034*x1*x3*x7*x9-9240*x1*x4*x5*x6-14440*x1*x4*x6*x7-2432*x1*x4*x7*x9+880*x1*x5*x6*x9+14630*x1*x6^3+2464*x1*x6^2*x9+9702*x1*x6*x7*x9+8330*x1*x7*x9^2+9216*x2^2*x4^2+2401*x2^2*x8^2+3137*x2^2*x9^2+10944*x2*x3*x4*x7-76*x2*x3*x4*x9+10640*x2*x3*x6*x9+6566*x2*x4^2*x8+168*x2*x4*x5*x9+19008*x2*x4*x6*x7+24832*x2*x4*x7*x9-16*x2*x5*x9^2-1372*x2*x6^2*x8-8624*x2*x6^2*x9+6076*x2*x7*x8*x9+1444*x3^2*x4^2+9025*x3^2*x6^2+3249*x3^2*x7^2-6384*x3*x4^2*x5+608*x3*x4*x5*x9+14440*x3*x4*x6*x7-14630*x3*x6^3+11286*x3*x6*x7^2+9690*x3*x7^2*x9+4489*x4^4+7056*x4^2*x5^2-1876*x4^2*x6^2+5776*x4^2*x7^2+8308*x4^2*x7*x9-1344*x4*x5^2*x9-11704*x4*x6^2*x7+64*x5^2*x9^2+6125*x6^4+9801*x6^2*x7^2-1736*x6^2*x7*x9+16830*x6*x7^2*x9+11069*x7^2*x9^2+10670*x1*x2*x6-2156*x1*x3*x9-194*x2^2*x9+3148*x2*x3*x4-16296*x2*x4*x5+1552*x2*x5*x9-2508*x3^2*x7-4356*x3*x6*x7-3740*x3*x7*x9-2080*x1*x3+9409*x2^2-1960*x2*x8+484*x3^2-2680*x4^2+560*x6^2-2480*x7*x9+400, 5096*x1*x3*x8+9408*x1*x4*x9-10750*x1*x6*x9-1792*x1*x9^2+18432*x2*x4^2+4802*x2*x8^2+6274*x2*x9^2+10944*x3*x4*x7-76*x3*x4*x9+10640*x3*x6*x9+6566*x4^2*x8+168*x4*x5*x9+19008*x4*x6*x7+24832*x4*x7*x9-16*x5*x9^2-1372*x6^2*x8-8624*x6^2*x9+6076*x7*x8*x9+10670*x1*x6-388*x2*x9+3148*x3*x4-16296*x4*x5+1552*x5*x9+18818*x2-1960*x8, 5408*x1^2*x3+5096*x1*x2*x8+6968*x1*x4^2+4180*x1*x4*x6-19506*x1*x6^2-3040*x1*x6*x9+12034*x1*x7*x9+10944*x2*x4*x7-76*x2*x4*x9+10640*x2*x6*x9+2888*x3*x4^2+18050*x3*x6^2+6498*x3*x7^2-6384*x4^2*x5+608*x4*x5*x9+14440*x4*x6*x7-14630*x6^3+11286*x6*x7^2+9690*x7^2*x9-2156*x1*x9+3148*x2*x4-5016*x3*x7-4356*x6*x7-3740*x7*x9-2080*x1+968*x3, 9408*x1*x2*x9+13936*x1*x3*x4+4180*x1*x3*x6-9240*x1*x5*x6-14440*x1*x6*x7-2432*x1*x7*x9+18432*x2^2*x4+10944*x2*x3*x7-76*x2*x3*x9+13132*x2*x4*x8+168*x2*x5*x9+19008*x2*x6*x7+24832*x2*x7*x9+2888*x3^2*x4-12768*x3*x4*x5+608*x3*x5*x9+14440*x3*x6*x7+17956*x4^3+14112*x4*x5^2-3752*x4*x6^2+11552*x4*x7^2+16616*x4*x7*x9-1344*x5^2*x9-11704*x6^2*x7+3148*x2*x3-16296*x2*x5-5360*x4, -9240*x1*x4*x6+880*x1*x6*x9+168*x2*x4*x9-16*x2*x9^2-6384*x3*x4^2+608*x3*x4*x9+14112*x4^2*x5-2688*x4*x5*x9+128*x5*x9^2-16296*x2*x4+1552*x2*x9, 24100*x1^2*x6+3040*x1^2*x9-10750*x1*x2*x9+4180*x1*x3*x4-39012*x1*x3*x6-3040*x1*x3*x9-9240*x1*x4*x5-14440*x1*x4*x7+880*x1*x5*x9+43890*x1*x6^2+4928*x1*x6*x9+9702*x1*x7*x9+10640*x2*x3*x9+19008*x2*x4*x7-2744*x2*x6*x8-17248*x2*x6*x9+18050*x3^2*x6+14440*x3*x4*x7-43890*x3*x6^2+11286*x3*x7^2-3752*x4^2*x6-23408*x4*x6*x7+24500*x6^3+19602*x6*x7^2-3472*x6*x7*x9+16830*x7^2*x9+10670*x1*x2-4356*x3*x7+1120*x6, 12034*x1*x3*x9-14440*x1*x4*x6-2432*x1*x4*x9+9702*x1*x6*x9+8330*x1*x9^2+10944*x2*x3*x4+19008*x2*x4*x6+24832*x2*x4*x9+6076*x2*x8*x9+6498*x3^2*x7+14440*x3*x4*x6+22572*x3*x6*x7+19380*x3*x7*x9+11552*x4^2*x7+8308*x4^2*x9-11704*x4*x6^2+19602*x6^2*x7-1736*x6^2*x9+33660*x6*x7*x9+22138*x7*x9^2-2508*x3^2-4356*x3*x6-3740*x3*x9-2480*x9, 5096*x1*x2*x3+4802*x2^2*x8+6566*x2*x4^2-1372*x2*x6^2+6076*x2*x7*x9-1960*x2, 3040*x1^2*x6+5314*x1^2*x9+9408*x1*x2*x4-10750*x1*x2*x6-3584*x1*x2*x9-3040*x1*x3*x6+12034*x1*x3*x7-2432*x1*x4*x7+880*x1*x5*x6+2464*x1*x6^2+9702*x1*x6*x7+16660*x1*x7*x9+6274*x2^2*x9-76*x2*x3*x4+10640*x2*x3*x6+168*x2*x4*x5+24832*x2*x4*x7-32*x2*x5*x9-8624*x2*x6^2+6076*x2*x7*x8+608*x3*x4*x5+9690*x3*x7^2+8308*x4^2*x7-1344*x4*x5^2+128*x5^2*x9-1736*x6^2*x7+16830*x6*x7^2+22138*x7^2*x9-2156*x1*x3-194*x2^2+1552*x2*x5-3740*x3*x7-2480*x7];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()