text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# ECOS.jl
# Contains the wrapper itself
#############################################################################

__precompile__()

module ECOS

import Compat: unsafe_string, unsafe_wrap

# Try to load the binary dependency
if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("ECOS not properly installed. Please run Pkg.build(\"ECOS\")")
end

# ver  [not exported]
# Returns the version of ECOS in use
function ver()
    ver_ptr = ccall((:ECOS_ver, ECOS.ecos), Ptr{UInt8}, ())
    return unsafe_string(ver_ptr)
end

macro ecos_ccall(func, args...)
    f = "ECOS_$(func)"
    quote
        @unix_only ret = ccall(($f,ECOS.ecos), $(args...))
        @windows_only ret = ccall(($f,ECOS.ecos), stdcall, $(args...))
        ret
    end
end

function __init__()
    vnum = VersionNumber(ver())
    if !(vnum.major == 2 && vnum.minor == 0)
        depsdir = realpath(joinpath(dirname(@__FILE__),"..","deps"))
        error("Current ECOS version installed is $(ver()), but we require version 2.0.*. On Linux and Windows, delete the contents of the `$depsdir` directory except for the files `build.jl` and `.gitignore`, then rerun Pkg.build(\"ECOS\"). On OS X, run `using Homebrew; Homebrew.update()` in Julia.")
    end
end

# @fieldtype
macro fieldtype(instance, field)
    return :(fieldtype(typeof($instance), $field))
end

include("ECOSSolverInterface.jl")  # MathProgBase interface
include("types.jl")  # All the types and constants defined in ecos.h

# setup  (direct interface)
# Provide ECOS with a problem in the form
# min c'x
# st  A*x = b
#     G*x in_K h,
#       or equivalently  h - G*x in K
#       or equivalently  G*x + s = h, s in K
# K is the product of R+ and multiple second-order cones, with the
# positive orthant first.
# Inputs:
#   n       Number of variables == length(x)
#   m       Number of inequality constraints == length(h)
#   p       Number of equality constraints == length(b)  (can be 0)
#   l       Dimension of the positive orthant,
#               i.e. the first l elements of s are >= 0
#   ncones  Number of second-order cones present in K
#   q       ncones-long vector of the length of the index sets of the SOCs
#               e.g. cone 1 is indices 4:6, cone 2 is indices 7:10
#                    ->  q = [3, 4]
#   e       Number of exponential cones present in problem
#   Gpr, Gjc, Gir
#           Non-zeros, column indices, and the row index arrays, respectively,
#           for the matrix G represented in column compressed storage (CCS) format
#           which is equivalent to Julia's SparseMatrixCSC
#   Apr, Ajc, Air
#           Equivalent to above for the matrix A.
#           Can be all nothing if no equalities are present.
#   c       Objective coefficients, length(c) == n
#   h       RHS for inequality constraints, length(h) == m
#   b       RHS for equality constraints, length(b) == b (can be nothing)
# Returns a pointer to the ECOS pwork structure (Cpwork in ECOS.jl). See
# types.jl for more information.
function setup(n::Int, m::Int, p::Int, l::Int, ncones::Int, q::Union{Vector{Int},Void}, e::Int,
                Gpr::Vector{Float64}, Gjc::Vector{Int}, Gir::Vector{Int},
                Apr::Union{Vector{Float64},Void}, Ajc::Union{Vector{Int},Void}, Air::Union{Vector{Int},Void},
                c::Vector{Float64}, h::Vector{Float64}, b::Union{Vector{Float64},Void}; kwargs...)
    # Convert to canonical forms
    q = (q == nothing) ? convert(Ptr{Clong}, C_NULL) : convert(Vector{Clong},q)
    Apr = (Apr == nothing) ? convert(Ptr{Cdouble}, C_NULL) : Apr
    Ajc = (Ajc == nothing) ? convert(Ptr{Cdouble}, C_NULL) : convert(Vector{Clong},Ajc)
    Air = (Air == nothing) ? convert(Ptr{Cdouble}, C_NULL) : convert(Vector{Clong},Air)
    b = (b == nothing) ? convert(Ptr{Cdouble}, C_NULL) : b
    problem_ptr = ccall((:ECOS_setup, ECOS.ecos), Ptr{Cpwork},
        (Clong, Clong, Clong, Clong, Clong, Ptr{Clong}, Clong,
         Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong},
         Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong},
         Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
        n, m, p, l, ncones, q, e,
        Gpr, convert(Vector{Clong},Gjc), convert(Vector{Clong},Gir),
        Apr, Ajc, Air,
        c, h, b)

    problem_ptr != C_NULL || error("ECOS failed to construct problem.")

    if !isempty(kwargs)
        problem = unsafe_load(problem_ptr)
        problem.stgs != C_NULL || error("ECOS returned a malformed settings struct.")
        settings = unsafe_load(problem.stgs)

        options = Dict{Symbol, Any}()
        for (k,v) in kwargs
            options[k] = v
        end

        new_settings = Csettings([
            setting in keys(options) ?
            convert(@fieldtype(settings, setting), options[setting]) :
            getfield(settings, setting)
            for setting in fieldnames(settings)]...)

        unsafe_store!(problem.stgs, new_settings)
    end

    problem_ptr
end

# setup  (more general interface)
# A more tolerant version of the above method that doesn't require
# user to fiddle with internals of the sparse matrix format
# User can pass nothing as argument for A, b, and q
function setup(n, m, p, l, ncones, q, e, G, A, c, h, b; options...)
    if A == nothing
        if b != nothing
            @assert length(b) == 0
            b = nothing
        end
        @assert p == 0
        Apr = nothing
        Ajc = nothing
        Air = nothing
    else
        numrow, numcol = size(A)
        @assert numcol == n
        @assert numrow == length(b)
        sparseA = sparse(A)
        Apr = convert(Vector{Float64}, sparseA.nzval)
        Ajc = sparseA.colptr - 1  # C is 0-based
        Air = sparseA.rowval - 1
    end

    sparseG = sparse(G)
    Gpr = convert(Vector{Float64}, sparseG.nzval)
    Gjc = sparseG.colptr - 1  # C is 0-based
    Gir = sparseG.rowval - 1

    setup(  n, m, p, l, ncones, q, e,
            Gpr, Gjc, Gir,
            Apr, Ajc, Air,
            c, h, b; options...)
end

# solve
# Solves the provided problem. Results are stored inside the structure,
# but currently there is no convenient interface-provided way to access
# this - use MathProgBase interface.
function solve(problem::Ptr{Cpwork})
    exitflag = ccall((:ECOS_solve, ECOS.ecos), Clong, (Ptr{Cpwork},), problem)
end

# cleanup
# Frees memory allocated by ECOS for the problem.
# The optional keepvars argument is number of variables to NOT free.
function cleanup(problem::Ptr{Cpwork}, keepvars::Int = 0)
    ccall((:ECOS_cleanup, ECOS.ecos), Void, (Ptr{Cpwork}, Clong), problem, keepvars)
end

end # module


text:
module Demo

include("function_1.jl")

export my_f

end

ext: jl


text:
using SafeTestsets

@safetestset "PythonCallChainRules.jl" begin
    @safetestset "torch" begin
        include("test_pytorch.jl")
    end
    @safetestset "jax" begin
        include("test_jax.jl")
    end
end


text:
@defcomp damages begin
    DAMAGES = Variable(index=[time])    # Damages (trillions 2005 USD per year)
    DAMFRAC = Variable(index=[time])    # Damages (fraction of gross output)

    TATM    = Parameter(index=[time])   # Increase temperature of atmosphere (degrees C from 1900)
    YGROSS  = Parameter(index=[time])   # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    a1      = Parameter()               # Damage coefficient
    a2      = Parameter()               # Damage quadratic term
    a3      = Parameter()               # Damage exponent
    damadj  = Parameter()               # Adjustment exponent in damage function
    usedamadj = Parameter{Bool}()       # Only the Excel version uses the damadj parameter

    function run_timestep(p, v, d, t)
        # Define function for DAMFRAC
        v.DAMFRAC[t] = p.a1 * p.TATM[t] + p.a2 * p.TATM[t]^p.a3

        # Define function for DAMAGES
        if p.usedamadj
            # Excel version
            v.DAMAGES[t] = p.YGROSS[t] * v.DAMFRAC[t] / (1 + v.DAMFRAC[t]^p.damadj)
        else
            # GAMS Version
            v.DAMAGES[t] = p.YGROSS[t] * v.DAMFRAC[t]
        end
    end
end


text:
# https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities
using JuMP
using DisjunctiveProgramming

m = Model()
@variable(m, 0<=x[1:2]<=10)

@constraint(m, con1[i=1:2], x[i]<=[3,4][i])
@constraint(m, con2, 0<=x[1])
@constraint(m, con3, 0<=x[2])
@constraint(m, con4[i=1:2], [5,4][i]<=x[i])
@constraint(m, con5, x[1] <= 9)
@constraint(m, con6, x[2] <= 6)

@disjunction(m,(con1,con2,con3),(con4,con5,con6),reformulation=:BMR)


text:
struct RobbinsMonro{T}
    κ::T
    τ::T
end

function RobbinsMonro(κ::Real=0.51, τ::Real=1)
    0.5 < κ <= 1 || error("κ should be in the interval (0.5,1]")
    τ > 0 || error("τ should be positive")
    return RobbinsMonro{promote_type(typeof(κ), typeof(τ))}(κ, τ)
end

Optimisers.init(::RobbinsMonro, ::Any) = 1

function Optimisers.apply(o::RobbinsMonro, st, x, Δ)
    κ = o.κ
    τ = o.τ
    n = st
    return (n + 1), Δ * 1 / (τ + n)^κ
end

"""
    ALRSVI(n_mc=10, \rho)

Adaptive Learning Rate for Stochastic Variational Inference
"""
struct ALRSVI{T}
    n_mc::Int # Number of initial MC steps
    ρ::T # Initial learning rate
end

function Optimisers.init(opt::ALRSVI{T}, x::AbstractArray) where {T}
    i = 1
    g = zero(x)
    h = norm(g)
    τ = zero(T)
    return (; i, g, h, ρ=opt.ρ, τ)
end

function apply(opt::ALRSVI, state, x::AbstractArray, Δx::AbstractArray)
    if state.i <= opt.n_mc
        g = state.g + Δx
        h = state.h + norm(Δx)
        ρ = opt.ρ
        τ = state.τ
        if state.i == opt.τ
            g = g / τ
            h = h / τ
            ρ = sum(abs2, g) / h
        end
    else
        g = (1 - 1 / opt.τ) * state.g + 1 / opt.τ * Δx
        h = (1 - 1 / opt.τ) * state.h + 1 / opt.τ * sum(abs2, Δx)
        ρ = sum(abs2, g) / h
        τ = state.τ * (1 - ρ) + 1.0
    end

    return (; i=(i + 1), g, h, ρ, τ), ρ * Δx
end


text:
using Test
using TestSetExtensions
using IntensityMetrics
using IntensityMetrics.Unitful

@testset "perfect sinusoid" begin
    # ultrasound parameters
    f0 = 5e6u"Hz"                           # center frequency
    pressure_amplitude_1 = 1u"MPa"
    pressure_amplitude_2 = 2u"MPa"
    pressure_amplitude_3 = 3u"MPa"
    duty_cycle = 1.0
    pulse_duration = 10u"µs"
    total_duration = pulse_duration
    fs = 40e6u"Hz"                          # sampling frequency of received signal
    medium = Medium();
    excitation = Excitation(f0, pulse_duration, duty_cycle, total_duration, fs);

    # waveform parameters
    dt = 1/fs
    num_samples = ceil(pulse_duration/dt)
    time = collect((0:num_samples-1)*dt);
    data1 = sin.(2*pi*f0*time) * pressure_amplitude_1;
    data2 = cos.(2*pi*f0*time) * pressure_amplitude_2;
    data3 = cos.(2*pi*f0*time) * pressure_amplitude_3;
    data_all = repeat(hcat(data1, data2, data3), 1, 1, 2);

    @testset "intensity tests" begin
        intensity_1_expected = uconvert(u"W/m^2", pressure_amplitude_1^2 / (medium.density * medium.c))

        intensity_1_calculated = intensity(pressure_amplitude_1, medium)
        @test intensity_1_calculated ≈ intensity_1_expected atol=0.001u"W/m^2"

        intensity_1_calculated = intensity(uconvert(u"Pa", pressure_amplitude_1), medium)
        @test intensity_1_calculated ≈ intensity_1_expected atol=0.001u"W/m^2"
    end

    @testset "intensity_sppa tests" begin
        intensity_sppa_1_expected = uconvert(u"W/cm^2", pressure_amplitude_1^2 / (2 * medium.density * medium.c))
        intensity_sppa_2_expected = uconvert(u"W/cm^2", pressure_amplitude_2^2 / (2 * medium.density * medium.c))
        intensity_sppa_3_expected = uconvert(u"W/cm^2", pressure_amplitude_3^2 / (2 * medium.density * medium.c))
        intensity_sppa_expected = repeat(vcat(intensity_sppa_1_expected, intensity_sppa_2_expected, intensity_sppa_3_expected), 1, 2)

        # output is expected to be a vector
        intensity_sppa_1_calculated = intensity_sppa(data1, medium, excitation)
        @test intensity_sppa_1_calculated ≈ [intensity_sppa_1_expected] atol=0.001u"W/cm^2"

        intensity_sppa_calculated = intensity_sppa(data_all, medium, excitation)
        @test intensity_sppa_calculated ≈ intensity_sppa_expected atol=0.001u"W/cm^2"
    end

    @testset "mechanical_index tests" begin
        mechanical_index_1_expected = ustrip(uconvert(u"MPa", pressure_amplitude_1) / sqrt(uconvert(u"MHz", f0)))
        mechanical_index_2_expected = ustrip(uconvert(u"MPa", pressure_amplitude_2) / sqrt(uconvert(u"MHz", f0)))
        mechanical_index_3_expected = ustrip(uconvert(u"MPa", pressure_amplitude_3) / sqrt(uconvert(u"MHz", f0)))
        mechanical_index_expected = repeat(vcat(mechanical_index_1_expected, mechanical_index_2_expected, mechanical_index_3_expected), 1, 2)

        mechanical_index_1_calculated = mechanical_index(data1, excitation)
        @test mechanical_index_1_calculated ≈ [mechanical_index_1_expected] atol=0.001

        mechanical_index_1_calculated = mechanical_index(uconvert.(u"Pa", data1), excitation)
        @test mechanical_index_1_calculated ≈ [mechanical_index_1_expected] atol=0.001

        mechanical_index_calculated = mechanical_index(uconvert.(u"Pa", data_all), excitation)
        @test mechanical_index_calculated ≈ mechanical_index_expected atol=0.001
    end

    @testset "peak pos pressure tests" begin
        ppp_1_expected = pressure_amplitude_1
        ppp_2_expected = pressure_amplitude_2
        ppp_3_expected = pressure_amplitude_3
        ppp_expected = repeat(vcat(ppp_1_expected, ppp_2_expected, ppp_3_expected), 1, 2)

        @test pp_pressure(data1) ≈ [ppp_1_expected] atol=1u"Pa"
        @test pp_pressure(data_all) ≈ ppp_expected atol=1u"Pa"
    end

    @testset "peak neg pressure tests" begin
        pnp_1_expected = abs(-pressure_amplitude_1)
        pnp_2_expected = abs(-pressure_amplitude_2)
        pnp_3_expected = abs(-pressure_amplitude_3)
        pnp_expected = repeat(vcat(pnp_1_expected, pnp_2_expected, pnp_3_expected), 1, 2)

        @test pn_pressure(data1) ≈ [pnp_1_expected] atol=1u"Pa"
        @test pn_pressure(data_all) ≈ pnp_expected atol=1u"Pa"
    end

    @testset "ptp pressure tests" begin
        ptp_pressure_1_expected = 2 * pressure_amplitude_1
        ptp_pressure_2_expected = 2 * pressure_amplitude_2
        ptp_pressure_3_expected = 2 * pressure_amplitude_3
        ptp_pressure_expected = repeat(vcat(ptp_pressure_1_expected, ptp_pressure_2_expected, ptp_pressure_3_expected), 1, 2)

        @test ptp_pressure(data1) ≈ [ptp_pressure_1_expected] atol=1u"Pa"
        @test ptp_pressure(data_all) ≈ ptp_pressure_expected atol=1u"Pa"
    end
end


text:
using ADCME
using PyPlot
matplotlib.use("agg")

n = 80
pml = 30
N = n + 2pml + 1
NT = 1500
Δt = 1.5/NT
Δx = 1/n
x0 = LinRange(0, 1, n+1)
xE = Array((0:n+2pml)*Δx .- pml*Δx)
xH = (xE[2:end]+xE[1:end-1])/2
C = 1000.0

## PML
σE = zeros(N)
for i = 1:pml
    d = i*Δx
    σE[pml + n + 1 + i] = C* (d/(pml*Δx))^3
    σE[pml+1-i] = C* (d/(pml*Δx))^3
end

σH = zeros(N-1)
for i = 1:pml
    d = (i-1/2)*Δx
    σH[pml + n + i] = C* (d/(pml*Δx))^3
    σH[pml+1-i] = C* (d/(pml*Δx))^3
end

## Source
function ricker(dt = 0.002, f0 = 3.0)
    nw = 2/(f0*dt)
    nc = floor(Int, nw/2)
    t = dt*collect(-nc:1:nc)
    b = (π*f0*t).^2
    w = @. (1 - 2b)*exp(-b)
    # w = @. exp(-b)
end
R = ricker() * 1e7
if length(R)<NT+1
    R = [R;zeros(NT+1-length(R))]
else
    R = R[1:NT+1]
end
R_ = constant(R)

idx_src = pml + Int(round(0.0/Δx))+1
Z = zeros(N)
Z[idx_src] = 1.0
Z = Z[2:end-1]

idx_rcv = [pml + Int(round(0.0/Δx))+2, pml + Int(round(1.0/Δx))+1]

## Anomaly
model = "step"
cH = ones(length(xH)) * 3.0
cH_init = copy(cH)
cH[pml + Int(round(0.6/Δx)):end] .= 1.0
# model = "slop"
# cH = ones(length(xH)) * 2.0
# cH_init = copy(cH)
# cH[pml + Int(round(0.3/Δx)):end] = range(2.0, 3.0, length=length(cH[pml + Int(round(0.3/Δx)):end]))

cH_true = copy(cH)
cE = (cH[1:end-1]+cH[2:end])/2

## Similation
function condition(i, E_arr, H_arr)
    i<=NT+1
end

function body(i, E_arr, H_arr)
    E = read(E_arr, i-1)
    H = read(H_arr, i-1)

    # second order
    # ΔH = cH * (E[2:end]-E[1:end-1])/Δx - σH*H
    # H += ΔH * Δt
    # ΔE = cE * (H[2:end]-H[1:end-1])/Δx - σE[2:end-1]*E[2:end-1] + R_[i] * Z
    # E = scatter_add(E, 2:N-1, ΔE * Δt)

    ## forth order
    ΔH = cH * (E[2:end]-E[1:end-1])/Δx - σH*H
    H += ΔH * Δt
    ΔH = 1/(24Δx) * cH[2:N-2] * (-E[4:end] + 3E[3:end-1] - 3E[2:end-2] + E[1:end-3])
    H = scatter_add(H, 2:N-2, ΔH * Δt)
    ΔE = cE * (H[2:end]-H[1:end-1])/Δx - σE[2:end-1]*E[2:end-1] + R_[i] * Z
    E = scatter_add(E, 2:N-1, ΔE * Δt)
    ΔE = 1/(24Δx) * cE[2:end-1] * (-H[4:end] + 3H[3:end-1] - 3H[2:end-2] + H[1:end-3])
    E = scatter_add(E, 3:N-2, ΔE*Δt)

    i+1, write(E_arr, i, E), write(H_arr, i, H)
end

E_arr = TensorArray(NT+1)
H_arr = TensorArray(NT+1)

E_arr = write(E_arr, 1, zeros(N))
H_arr = write(H_arr, 1, zeros(N-1))

i = constant(2, dtype = Int32)

_, E, H = while_loop(condition, body, [i, E_arr, H_arr])

E = stack(E)
H = stack(H)

sess = Session(); init(sess)
E_true, H_true = run(sess, [E, H])

## data
figure(figsize = (14,4))
subplot(132)
pcolormesh(xE, (0:NT)*Δt, E_true, cmap="seismic", vmax=maximum(E_true), vmin=-maximum(E_true), rasterized=true)
plot(xE[idx_src]*ones(NT+1), (0:NT)*Δt, "r", label="Source", alpha=0.8)
plot(xE[idx_rcv[1]]*ones(NT+1), (0:NT)*Δt, "--", label="Receiver #1", alpha=0.8)
plot(xE[idx_rcv[2]]*ones(NT+1), (0:NT)*Δt, "--", label="Receiver #2", alpha=0.8)
autoscale(enable=true, axis="x", tight=true)
legend(loc="upper center")
xlabel("Localtion (km)")
ylabel("Time (s)")
subplot(133)
plot((0:NT)*Δt, E_true[:,idx_rcv[1]], label="Receiver #1", alpha=0.8)
plot((0:NT)*Δt, E_true[:,idx_rcv[2]], label="Receiver #2", alpha=0.8)
legend()
gca().ticklabel_format(axis="y", style="sci", scilimits=(0,0))
autoscale(enable=true, axis="x", tight=true)
xlabel("Time (s)")
ylabel("Amplitude")
subplot(131)
plot(xH, cH_true)
autoscale(enable=true, axis="x", tight=true)
xlabel("Location (km)")
ylabel("Velocity (km/s)")
tight_layout()
savefig("data.png", bbox_inch="tight")
savefig("data.pdf", bbox_inch="tight", dpi=300)


text:
type PsPot_HGH
  itype::Int
  atsymb::ASCIIString
  zval::Float64
  lloc::Int
  lmax::Int
  rloc::Float64
  rc::Array{Float64}   # indexed (l+1), l=0,1,2,3
  c::Array{Float64}    # indexed 1,2,3,4
  h::Array{Float64,3}  # originally indexed [0:3,1:3,1:3]
  k::Array{Float64,3}  # indexed [0:3,1:3,1:3]
  nprj::Array{Int}
  snprj::Int           # snprj = sum(nprj)
  lll::Array{Int}
  ipr::Array{Int}
end

# Constructor
function PsPot_HGH( itype::Int, atsymb::ASCIIString, filename::ASCIIString )

  @printf("\nFile = %s\n", filename)
  file = open( filename )

  comment = readline(file)

  lines = split( readline(file) )
  zval = parse( Float64, lines[1] )
  zion = parse( Float64, lines[2] )

  lines  = split( readline(file) )
  pspcod = parse( Int, lines[1] )
  pspxc  = parse( Int, lines[2] )
  lmax = parse( Int, lines[3] )
  lloc = parse( Int, lines[4] )
  mmax   = parse( Int, lines[5] )
  r2well = parse( Int, lines[6] )

  rc = zeros(Float64,4)
  c = zeros(Float64,4)
  h = zeros(Float64,4,3,3)
  k = zeros(Float64,4,3,3)

  lines = split( readline(file) )
  rloc = parse( Float64, lines[1] )
  c[1] = parse( Float64, lines[2] )
  c[2] = parse( Float64, lines[3] )
  c[3] = parse( Float64, lines[4] )
  c[4] = parse( Float64, lines[5] )

  const ANGMOM = ["s", "p", "d", "f"]

  l = 0  # s
  lines = split( readline(file) )
  #
  rc[1] = parse( lines[1] )
  #
  h[1,1,1] = parse( lines[2] )
  h[1,2,2] = parse( lines[3] )
  h[1,3,3] = parse( lines[4] )

  for l = 1:3
    lines = split( readline(file) )
    #
    rc[l+1] = parse( lines[1] )
    #
    h[l+1,1,1] = parse( lines[2] )
    h[l+1,2,2] = parse( lines[3] )
    h[l+1,3,3] = parse( lines[4] )
    #
    lines = split( readline(file) )
    k[l+1,1,1] = parse( lines[1] )
    k[l+1,2,2] = parse( lines[2] )
    k[l+1,3,3] = parse( lines[3] )
    #
  end

  # additional info
  lines = split(readline(file))
  nprj = zeros(Int,4)
  nprj[1] = parse( Int, lines[1] )
  nprj[2] = parse( Int, lines[2] )
  nprj[3] = parse( Int, lines[3] )
  nprj[4] = parse( Int, lines[4] )

  snprj  = sum(nprj)
  lll = zeros(Int,snprj)
  ipr = zeros(Int,snprj)

  iv = 0
  for il = 0:lmax
    for ii = 1:nprj[il+1]
      iv = iv + 1
      lll[iv] = il
      ipr[iv] = ii
    end
  end

  close(file)

  const M_HALF = 0.5
  const M_THREE = 3.0
  const M_FIVE = 5.0
  const M_ONE = 1.0

  # from Octopus code, see also the appendix of HGH paper

  h[0+1, 1, 2] = -M_HALF    * sqrt(M_THREE/M_FIVE) * h[0+1, 2, 2]
  h[0+1, 1, 3] =  M_HALF    * sqrt(M_FIVE/21.0)    * h[0+1, 3, 3]
  h[0+1, 2, 3] = -M_HALF    * sqrt(100.0/63.0)     * h[0+1, 3, 3]
  h[1+1, 1, 2] = -M_HALF    * sqrt(M_FIVE/7.0)     * h[1+1, 2, 2]
  h[1+1, 1, 3] =  M_ONE/6.0 * sqrt(35.0/11.0)      * h[1+1, 3, 3]
  h[1+1, 2, 3] = -M_ONE/6.0 * ( 14.0 / sqrt(11.0)) * h[1+1, 3, 3]
  h[2+1, 1, 2] = -M_HALF    * sqrt(7.0/9.0)        * h[2+1, 2, 2]
  h[2+1, 1, 3] =  M_HALF    * sqrt(63.0/143.0)     * h[2+1, 3, 3]
  h[2+1, 2, 3] = -M_HALF    * (18.0/sqrt(143.0))   * h[2+1, 3, 3]

  k[0+1, 1, 2] = -M_HALF    * sqrt(M_THREE/M_FIVE) * k[0+1, 2, 2]
  k[0+1, 1, 3] =  M_HALF    * sqrt(M_FIVE/21.0)    * k[0+1, 3, 3]
  k[0+1, 2, 3] = -M_HALF    * sqrt(100.0/63.0)     * k[0+1, 3, 3]
  k[1+1, 1, 2] = -M_HALF    * sqrt(M_FIVE/7.0)     * k[1+1, 2, 2]
  k[1+1, 1, 3] =  M_ONE/6.0 * sqrt(35.0/11.0)      * k[1+1, 3, 3]
  k[1+1, 2, 3] = -M_ONE/6.0 * (14.0 / sqrt(11.0))  * k[1+1, 3, 3]
  k[2+1, 1, 2] = -M_HALF    * sqrt(7.0/9.0)        * k[2+1, 2, 2]
  k[2+1, 1, 3] =  M_HALF    * sqrt(63.0/143.0)     * k[2+1, 3, 3]
  k[2+1, 2, 3] = -M_HALF    * (18.0/sqrt(143.0))   * k[2+1, 3, 3]

  # Parameters are symmetric.
  # be careful of index !!!
  for ik = 1:4
    for i = 1:3
      for j = i+1:3
        h[ik, j, i] = h[ik, i, j]
        k[ik, j, i] = k[ik, i, j]
      end
    end
  end

  psp = PsPot_HGH( itype, atsymb, zval, lloc, lmax, rloc,
                   rc, c, h, k, nprj, snprj, lll, ipr )

  return psp
end

function info_PsPot_HGH( psp::PsPot_HGH )

  const ANGMOM = ["s", "p", "d", "f"]

  @printf("Local pseudopotential info:\n")
  @printf("rloc: %f, c: %f, %f, %f, %f\n", psp.rloc, psp.c[1], psp.c[2], psp.c[3], psp.c[4])

  @printf("\nNonlocal pseudopotential info:\n")
  for i=1:4
    @printf("Angular momentum: %s, rc = %f\n", ANGMOM[i], psp.rc[i])
    @printf("h = \n")
    PrintMatrix( reshape(psp.h[i,:,:],(3,3) ) )
  end

end

# Evaluate HGH projector function in G-space
function eval_HGH_Vloc_G( psp, G2, Ω )

  Ng = size(G2)[1]
  Vg = zeros(Ng)

  rloc = psp.rloc
  zval = psp.zval
  c1 = psp.c[1]
  c2 = psp.c[2]
  c3 = psp.c[3]
  c4 = psp.c[4]

  pre1 = -4*pi*zval/Ω
  pre2 = sqrt(8*pi^3)*rloc^3/Ω
  #
  for ig=2:Ng
    Gr = sqrt(G2[ig])*rloc
    expGr2 = exp(-0.5*Gr^2)
    Vg[ig] = pre1/G2[ig]*expGr2 + pre2*expGr2 * (c1 + c2*(3-Gr^2) +
             c3*(15 - 10*Gr^2 + Gr^4) + c4*(105 - 105*Gr^2 + 21*Gr^4 - Gr^6) )
  end
  # limiting value, with minus sign ?
  Vg[1] = 2*pi*zval*rloc^2 + (2*pi)^1.5 * rloc^3 * (c1 + 3.0*c2 + 15*c3 + 105*c4)

  return Vg
end

# Evaluate HGH projector function in G-space
function eval_HGH_proj_G( psp, l, iproj, Gm, Ω )

  # NOTICE that G is the magnitudes of G-vectors
  Ng = size(Gm)[1]

  Vprj = zeros(Ng)

  rrl = psp.rc[l+1]

  # s-channel
  if l == 0

    if iproj==1

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = exp( -0.5*Gr2 )
      end

    elseif iproj==2

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = 2.0/sqrt(15.0) * exp( -0.5*Gr2 ) * ( 3.0 - Gr2 )
      end

    elseif iproj==3

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = (4.0/3.0)/sqrt(105.0) * exp( -0.5*Gr2 ) * (15.0 - 10.*Gr2 + Gr2^2)
       end

    end  # if iproj

  # p-channel
  elseif l == 1

    if iproj == 1

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = (1.0/sqrt(3.0)) * exp(-0.5*Gr2) * Gm[ig]
      end

    elseif iproj == 2

      for ig = 1:Ng
        Gr2 = (Gm[ig]*rrl)^2
        Vprj[ig] = (2./sqrt(105.)) * exp(-0.5*Gr2) * Gm[ig]*(5. - Gr2)
      end

    elseif iproj == 3

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl)^2
        Vprj[ig] = (4./3.)/sqrt(1155.) * exp(-0.5*Gr2) * Gm[ig] * (35. - 14.*Gr2 + Gr2^2)
      end

    end # if iproj

  # d-channel
  elseif l == 2

    if iproj == 1

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = (1.0/sqrt(15.0)) * exp(-0.5*Gr2) * Gm[ig]^2
      end

    elseif iproj == 2

      for ig = 1:Ng
        Gr2 = (Gm[ig]*rrl)^2
        Vprj[ig] = (2./3.)/sqrt(105.) * exp(-0.5*Gr2) * Gm[ig]^2 * (7.-Gr2)
      end

    end # if iproj

  # f-channel
  elseif l == 3

    for ig = 1:Ng
      Gr2 = ( Gm[ig]*rrl )^2
      Vprj[ig] = Gm[ig]^3 * exp(-0.5*Gr2)
    end

  end  # if l

  pre =  4 * pi^(5./4.) * sqrt( 2.^(l+1) * rrl^(2*l+3) / Ω )

  Vprj[:] = pre * Vprj[:]

  return Vprj

end


text:
"""
    multiple_target_state_extraction(b::GaussianMixture, th::Real)

Extracts targets whose weights (x.w) are above threshold.

Arguments:
    `b::GaussianMixture`: Set of Gaussian Mixtures
    `th::Real`: Threshold on weights. Above this threshold,
    state estimate is extracted

Returns:
    X: Multi-Target State Estimate
"""
function multiple_target_state_extraction(b::GaussianMixture, th::Real)
    inds = b.w .> th
    N = sum(inds)
    m = b.μ[inds]
    w = b.w[inds]

    X = []
    for i = 1:N
        tmp = [m[i] for j = 1:round(w[i])]
        X = vcat(X, tmp)
    end

    return X
end


text:
"""
    Region(SPParams, TMParams; recurrent= true)
A Region is a set of neurons arranged in minicolumns with 2 input gates:
- _proximal_ (size: number of minicolumns Nc): feedforward input that activates some of the region's minicolumns
- _distal_ (size: number of neurons Nn): contextual input that causes individual neurons to become predictive and fire preferentially on the next time step

Each gate accepts 1 or more input SDRs of the correct size, which will be ORed together.
If `recurrent= true`, the distal gate implicitly receives the same layer's state from the last time step as input (in addition to any SDR given explicitly).

The activation of minicolumns resulting from the proximal synapses is defined by the [`SpatialPooler`](@ref) algorithm,
while the activation of individual neurons is defined by the [`TemporalMemory`](@ref).

# Example

```julia
r= Region(SPParams(), TMParams(), recurrent= true)
feedforward_input= bitrand(Nc(r))

# activate based on feedforward input and the region's previous activation
active, predictive= r(feedforward_input)
# learn (adapt synapses)
step!(r, feedforward_input)

# use explicit contextual input eg from another region, that will be ORed to the region's state
distal_input= sprand(Bool,Nn(r), 0.003)|> BitVector
active, predictive= r(feedforward_input, distal_input)
```
"""
struct Region
  sp::SpatialPooler
  tm::TemporalMemory
end
Region(spp::SPParams,tmp::TMParams; recurrent=true, distal_input_size=0)= begin
  @assert prod(spp.szₛₚ) == tmp.Nc "Number of minicolumns in Spatial Pooler and Temporal Memory in the same region must match"
  Region(SpatialPooler(spp), TemporalMemory(tmp, recurrent=recurrent, distal_input_size=distal_input_size))
end
Nc(r::Region)= r.tm.params.Nc
Nₙ(r::Region)= Nₙ(r.tm)

# distal:= implicit recurrent connections only
(r::Region)(proximal, distal=falses(0))= @chain gateCombine(proximal) begin
  r.sp
  r.tm(_, gateCombine(distal))
end

step!(r::Region, proximal, distal=falses(0))= @chain gateCombine(proximal) begin
  step!(r.sp, _)
  step!(r.tm, _, gateCombine(distal))
end

reset!(r::Region)= begin
  reset!(r.sp)
  reset!(r.tm)
end

"""
    connectContext!(output::Region, input::Region)

Connect contextually 2 regions `output -> input` by wiring the activation of
`output` to the context gate of `input`. Only `input` changes.
"""
#connectContext!(input::Region, output::Region)

bitVcat(a::BitArray, b)= Base.typed_vcat(typeof(a), a,b)
gateCombine(x::Vector{CellActivity})= reduce((a,b)-> bitVcat(a,b), x)
gateCombine(x)= x


text:
mutable struct ShowGraphvizConfig
    dot::Cmd
    dot_option::Cmd
    pygmentize::Cmd
    pygmentize_option::Cmd
end

const CONFIG = ShowGraphvizConfig(`dot`, ``, `pygmentize`, ``)

# `Graphviz_jll.dot`-like interface
function dot(f)
    f(CONFIG.dot)
end

const SUPPORTS_GRAPHVIZ = Dict{Cmd,Bool}()
supports_graphviz(cmd) =
    get!(SUPPORTS_GRAPHVIZ, cmd) do
        _supports_graphviz(cmd)
    end

function _supports_graphviz(cmd)
    io = IOBuffer()
    try
        communicate(`$cmd -L`; stdin = devnull, stdout = io)
    catch err
        @error "Igoring error from `$cmd`" exception = (err, catch_backtrace())
        return false
    end
    occursin("graphviz", String(take!(io)))
end

highlight_dot(io::IO, input::AbstractString) = highlight_dot(io, IOBuffer(input))
function highlight_dot(io::IO, input::IO)
    cmd = CONFIG.pygmentize
    if !supports_graphviz(cmd) || get(io, :color, false) === false
        write(io, input)
    else
        cmd = `$cmd $(CONFIG.pygmentize_option)`
        cmd = `$cmd -l dot`
        communicate(cmd; stdin = input, stdout = io)
    end
end


text:
using BilevelJuMP
using Documenter
using Literate
using Test

const _EXAMPLE_DIR = joinpath(@__DIR__, "src", "examples")

"""
    _include_sandbox(filename)
Include the `filename` in a temporary module that acts as a sandbox. (Ensuring
no constants or functions leak into other files.)
"""
function _include_sandbox(filename)
    mod = @eval module $(gensym()) end
    return Base.include(mod, filename)
end

function _file_list(full_dir, relative_dir, extension)
    return map(
        file -> joinpath(relative_dir, file),
        filter(file -> endswith(file, extension), sort(readdir(full_dir))),
    )
end

function link_example(content)
    edit_url = match(r"EditURL = \"(.+?)\"", content)[1]
    footer = match(r"^(---\n\n\*This page was generated using)"m, content)[1]
    content = replace(
        content, footer => "!!! info\n    [View this file on Github]($(edit_url)).\n\n" * footer
    )
    return content
end

function literate_directory(dir)
    rm.(_file_list(dir, dir, ".md"))
    for filename in _file_list(dir, dir, ".jl")
        # `include` the file to test it before `#src` lines are removed. It is
        # in a testset to isolate local variables between files.
        @testset "$(filename)" begin
            _include_sandbox(filename)
        end
        Literate.markdown(
            filename,
            dir;
            documenter = true,
            postprocess = link_example
        )
    end
    return
end

literate_directory(_EXAMPLE_DIR)

makedocs(
    modules = [BilevelJuMP],
    doctest  = false,
    clean    = true,
    format   = Documenter.HTML(
        mathengine = Documenter.MathJax2(),
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    sitename = "BilevelJuMP.jl",
    authors  = "Joaquim Garcia",
    pages   = [
        "Home" => "index.md",
        "Manual" => "manual.md",
        "Examples" => [
            joinpath("examples", f) for
            f in readdir(_EXAMPLE_DIR) if endswith(f, ".md")
        ]
    ]
)

deploydocs(
    repo = "github.com/joaquimg/BilevelJuMP.jl.git",
    push_preview = true,
)


text:
# This file is for functions that just need simple standard overloading.

## Helpers:

function nameddimsarray_result(original_nda, reduced_data, reduction_dims)
    L = dimnames(original_nda)
    return NamedDimsArray{L}(reduced_data)
end

# if reducing over `:` then results is a scalar
function nameddimsarray_result(original_nda, reduced_data, reduction_dims::Colon)
    return reduced_data
end

################################################
# Overloads

# 1 Arg
for (mod, funs) in (
    (:Base, (:sum, :prod, :maximum, :minimum, :extrema, :argmax, :argmin)),
    (:Statistics, (:mean, :std, :var, :median)),
)
    for fun in funs
        @eval function $mod.$fun(a::NamedDimsArray; dims=:, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(parent(a); dims=numerical_dims, kwargs...)
            return nameddimsarray_result(a, data, numerical_dims)
        end
    end
end

# 1 Arg, 2 Results
for (mod, funs) in ((:Base, (:findmax, :findmin)),)
    for fun in funs
        @eval function $mod.$fun(a::NamedDimsArray; dims=:, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(parent(a); dims=numerical_dims, kwargs...)
            map(data) do item  # value and index
                nameddimsarray_result(a, item, numerical_dims)
            end
        end
    end
end

# 1 Arg - no default for `dims` keyword
for (mod, funs) in ((:Base, (:cumsum, :cumprod, :sort, :sort!, :sortslices)),)
    for fun in funs
        @eval function $mod.$fun(a::NamedDimsArray; dims, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(parent(a); dims=numerical_dims, kwargs...)
            return nameddimsarray_result(a, data, numerical_dims)
        end

        # Vector case
        @eval function $mod.$fun(a::NamedDimsArray{L,T,1}; kwargs...) where {L,T}
            data = $mod.$fun(parent(a); kwargs...)
            return NamedDimsArray{dimnames(a)}(data)
        end
    end
end

if VERSION > v"1.1-"
    function Base.eachslice(a::NamedDimsArray{L}; dims, kwargs...) where {L}
        numerical_dims = dim(a, dims)
        slices = eachslice(parent(a); dims=numerical_dims, kwargs...)
        return Base.Generator(slices) do slice
            # For unknown reasons (something to do with hoisting?) having this in the
            # function passed to `Generator` actually results in less memory being allocated
            names = remaining_dimnames_after_dropping(L, numerical_dims)
            return NamedDimsArray(slice, names)
        end
    end
end

# 1 arg before - no default for `dims` keyword
for (mod, funs) in ((:Base, (:mapslices,)),)
    for fun in funs
        @eval function $mod.$fun(f, a::NamedDimsArray; dims, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(f, parent(a); dims=numerical_dims, kwargs...)
            return nameddimsarray_result(a, data, numerical_dims)
        end
    end
end

# 2 arg before
for (mod, funs) in ((:Base, (:mapreduce,)),)
    for fun in funs
        @eval function $mod.$fun(f1, f2, a::NamedDimsArray; dims=:, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(f1, f2, parent(a); dims=numerical_dims, kwargs...)
            return nameddimsarray_result(a, data, numerical_dims)
        end
    end
end

################################################
# Non-dim Overloads

for fun in (:(==), :isequal, :isapprox)
    @eval function Base.$fun(
        a::NamedDimsArray{La}, b::NamedDimsArray{Lb}; kw...,
    ) where {La,Lb}
        names_are_unifiable(La, Lb) || return false
        return $fun(parent(a), parent(b); kw...)
    end
end

# Array then perhaps other args
for (mod, funs) in ((:Base, (:zero, :one, :copy, :empty!, :push!, :pushfirst!)),)
    for fun in funs
        @eval function $mod.$fun(a::NamedDimsArray{L}, x...) where {L}
            data = $mod.$fun(parent(a), x...)
            return NamedDimsArray{L}(data)
        end
    end
end

# Two arrays
for (mod, funs) in ((:Base, (:sum!, :prod!, :maximum!, :minimum!)),)
    for fun in funs
        @eval begin

            function $mod.$fun(a::NamedDimsArray{L}, b::AbstractArray) where {L}
                data = $mod.$fun(parent(a), b)
                return NamedDimsArray{L}(data)
            end

            function $mod.$fun(a::AbstractArray, b::NamedDimsArray{L}) where {L}
                data = $mod.$fun(a, parent(b))
                newL = unify_names_shortest(L, ntuple(_ -> :_, ndims(a)))
                return NamedDimsArray{newL}(data)
            end

            function $mod.$fun(a::NamedDimsArray{La}, b::NamedDimsArray{Lb}) where {La,Lb}
                newL = unify_names_shortest(La, Lb)
                data = $mod.$fun(parent(a), parent(b))
                return NamedDimsArray{newL}(data)
            end

        end
    end
end

Base.pop!(A::NamedDimsArray) = pop!(parent(A))
Base.popfirst!(A::NamedDimsArray) = popfirst!(parent(A))

function Base.append!(A::NamedDimsArray{L,T,1}, B::AbstractVector) where {L,T}
    newL = unify_names(L, dimnames(B))
    data = append!(parent(A), unname(B))
    return NamedDimsArray{newL}(data)
end

################################################
# map, collect

Base.map(f, A::NamedDimsArray) = NamedDimsArray(map(f, parent(A)), dimnames(A))

for (T, S) in [
    (:NamedDimsArray, :AbstractArray),
    (:AbstractArray, :NamedDimsArray),
    (:NamedDimsArray, :NamedDimsArray),
]
    for fun in [:map, :map!]

        # Here f::F where {F} is needed to avoid ambiguities in Julia 1.0
        @eval function Base.$fun(f::F, a::$T, b::$S, cs::AbstractArray...) where {F}
            data = $fun(f, unname(a), unname(b), unname.(cs)...)
            new_names = unify_names(dimnames(a), dimnames(b), dimnames.(cs)...)
            return NamedDimsArray(data, new_names)
        end

    end

    @eval function Base.foreach(f::F, a::$T, b::$S, cs::AbstractArray...) where {F}
        data = foreach(f, unname(a), unname(b), unname.(cs)...)
        unify_names(dimnames(a), dimnames(b), dimnames.(cs)...)
        return nothing
    end
end

function Base.filter(f, A::NamedDimsArray{L,T,1}) where {L,T}
    return NamedDimsArray(filter(f, parent(A)), L)
end
Base.filter(f, A::NamedDimsArray{L,T,N}) where {L,T,N} = filter(f, parent(A))

# We overload collect on various kinds of `Generators` so that that can keep names.
function Base.collect(x::Base.Generator{<:NamedDimsArray{L}}) where {L}
    data = collect(Base.Generator(x.f, parent(x.iter)))
    return NamedDimsArray(data, L)
end

function Base.collect(
    x::Base.Generator{<:Iterators.Enumerate{<:NamedDimsArray{L}}},
) where {L}
    data = collect(Base.Generator(x.f, enumerate(parent(x.iter.itr))))
    return NamedDimsArray(data, L)
end

for Tup in (
    Tuple{<:NamedDimsArray,Vararg{Any}},
    Tuple{<:Any,<:NamedDimsArray,Vararg{Any}},
    Tuple{<:NamedDimsArray,<:NamedDimsArray,Vararg{Any}},
)
    @eval function Base.collect(x::Base.Generator{<:Iterators.ProductIterator{<:$Tup}})
        return collect_product(x)
    end
end

function collect_product(x)
    data = collect(Base.Generator(x.f, Iterators.product(unname.(x.iter.iterators)...)))
    all_names = tuple_cat(dimnames.(x.iter.iterators)...)
    return NamedDimsArray(data, all_names)
end


text:
# Tested in Julia 1.5.1
using DelimitedFiles
data = readdlm("./data/adventInput02.txt")
function PasswordPhilosophy(data; position = false)
    validPasswords = 0
    for (policy, letter, password) in eachrow(data)
        policy_ranges = parse.(Int, split(policy, "-"))
        policy_letter = letter[1:1]
        if position == false
            countLetter = count(x-> (x == only("$(policy_letter)")), password)
            if policy_ranges[1]<= countLetter <=policy_ranges[2]
                validPasswords += 1
            end
        elseif position == true
            positions = findall(x -> (x == only("$(policy_letter)")), password)
            isone = sum([n in positions for n in policy_ranges])
            if isone == 1
                validPasswords += 1
            end
        end
    end
    validPasswords
end
PasswordPhilosophy(data; position=false)
PasswordPhilosophy(data; position=true)


text:
using Base.Threads

type Distribution
    G::Tuple{Vararg{Array{Float64}}}
    d::Array{Float64}
    T::SparseMatrixCSC
end
type ModelSummary
    equations::Expr
    parameters::Dict
    displays#::Vector{Window}
end

type Model
    G::NGrid
    SP::Array{Float64, 2}
    ProbWeights::Array{Float64, 2}
    X::Array{Float64, 2}
    XP::Array{Float64, 2}
    F::Function
    Fall::Function
    Fouter::Function
    Finner::Function
    Fval::Array{Float64, 2}
    variables::Vector{Variable}
    distribution::Distribution
    temp
    summary::ModelSummary
end

display(M::Model) = display(M.variables)
length(M::Model) =  length(M.G)

function Model(foc, states, vars, params, B = Quadratic)
    foc1 = deepcopy(foc)
    variables, parameters, F, J = parsevars(foc, states, vars, params)
    G, S, SP, ProbWeights       = initS(State(variables), B)
    F1      = deepcopy(F)
    Fall, outfunc, infunc    = buildSolver(F, J, G, variables)
    X       = [S zeros(length(G), length(Policy(variables)) + length(Static(variables)))]
    XP      = zeros(length(G) * size(ProbWeights, 2), length(Future(variables)))

    tG      = ndgrid(Vector{Float64}[s.x for s in State(variables)]...)
    distribution = Distribution(tG, zeros(size(tG[1])) + 1 / length(tG[1]), spzeros(length(tG[1]), length(tG[1])))

    F1       = :($(gensym(:F))(M::Model) = @fastmath $(buildfunc(F1, :(M.Fval))))

    M = Model(G, SP, ProbWeights, X, XP, eval(F1), eval(Fall), eval(outfunc), eval(infunc), zeros(length(G), length(Policy(variables))), variables, distribution, (Fall, outfunc, infunc), ModelSummary(foc1, parameters, []))

    precompile(M.Fall, (Model,))

    initX(M)
    initD(M)

    return M
end

(M::Model)(s::Symbol, x::Array{Float64, 2}) = M.G(M[s, 0], x)

function getindex(M::Model, x::Symbol, t::Int)
    if t ≤ 0
        for i ∈ 1:length(M.variables)
            if M.variables[i].name == x && timeof(M.variables[i]) == t
                return M.X[:, i]
            end
        end
    else
        for i ∈ 1:length(Stochastic(M.variables))
            if Stochastic(M.variables)[i].name == x
                return M.SP[:, length(Endogenous(M.variables)) + i]
            end
        end
        for i ∈ 1:length(Future(M.variables))
            if Future(M.variables)[i].name == x && timeof(Future(M.variables)[i]) == t
                return M.XP[:, i]
            end
        end
    end
end

Expect(M::Model, x::Symbol) = vec(sum(reshape(M[x, 1], size(M.ProbWeights)) .* M.ProbWeights, 2))

function setindex!(M::Model, x, v::Symbol, t::Int)
    for i ∈ 1:length(M.variables)
        if M.variables[i].name == v && timeof(M.variables[i]) == t
            M.X[:, i] = x
        end
    end
end

getindex(M::Model, v::Symbol) = M.variables[v]

function buildfunc(ex::Expr, targ, t = 0)
    F = :(for i = 1:length(M) end)
    for i = 1:length(ex.args)
        push!(F.args[2].args, :($targ[i, $(i + t)] = $(ex.args[i])))
    end
    return F
end

@inline pow(a, b) = exp(b * log(a))
function subspow!(x)
    if isa(x, Expr)
        if x.head == :call && x.args[1] == :^
            x.args[1] = :pow
            for i = 2:length(x.args)
                x.args[i] = subspow!(x.args[i])
            end
        else
            for i = 1:length(x.args)
                x.args[i] = subspow!(x.args[i])
            end
        end
    end
    return x
end
# function to update policy variables
# 1 - defineds equation errors
# 2 - finds jacobian
# 3 - 1 Newton step (x = x-ϕ*J(x)\F(x))
# 4 - strip ~long repeated expressions

function buildSolver(F, J, G, variables)
    Fexpr = Expr(:function, :($(gensym("updateX"))(M::Model, ϕ = 0.8)), Expr(:block, :(Fval = zeros($(length(Policy(variables))))), :(Jval = zeros($(length(Policy(variables))), $(length(Policy(variables)))))))
    preloopblock = Expr(:block)
    mainloopblock = Expr(:block)
    for i = 1:length(F.args)
        push!(mainloopblock.args, :(Fval[$i] = $(F.args[i])))
    end

    for i = 1:length(J.args)
        if isa(J.args[i], Number) && J.args[i] == 0
            J.args[i] != 0 && push!(Fexpr.args[2].args, :(Jval[$i] = $(J.args[i])))
        else
            push!(mainloopblock.args, :(Jval[$i] = $(J.args[i])))
        end
    end

    push!(mainloopblock.args, :(A_ldiv_B!(lufact(Jval), Fval)))
    push!(mainloopblock.args, Expr(:for, :(j = 1:$(length(Policy(variables)))), Expr(:macrocall, Symbol("@inbounds"), :(M.X[i, $(length(State(variables))) + j] -= ϕ * Fval[j]))))
    mainloop = Expr(:macrocall, Symbol("@threadsfixed"), :([Fval,Jval]), Expr(:for, :(i = 1:(length(M))), mainloopblock))

    push!(Fexpr.args[2].args, :(N = length(M)))
    subs!(mainloop, :(length(M)) => (:N))
    push!(Fexpr.args[2].args, mainloop)

    outfunc = Expr(:function, :($(gensym("OuterFunc"))(M::Model, ϕ = 0.8)), Expr(:block, :(Fval = zeros($(length(Policy(variables))))), :(Jval = zeros($(length(Policy(variables))), $(length(Policy(variables))))), :(N = length(M)), Expr(:macrocall, Symbol("@threadsfixed"), :([Fval,Jval]), Expr(:for, :(i = 1:(length(M))), :(M.Finner(M, i, N, Fval, Jval, ϕ))))))
    infunc  = Expr(:function, :($(gensym("InnerFunc"))(M, i, N, Fval, Jval, ϕ)), mainloopblock)

    return Fexpr, outfunc, infunc
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
# SF6 defined as "linear gas" or gas 4 in PAGE 2009; equations are the same as for gas 3 (CH4) in PAGE2002
@defcomp LGemissions begin
    region = Index()

    # global emissions
    e_globalLGemissions = Variable(index=[time], unit="Mtonne/year")
    # baseline emissions
    e0_baselineLGemissions = Parameter(index=[region], unit="Mtonne/year")
    # regional emissions
    e_regionalLGemissions = Variable(index=[time,region], unit="Mtonne/year")
    # growth rate by region
    er_LGemissionsgrowth = Parameter(index=[time,region], unit="%")

    # read in counterfactual GDP in absence of growth effects (gdp_leveleffects) and actual GDP
    gdp = Parameter(index=[time, region], unit="\$M")
    gdp_leveleffect   = Parameter(index=[time, region], unit="\$M")
    emfeed_emissionfeedback = Parameter(unit="none", default=1.)

    function run_timestep(p, v, d, t)

        # eq.4 in Hope (2006) - regional LG emissions as % change from baseline
        for r in d.region
            v.e_regionalLGemissions[t,r] = p.er_LGemissionsgrowth[t,r] * p.e0_baselineLGemissions[r] / 100

            # rescale emissions based on GDP deviation from original scenario pathway
            if p.emfeed_emissionfeedback == 1.
                v.e_regionalLGemissions[t,r] = v.e_regionalLGemissions[t,r] * (p.gdp[t,r] / p.gdp_leveleffect[t,r])
            end
        end

        # eq. 5 in Hope (2006) - global LG emissions are sum of regional emissions
        v.e_globalLGemissions[t] = sum(v.e_regionalLGemissions[t,:])
    end
end


text:
export rnnOrderGraph

function rnnOrderGraph(L,T)
  orderGraph=zeros(Int64,L*T,L*T)
  for i=1:L*T
    if i%L==1
      orderGraph[i,i+1]=1
    elseif i%L!=0
      orderGraph[i,i+1]=1
      if i+L<L*T
        orderGraph[i,i+L]=1
      end
    end
  end
  return orderGraph
end


text:
# ---
# title: Dynamic Range Compression
# id: demo_ppdrc
# cover: assets/ppdrc.png
# author: Peter Kovesi
# date: 2018-10-26
# ---

# An example using the 16 bit M51 image.  Phase preserving dynamic range compression allows
# the scale of analysis to be controlled.  Here we process the image at wavelengths up to 50
# pixels and up to 200 pixels.  Longer wavelengths allow larger structures to be seen. Small
# wavelengths allow fine structures to be seen.  Note the image size is (510, 320).

using TestImages
using Images
using ImageContrastAdjustment
using ImagePhaseCongruency

img = float64.(testimage("m51"))

## Histogram equalization for reference (with a very large number of bins!)
img_histeq = histeq(img, 100_000)

## Phase presserving dynamic range compression at cutoff wavelengths of 50 and
## 200 pixels.  Note we scale the image because its raw values are between 0 and
## 1, see the help information for ppdrc() for details.
scale = 1e4
img_ppdrc1 = ppdrc(img*scale, 50)
img_ppdrc2 = ppdrc(img*scale, 200)

mosaic(
    adjust_histogram(img, LinearStretching()),
    adjust_histogram(img_histeq, LinearStretching()),
    adjust_histogram(img_ppdrc1, LinearStretching()),
    adjust_histogram(img_ppdrc2, LinearStretching()),
    nrow=1
)

# save cover image #src
isdir("assets") || mkdir("assets") #src
cropped_cover = adjust_histogram(centered(img_ppdrc1)[-128:127, -128:127], LinearStretching()) #src
save(joinpath("assets", "ppdrc.png"), cropped_cover) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
@testset "DenseBE GPU" begin
    ensemble_size = 4
    samples_per_model = 4
    input_dim = 5
    output_dim = 5
    rank = 1
    inputs = rand(Float32, input_dim, samples_per_model)
    layer = DenseBE(
        input_dim,
        output_dim,
        rank,
        ensemble_size;
        alpha_init = ones,
        gamma_init = ones,
    )
    layer = layer |> gpu
    batch_inputs = gpu(repeat(inputs, 1, ensemble_size))
    batch_outputs = layer(batch_inputs)
    # Do the computation in for loop to compare outputs
    layer = layer |> cpu
    loop_outputs = []
    for i = 1:ensemble_size
        perturbed_inputs = inputs .* layer.alpha[i]
        outputs = layer.layer(perturbed_inputs) .* layer.gamma[i]
        outputs = layer.ensemble_act.(outputs .+ layer.ensemble_bias[i])
        push!(loop_outputs, outputs)
    end
    loop_outputs = Flux.batch(loop_outputs)
    loop_outputs = reshape(loop_outputs, (output_dim, samples_per_model * ensemble_size))
    @test batch_outputs isa CuArray
    @test size(batch_outputs) == size(loop_outputs)
    @test isapprox(cpu(batch_outputs), loop_outputs, atol = 0.05)

    # Test gradients
    layer = gpu(DenseBE(2, 5, 1, 2))
    i = gpu(rand(2, 4))
    y = gpu(ones(5, 4))
    grads = gradient(params(layer)) do
        ŷ = layer(i)
        return Flux.logitcrossentropy(ŷ, y)
    end
    for param in params(layer)
        @test size(param) == size(grads[param])
        @test grads[param] isa CuArray
    end
end

@testset "ConvBE GPU" begin
    ensemble_size = 4
    samples_per_model = 4
    input_dim = 5
    output_dim = 10
    rank = 1
    inputs = rand(Float32, 10, 10, input_dim, samples_per_model)
    beconv = ConvBE(
        (5, 5),
        5 => 10,
        rank,
        ensemble_size,
        relu;
        alpha_init = ones,
        gamma_init = ones,
    )
    beconv = beconv |> gpu
    batch_inputs = gpu(repeat(inputs, 1, 1, 1, ensemble_size))
    batch_outputs = beconv(batch_inputs)
    # Do the computation in for loop to compare outputs
    beconv = beconv |> cpu
    loop_outputs = []
    for i = 1:ensemble_size
        perturbed_inputs = inputs .* beconv.alpha[i]
        outputs = beconv.layer(perturbed_inputs) .* beconv.gamma[i]
        outputs = beconv.ensemble_act.(outputs .+ beconv.ensemble_bias[i])
        push!(loop_outputs, outputs)
    end
    loop_outputs = Flux.batch(loop_outputs)
    loop_outputs_size = size(batch_outputs)
    loop_outputs = reshape(
        loop_outputs,
        (
            loop_outputs_size[1],
            loop_outputs_size[2],
            output_dim,
            samples_per_model * ensemble_size,
        ),
    )
    @test batch_outputs isa CuArray
    @test size(batch_outputs) == size(loop_outputs)
    @test isapprox(cpu(batch_outputs), loop_outputs, atol = 0.05)

    # Test gradients
    layer = gpu(ConvBE((5, 5), 3 => 6, 1, 4, relu))
    i = gpu(rand(32, 32, 3, 4))
    y = gpu(rand(28, 28, 6, 4))
    grads = gradient(params(layer)) do
        ŷ = layer(i)
        return Flux.logitcrossentropy(ŷ, y)
    end
    for param in params(layer)
        @test size(param) == size(grads[param])
        @test grads[param] isa CuArray
    end
end


text:
##############################################################################
##
## Solution method: coordinate gradient descent
## (also called alternative regressions, Gauss Seidel)
##
##############################################################################

##############################################################################
##
## Factor Model (no regressors)
##
##############################################################################

function fit!(::Type{Val{:gauss_seidel}},
    fp::FactorModel,
    fs::FactorSolution;
    maxiter::Integer  = 100_000,
    tol::Real = 1e-9,
    lambda::Real = 0.0
    )
    lambda == 0.0 || error("The alternative regression method only works with lambda = 0.0")

    # initialize
    iter = 0
    converged = true

    fp = FactorModel(deepcopy(fp.y), fp.sqrtw, fp.idrefs, fp.timerefs, rank(fp))
    idscale = Array(Float64, size(fs.idpool, 1))
    timescale = Array(Float64, size(fs.timepool, 1))

    for r in 1:rank(fp)
        fsr = view(fs, :, r)
        oldf_x = ssr(fp, fsr)
        iter_inner = 0
        while iter_inner < maxiter
            iter_inner += 1
            update!(Val{:gauss_seidel}, fp.y, fp.sqrtw, fp.idrefs, fp.timerefs, fsr.idpool,
                idscale, fsr.timepool)
            update!(Val{:gauss_seidel}, fp.y, fp.sqrtw, fp.timerefs, fp.idrefs, fsr.timepool,
                timescale, fsr.idpool)
            f_x = ssr(fp, fsr)
            if abs(f_x - oldf_x) < (abs(f_x) + tol) * tol
                break
            else
                oldf_x = f_x
            end
        end
        iter = max(iter_inner, iter)
        converged = (iter_inner < maxiter) && converged
        if r < rank(fp)
            rescale!(fsr)
            subtract_factor!(fp, fsr)
        end
    end
    return rescale(fs), iter, converged
end

##############################################################################
##
## Interactive FixedEffectModel
##
##############################################################################

function fit!(::Type{Val{:gauss_seidel}},
    fp::InteractiveFixedEffectsModel,
    fs::InteractiveFixedEffectsSolution;
    maxiter::Integer = 100_000,
    tol::Real = 1e-9,
    lambda::Real = 0.0)
    lambda == 0.0 || error("The alternative regression method only works with lambda = 0.0")

    N = size(fs.idpool, 1)
    T = size(fs.timepool, 1)

    #qr fact factorization cannot divide in place for now
    crossx = cholfact!(At_mul_B(fp.X, fp.X))
    M = crossx \ fp.X'

    yoriginal = deepcopy(fp.y)
    fp = InteractiveFixedEffectsModel(deepcopy(fp.y), fp.sqrtw, fp.X, fp.idrefs, fp.timerefs, rank(fp))

    idscale = Array(Float64, size(fs.idpool, 1))
    timescale = Array(Float64, size(fs.timepool, 1))

    # starts loop
    converged = false
    iterations = maxiter
    iter = 0
    f_x = Inf
    oldf_x = Inf
    while iter < maxiter
        iter += 1
        (f_x, oldf_x) = (oldf_x, f_x)

        # Given beta, compute incrementally an approximate factor model
        copy!(fp.y, yoriginal)
        BLAS.gemm!('N', 'N', -1.0, fp.X, fs.b, 1.0, fp.y)
        for r in 1:rank(fp)
            fsr = view(fs, :, r)
            update!(Val{:gauss_seidel}, fp.y, fp.sqrtw, fp.idrefs, fp.timerefs, fsr.idpool, idscale, fsr.timepool)
            update!(Val{:gauss_seidel}, fp.y, fp.sqrtw, fp.timerefs, fp.idrefs, fsr.timepool, timescale, fsr.idpool)
            subtract_factor!(fp, fsr)
        end

        # Given factor model, compute beta
        copy!(fp.y, yoriginal)
        subtract_factor!(fp, fs)
        ## corresponds to Gauss Niedel with acceleration
        scale!(fs.b, -0.5)
        BLAS.gemm!('N', 'N', 1.5, M, fp.y, 1.0, fs.b)

        # Check convergence
        BLAS.gemm!('N', 'N', -1.0, fp.X, fs.b, 1.0, fp.y)
        f_x = sumabs2(fp.y)
        if abs(f_x - oldf_x) < (abs(f_x) + tol) * tol
            converged = true
            iterations = iter
            break
        end
    end
    return rescale(fs), iterations, converged
end

##############################################################################
##
## update! by alternative regressions
##
##############################################################################

function update!{R1, R2}(::Type{Val{:gauss_seidel}},
    y::Vector{Float64},
    sqrtw::AbstractVector{Float64},
    p1refs::Vector{R1},
    p2refs::Vector{R2},
    p1::AbstractVector{Float64},
    p1scale::Vector{Float64},
    p2::AbstractVector{Float64}
    )
    fill!(p1, zero(Float64))
    fill!(p1scale, zero(Float64))
    @inbounds @simd for i in 1:length(y)
        p1i = p1refs[i]
        xi = sqrtw[i] * p2[p2refs[i]]
        p1[p1i] += xi * y[i]
        p1scale[p1i] += abs2(xi)
    end
    @inbounds @simd for i in 1:length(p1)
        s = p1scale[i]
        if s > 0
            p1[i] /= s
        end
    end
end

## compute sum of squared residuals
function ssr(fm::FactorModel, fs::FactorSolution{1})
    out = zero(Float64)
    @inbounds @simd for i in 1:length(fm.y)
        out += abs2(fm.y[i] - fm.sqrtw[i] * fs.idpool[fm.idrefs[i]] * fs.timepool[fm.timerefs[i]])
    end
    return out
end


text:
# copied from https://github.com/matago/TSPLIB.jl

function readCVRPLIB(name::String; add_dummy=false)
    vrp, sol = download_cvrp(name)
    raw = read(vrp, String)
    return _generateCVRP(raw; add_dummy=add_dummy), vrp, sol
end

function readCVRP(path::String; add_dummy=false)
    raw = read(path, String)
    return _generateCVRP(raw; add_dummy=add_dummy)
end

function _generateCVRP(raw::String; add_dummy=false)
    _dict = TSPLIB.keyextract(raw, cvrp_keys)

    name = _dict["NAME"]
    dimension = parse(Int, _dict["DIMENSION"])
    weight_type = _dict["EDGE_WEIGHT_TYPE"]

    if haskey(_dict, "DEPOT_SECTION")
        depot = parse.(Int, split(_dict["DEPOT_SECTION"]))[1]
    else
        depot = 1
        @warn "$name does not have 'DEPOT_SECTION'. Depot is set to node 1."
    end
    capacity = parse.(Int, _dict["CAPACITY"])

    check_dimension = add_dummy ? dimension + 1 : dimension
    dummy = add_dummy ? dimension + 1 : depot

    customers = setdiff(1:(dimension+1), [depot, dummy])

    if weight_type == "EXPLICIT" && haskey(_dict, "EDGE_WEIGHT_SECTION")
        explicits = parse.(Float64, split(_dict["EDGE_WEIGHT_SECTION"]))
        weights = TSPLIB.explicit_weights(_dict["EDGE_WEIGHT_FORMAT"], explicits)
        #Push display data to nodes if possible
        if haskey(_dict,"DISPLAY_DATA_SECTION")
            coords = parse.(Float64, split(_dict["DISPLAY_DATA_SECTION"]))
            n_r = convert(Integer, length(coords) / dimension)
            nodes = reshape(coords, (n_r, dimension))'[:,2:end]
            dxp = true
        else
            nodes = zeros(dimension, 2)
        end

        # Expanding weights for dummy depot
        if add_dummy
            weights = [weights weights[:, depot]]
            weights = [weights; weights[depot, :]']
        end
        # No coordinate information
        coordinates = Matrix{Float64}(undef, 0, 0)

    elseif haskey(_dict, "NODE_COORD_SECTION")
        coords = parse.(Float64, split(_dict["NODE_COORD_SECTION"]))
        n_r = convert(Integer, length(coords) / dimension)
        coordinates = reshape(coords, (n_r, dimension))'[:, 2:end]

        # copy depot and create a dummy
        if add_dummy
            coordinates = [coordinates; coordinates[depot, :]']
        end

        weights = Int.(TSPLIB.calc_weights(_dict["EDGE_WEIGHT_TYPE"], coordinates))
        @assert size(coordinates) == (check_dimension, 2)
    end

    @assert size(weights) == (check_dimension, check_dimension)

    distance_limit = Inf
    if haskey(_dict, "DISTANCE")
        distance_limit = parse(Float64, _dict["DISTANCE"])
    end

    service_time = 0.0
    if haskey(_dict, "SERVICE_TIME")
        service_time = parse(Float64, _dict["SERVICE_TIME"])
    end

    if haskey(_dict, "DEMAND_SECTION")
        demand_data = parse.(Int, split(_dict["DEMAND_SECTION"]))
        n_r = convert(Integer, length(demand_data) / dimension)
        demands_data = reshape(demand_data, (n_r, dimension))'[:, 2:end]
        demands = dropdims(demands_data, dims=2)

        # copy depot and create a dummy
        if add_dummy
            demands = [demands; demands[depot]]
        end

        @assert length(demands) == check_dimension
    end

    return CVRP(
        _dict["NAME"],
        dimension,
        _dict["EDGE_WEIGHT_TYPE"],
        weights,
        capacity,
        distance_limit,
        service_time,
        coordinates,
        demands,
        depot,
        dummy,
        customers
    )
end


text:
"""

    readRequiredFilenameARG(myARGS, currfilename, callingfunc)

Read from the list of command-line arguments `myARGS` the file name prefixed by
`currfilename` and an equals sign. Throw an error (using `callingfunc` to emphasize where it
is needed) if that argument is missing.

See also: [`readOptionalFilenameARG`](@ref), [`readGenvstr`](@ref),
[`readOptionalBoolcondARG`](@ref), [`readOptionalFloatType`](@ref),
[`readfreqlist!`](@ref), [`readklist`](@ref)

"""
function readRequiredFilenameARG(myARGS, currfilename::AbstractString,
                                 callingfunc::AbstractString)

    curridx = findfirst(map(u -> startswith(u, string(currfilename, "=")), myARGS));
    if (curridx == nothing)
        error(string("In", callingfunc, ", ", currfilename, "=\${", currfilename, "} is a ",
                     " REQUIRED argument"));
    end
    return replace(myARGS[curridx], string(currfilename, "=") => "");

end

"""

    readOptionalFilenameARG(myARGS, currfilename)

Read from the list of command-line arguments `myARGS` the file name prefixed by
`currfilename` and an equals sign. Return a blank string if that argument is missing.

See also: [`readRequiredFilenameARG`](@ref), [`readGenvstr`](@ref),
[`readOptionalBoolcondARG`](@ref), [`readOptionalFloatType`](@ref),
[`readfreqlist!`](@ref), [`readklist`](@ref)

"""
function readOptionalFilenameARG(myARGS, currfilename::AbstractString)

    curridx = findfirst(map(u -> startswith(u, string(currfilename, "=")), myARGS));
    if (curridx != nothing)
        return replace(myARGS[curridx], string(currfilename, "=") => "");
    end
    return "";

end

"""

    readGenvstr(myARGS)

Read from the list of command-line arguments `myARGS` the environmental Green's function
prefixed by `\"Genv=\"`. Return `\"PEC\"` (representing a PEC plane) only if that is
explicitly found, otherwise return `\"VAC\"` by default.

See also: [`readRequiredFilenameARG`](@ref), [`readOptionalFilenameARG`](@ref),
[`readOptionalBoolcondARG`](@ref), [`readOptionalFloatType`](@ref),
[`readfreqlist!`](@ref), [`readklist`](@ref)

"""
function readGenvstr(myARGS)

    if (findfirst(myARGS .== "Genv=PEC") != nothing)
        return "PEC";
    end
    return "VAC";

end

"""

    readOptionalBoolcondARG(myARGS, currvarname, outfilename)

Read from the list of command-line arguments `myARGS` optional boolean arguments labeled by
`currvarname`, and modify `outfilename` as appropriate.

See also: [`readRequiredFilenameARG`](@ref), [`readOptionalFilenameARG`](@ref),
[`readGenvstr`](@ref), [`readOptionalFloatType`](@ref),
[`readfreqlist!`](@ref), [`readklist`](@ref)

"""
function readOptionalBoolcondARG(myARGS, currvarname::AbstractString,
                                 outfilename::AbstractString)

    currboolvar = (findfirst(myARGS .== currvarname) != nothing);
    outfilename = currboolvar ? string(currvarname, "_", outfilename) : outfilename;
    return currboolvar, outfilename;

end

"""

    readOptionalFloatType(myARGS)

Read from the list of command-line arguments `myARGS` optional arguments labeling the
floating-point type. Return `Float64` by default if nothing else is found.

See also: [`readRequiredFilenameARG`](@ref), [`readOptionalFilenameARG`](@ref),
[`readGenvstr`](@ref), [`readOptionalBoolcondARG`](@ref),
[`readfreqlist!`](@ref), [`readklist`](@ref)

"""
function readOptionalFloatType(myARGS)

    FT = Float64;
    curridx = findfirst(map(u -> startswith(u, "FloatType="), myARGS));
    if (curridx != nothing)
        currsymb = Symbol(replace(myARGS[curridx], "FloatType=" => ""));
        if (isdefined(Main, currsymb))
            currFT = getfield(Main, currsymb);
            if (isa(currFT, DataType) && currFT<:AbstractFloat)
                FT = currFT;
            end
        end
    end
    return FT;

end

"""

    readfreqlist!(freqlist, myARGS, callingfunc)

Read from the list of command-line arguments `myARGS` the list of frequencies, modifying
`freqlist` \"in-place\" (actually appending as needed), and using the name `callingfunc`
to be more informative when throwing errors.

!!! note

    Any frequencies will be converted to their absolute values. Frequencies are not allowed
    to be exactly zero.

See also: [`readRequiredFilenameARG`](@ref), [`readOptionalFilenameARG`](@ref),
[`readGenvstr`](@ref), [`readOptionalBoolcondARG`](@ref),
[`readOptionalFloatType`](@ref), [`readklist`](@ref)

"""
function readfreqlist!(freqlist::Array{FT, 1}, myARGS,
                       callingfunc::AbstractString) where FT<:AbstractFloat

    freq = zero(eltype(freqlist));
    if (findfirst(map(u -> startswith(u, "freqlistfilename="), myARGS)) != nothing)

        freqliststrarray = strip.(chomp.(readlines(
                                             replace(myARGS[findfirst(map(u ->
                                                                  startswith(u, "freqlistfilename="),
                                                                  myARGS))],
                                             "freqlistfilename=" => ""))));
        filter!(x -> !isempty(x), freqliststrarray);
        append!(freqlist, abs.(parse.(Complex{FT}, freqliststrarray)));
        if (findfirst(freqlist .== zero(freqlist)) != nothing)
            error("In ", callingfunc, ", argument freqlistfilename=\${freqlistfilename}",
                  " cannot have any entries that are exactly 0.0");
        end
        print("Warning: In ", callingfunc, ", all entries in argument",
              " freqlistfilename=\${freqlistfilename} are converted to their");

    elseif (findfirst(map(u -> startswith(u, "freq="), myARGS)) != nothing)

        freq = abs(parse(Complex{FT},
                         replace(myARGS[findfirst(map(u -> startswith(u, "freq="),
                                                      myARGS))], "freq=" => "")));
        if (freq == zero(freq))
            error("In ", callingfunc, ", argument freq=\${freq} cannot be exactly 0.0");
        end
        print("Warning: In ", callingfunc, ", argument freq=\${freq} is converted to its");
        push!(freqlist, freq);

    else
        error("In ", callingfunc, ", EITHER arguments freq=\${freq} OR",
              " freqlistfilename=\${freqlistfilename} are REQUIRED");
    end
    print(" absolute value(s)");
    if (callingfunc == "AARMBEM_vdW.jl")
        println(", then multiplied by 1im");
    else
        println("");
    end

end

"""

    readklist(myARGS, callingfunc, numdims)

Read from the list of command-line arguments `myARGS` the list of Bloch wavevectors,
returning `zeros(FT, 1, 3)` if `numdims` is zero.

See also: [`readRequiredFilenameARG`](@ref), [`readOptionalFilenameARG`](@ref),
[`readGenvstr`](@ref), [`readOptionalBoolcondARG`](@ref),
[`readOptionalFloatType`](@ref), [`readfreqlist!`](@ref)

"""
function readklist(myARGS, callingfunc::AbstractString, numdims::Integer,
                   FT::DataType=Float64)

    if (!(FT <: AbstractFloat))
        error("FT must be a subtype of AbstractFloat");
    end
    klist = zeros(FT, 1, 3);
    if (numdims <= 0)

        return klist;

    elseif (findfirst(map(u -> startswith(u, "klistfilename="), myARGS)) != nothing)

        return readdlm(replace(myARGS[findfirst(map(u ->
                                                    startswith(u, "klistfilename="),
                                                    myARGS))],
                               "klistfilename=" => ""), ' ', FT);

    elseif (findfirst(map(u -> startswith(u, "kx="), myARGS)) != nothing ||
                findfirst(map(u -> startswith(u, "ky="), myARGS)) != nothing ||
                findfirst(map(u -> startswith(u, "kz="), myARGS)) != nothing)

        if (findfirst(map(u -> startswith(u, "kx="), myARGS)) != nothing)

            klist[1, 1] = parse(FT,
                                replace(myARGS[findfirst(map(u ->
                                                             startswith(u, "kx="),
                                                             myARGS))],
                                        "kx=" => ""));

        end
        if (findfirst(map(u -> startswith(u, "ky="), myARGS)) != nothing)

            klist[1, 2] = parse(FT,
                                replace(myARGS[findfirst(map(u ->
                                                             startswith(u, "ky="),
                                                             myARGS))],
                                        "ky=" => ""));

        end
        if (findfirst(map(u -> startswith(u, "kz="), myARGS)) != nothing)

            klist[1, 3] = parse(FT,
                                replace(myARGS[findfirst(map(u ->
                                                             startswith(u, "kz="),
                                                             myARGS))],
                                        "kz=" => ""));

        end

        return klist;

    else
        error("In ", callingfunc, ", if there are nontrivial periodic dimensions in",
              " \${periodicfilename}, EITHER arguments",
              " klistfilename=\${klistfilename} OR AT LEAST ONE OF kx=\${kx} OR",
              " ky=\${ky} OR kz=\${kz} are REQUIRED");
    end

end


text:
tests = ["utils",
         "classification",
         "perfeval",
         "crossval",
         "modeltune"]

for t in tests
    fp = "$t.jl"
    println("* running $fp ...")
    include(fp)
end


text:
module LogisticRegression1

using Plots
using Distributions
using LinearAlgebra

# Generate dummy data
n = 100;
x1 = transpose(rand(MvNormal([1., 2.], I), n));
x2 = transpose(rand(MvNormal([-2., -1.5], I), n));
y1 = repeat([1], n);
y2 = repeat([-1], n);
p = scatter(x1[:, 1], x1[:, 2], color = :blue)
scatter!(p, x2[:, 1], x2[:, 2], color = :red)

using Flux
using Flux.Tracker: gradient, update!, data
using Random
using Formatting: printfmt

# Logistic regression
Random.seed!(666);
α = param(rand());   # Bias
β = param(rand(2, 1));  # Weight
predict(X) = σ.(α .+ X*β)
X = [x1; x2];
y = [y1; y2];

function plot_contour(xlim, ylim, f)
    x = xlim[1]:0.5:xlim[2]
    y = ylim[1]:0.5:ylim[2]
    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    contour(x, y, f, fill=true, color=:bluesreds)
    scatter!(x1[:, 1], x1[:, 2],
            color = :red, label="Class 1")
    scatter!(x2[:, 1], x2[:, 2],
            color = :blue, label="Class 2")
end

xlim = ylim = (-5, 5)``
plot_contour(xlim, ylim,
                 (x, y) -> (first ∘ data ∘ predict)([x y]))

function loss(x, y)
    ŷ = predict(x)
    mean(Flux.binarycrossentropy.(ŷ, y))
end

printfmt("Initial loss: {:.4f}", (data ∘ loss)(X, y));

η = 0.01 # learning rate
epochs = 200
plot_interval = 10
anim = @animate for i=1:epochs
    for i in 1:plot_interval
        gs = gradient(() -> loss(X, y), Flux.params(α, β))
        update!(β, -η*gs[β])
        #update!(α, -η*gs[α])
    end
    display(plot_contour(xlim, ylim, (x, y) ->
                 (first ∘ data ∘ predict)([x y])))
end
gif(anim, "bluesreds.gif", fps = 60)

printfmt("Loss after training: {:.4f}", (data ∘ loss)(X, y));

end

max_stars_repo_name: SebastianCallh/julia-baby-steps

max_issues_repo_name: SebastianCallh/julia-baby-steps

max_forks_repo_name: SebastianCallh/julia-baby-steps


text:
using GraphMLDatasets
using LightGraphs
using SparseArrays
using Test

tests = [
    "planetoid",
    "cora",
    "ppi",
    "reddit",
    "qm7b",
    "entities",
    "utils",
]

@testset "GraphMLDatasets.jl" begin
    for t in tests
        include("$(t).jl")
    end
end


text:
"""
    raycolor(ray, scene, depth)

Trace a single ray and return color of ray as vector of RGB value.

# Arguments:
- `ray::Ray`: to-be traced ray
- `scene::Surface`: geometric scene as surface
- `depth::Integer`: how often the ray is allowed to scatter
"""
function raycolor(ray::Ray, scene::Surface, depth::Integer)
    rec = hit(scene, ray, 0.001, 1e6)[1]
    if rec != nothing
        scattered, attenuation = scatter(rec.material, ray, rec)
        if depth > 0 && scattered != nothing
            # pointwise multiplication between attenuation and
            # return value of recursive function call
            return attenuation .* raycolor(scattered, scene, depth - 1)
        else
            return zeros(3)
        end
    else
        # blue background sky
        unitdir = unitvector(ray.direction)
        t = 0.5*(unitdir[2] + 1)
        return (1 - t)*[1.0, 1.0, 1.0] + t*[0.5, 0.7, 1.0]
    end
end

"""
    renderimage(nx, ny, ns, scene, camera)

Render an image via raytracing; the image is returned as array of RGB values.

# Arguments
- `nx::Integer`: width of rendered image (pixels)
- `ny::Integer`: height of rendered image (pixels)
- `ns::Integer`: number of samples (rays) per pixel
- `scene::Surface`: geometric scene as surface
- `camera::Camera`: camera for generating rays
"""
function renderimage(nx::Integer, ny::Integer, ns::Integer, scene::Surface, camera::Camera)
    # fill image pixels
    img = zeros(RGB{Float32}, (nx, ny))
    for i in 1:nx
        for j in 1:ny
            color = zeros(3)
            for s in 1:ns
                # add a random offset for antialiasing
                u = (i-1 + rand()) / nx
                v = (j-1 + rand()) / ny
                ray = getray(camera, u, v)
                color += raycolor(ray, scene, 50)
            end
            color /= ns

            # take sqrt for gamma correction
            color = sqrt.(color)

            img[i, j] = RGB{Float32}(color[1], color[2], color[3])
        end
    end
    # up-down flip
    img = reverse(img, dims=2)

    return img
end


text:
using Color

export noborder,
       dotted,
       dashed,
       solid,
       strokewidth,
       border,
       bordercolor,
       roundcorner,
       shadow,
       fillcolor

render_color(c) = string("#" * hex(c))

## Borders

abstract BorderProperty

abstract StrokeStyle <: BorderProperty

@terms StrokeStyle begin
    noborder => NoStroke
    dotted => Dotted
    dashed => Dashed
    solid => Solid
end

propname(::StrokeStyle) = "Style"

name(::NoStroke) = "none"
name(::Solid) = "solid"
name(::Dotted) = "dotted"
name(::Dashed) = "dashed"

@api strokewidth => StrokeWidth <: BorderProperty begin
    arg(thickness::Length)
end

propname(::StrokeWidth) = "Width"
name(p::StrokeWidth) = p.thickness

@api bordercolor => BorderColor <: BorderProperty begin
    arg(color::ColorValue)
end

propname(::BorderColor) = "Color"
name(p::BorderColor) = render_color(p.color)

@api border => WithBorder{P <: BorderProperty} <: Tile begin
    typedarg(sides::AbstractArray=Side[])
    typedarg(prop::P)
    curry(tile::Tile)
end

render(t::WithBorder) =
    render(t.tile) &
        (isempty(t.sides) ? # Apply padding to all sides if none specified
                [:style => ["border" * propname(t.prop) => name(t.prop)]] :
                [:style => ["border" * name(s) * propname(t.prop) => name(t.prop) for s=t.sides]])

# Autopromote some types to border property
border(sides::AbstractArray, p::Length, x) =
    WithBorder(sides, strokewidth(p), x)
border(sides::AbstractArray, p::ColorValue, x) =
    WithBorder(sides, bordercolor(p), x)
border(p::Union(ColorValue, Length), x) =
    border(Side[], p, x)
border(sides::AbstractArray, ps::BorderProperty...) =
    t -> foldl((acc, p) -> border(sides, p, acc), t, ps)
border(ps::BorderProperty...) =
    border(Side[], ps...)

## RoundRects

@api roundcorner => RoundedRect <: Tile begin
    arg(radius::Length)
    curry(tile::Tile)
    kwarg(corners::AbstractArray=Corner[])
end

render(t::RoundedRect) =
    render(t.tile) &
        (isempty(t.corners) ? # Apply padding to all sides if none specified
                [:style => ["borderRadius" => t.radius]] :
                [:style => ["borderRadius" * name(c) => t.radius for c=t.corners]])

## Box shadow

@api shadow => Shadow <: Tile begin
    curry(tile::Tile)
    kwarg(inset::Bool=false)
    kwarg(offset::(Length, Length)=(0px, 0px))
    kwarg(blur_radius::Length=5px)
    kwarg(spread_radius::Length=5px)
    kwarg(color::ColorValue=color("black"))
end

## Fill color

@api fillcolor => FillColor <: Tile begin
    arg(color::ColorValue)
    curry(tile::Tile)
end

render(t::FillColor) =
    render(t.tile) & [:style => [:backgroundColor => render_color(t.color)]]


text:
# This file is a part of LegendDataTypes.jl, licensed under the MIT License (MIT).

export Geant4HDF5Input
struct Geant4HDF5Input <: AbstractLegendInput
    hdf5file::HDF5.File
end

struct GEARS_HDF5Input <: AbstractLegendInput
    hdf5file::HDF5.File
end

struct G4SIMPLE_HDF5Input <: AbstractLegendInput
    hdf5file::HDF5.File
end

Base.open(filename::AbstractString, ::Type{Geant4HDF5Input}) =
    Geant4HDF5Input(HDF5.h5open(filename, "r"))

Base.close(input::Geant4HDF5Input) = close(input.hdf5file)

function Base.getindex(input::Geant4HDF5Input, ::Colon)
    read(input)
end

function Base.read(input::Geant4HDF5Input)
    if haskey(input.hdf5file, "/default_ntuples/t/")
        read(GEARS_HDF5Input(input.hdf5file))
    elseif haskey(input.hdf5file, "/default_ntuples/g4sntuple/")
        read(G4SIMPLE_HDF5Input(input.hdf5file))
    end
end

function Base.read(input::GEARS_HDF5Input)

    function _get_indices_raw(g)
        len = filter(x->x!=0x00, g["e/pages"][:])
        indices_raw = []
        start = 0
        for i in len
            start+=1
            stop = start+i-1
            push!(indices_raw, collect(start:stop))
            start = stop
        end
        indices_raw
    end

    h5f = input.hdf5file
    g = h5f["default_ntuples/t/"]
    indices_raw = _get_indices_raw(g)
    h5_edep = g["e_data/pages"][:]
    indices_non_zero_energies  = filter(i-> sum(h5_edep[i]) > 0.0, indices_raw) # Filters events that don't deposit any energy
    indices = vcat(indices_non_zero_energies...)
    n_ind = length(indices)

    evtno = Int32.(vcat([ [j for i in 1:length(indices_non_zero_energies[j])] for j in 1:length(indices_non_zero_energies) ]...))
    detno = try Int32.(g["vlm_data/pages"][:][indices]) catch ; ones(Int32, n_ind) end
    thit = ( try Float32.(g["t0_data/pages"][:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"s"
    edep = ( try Float32.(g["e_data/pages"][:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"keV"

    x0 = try Float32.(g["x0_data/pages"][:][indices]) catch ; zeros(Float32, n_ind) end
    y0 = try Float32.(g["y0_data/pages"][:][indices]) catch ; zeros(Float32, n_ind) end
    z0 = try Float32.(g["z0_data/pages"][:][indices]) catch ; zeros(Float32, n_ind) end

    pos = [ SVector{3}(([ x0[i], y0[i], z0[i] ] .* u"mm")...) for i in 1:n_ind ]

    ekin = ( try Float32.(g["k_data/pages"][:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"keV"# kinetic Energy of track [keV]
    stp = try Int32.(g["stp_data/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # step number
    l = ( try Float32.(g["l_data/pages"][:][indices]) catch ; zeros(Float32, n_ind) end  ) .* u"mm"# step length
    mom = try Int32.(g["mom_data/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # parent id
    trk = try Int32.(g["trk_data/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # trk id
    pdg = try Int32.(g["pdg_data/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # particle id
    pro = try Int32.(g["pro_data/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # process id

    hits = TypedTables.Table(
        evtno = evtno,
        detno = detno,
        thit = thit,
        edep = edep,
        pos = pos,
        ekin = ekin,
        stp = stp,
        l = l,
        mom = mom,
        trk = trk,
        pdg = pdg,
        pro = pro,
    )
    return hits
end

function Base.read(input::G4SIMPLE_HDF5Input)
    function _get_indices_raw(g)
        len = filter(x->x!=0x00, g["e/pages"][:])
        indices_raw = []
        start = 0
        for i in len
            start+=1
            stop = start+i-1
            push!(indices_raw, collect(start:stop))
            start = stop
        end
        indices_raw
    end

    h5f = input.hdf5file
    g = h5f["default_ntuples/g4sntuple/"]
    # indices_raw = _get_indices_raw(g)
    h5_edep = g["event/pages"][:]
    # indices_non_zero_energies  = filter(i-> sum(h5_edep[i]) > 0.0, indices_raw) # Filters events that don't deposit any energy
    # indices = vcat(indices_non_zero_energies...)
    n_ind = length(h5_edep)
    indices = [i for i in 1:n_ind]

    evtno = try Int32.(g["event/pages"][:][indices]) catch ; ones(Int32, n_ind) end
    irep = try Int32.(g["iRep/pages"][:][indices]) catch ; ones(Int32, n_ind) end
    detno = try Int32.(g["volID/pages"][:][indices]) catch ; ones(Int32, n_ind) end
    thit = ( try Float32.(g["t/pages"][:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"ns"
    edep = ( try Float32.(g["Edep/pages"][:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"MeV"

    x0 = try Float32.(g["x/pages"][:][indices]) catch ; zeros(Float32, n_ind) end
    y0 = try Float32.(g["y/pages"][:][indices]) catch ; zeros(Float32, n_ind) end
    z0 = try Float32.(g["z/pages"][:][indices]) catch ; zeros(Float32, n_ind) end

    pos = [ SVector{3}(([ x0[i], y0[i], z0[i] ] .* u"mm")...) for i in 1:n_ind ]

    ekin = ( try Float32.(g["KE/pages"][:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"MeV" # kinetic Energy of track
    stp = try Int32.(g["step/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # step number
    mom = try Int32.(g["parentID/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # parent id
    trk = try Int32.(g["trackID/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # trk id
    pdg = try Int32.(g["pid/pages"][:][indices]) catch ; zeros(Int32, n_ind) end # particle id

    hits = TypedTables.Table(
        evtno = evtno,
        detno = irep,
        thit = thit,
        edep = edep,
        pos = pos,
        ekin = ekin,
        volID = detno,
        stp = stp,
        mom = mom,
        trk = trk,
        pdg = pdg,
    )
    return hits
end


text:
using Test
using OpenQuantumSystems
using LinearAlgebra
using SparseArrays

# https://github.com/JuliaIO/Suppressor.jl/blob/master/src/Suppressor.jl
macro suppress(block)
    quote
        if ccall(:jl_generating_output, Cint, ()) == 0
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            out_reader = @async read(out_rd, String)

            original_stderr = stderr
            err_rd, err_wr = redirect_stderr()
            err_reader = @async read(err_rd, String)

            # approach adapted from https://github.com/JuliaLang/IJulia.jl/pull/667/files
            logstate = Base.CoreLogging._global_logstate
            logger = logstate.logger
            if logger.stream == original_stderr
                new_logstate =
                    Base.CoreLogging.LogState(typeof(logger)(err_wr, logger.min_level))
                Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), new_logstate))
            end
        end

        try
            $(esc(block))
        finally
            if ccall(:jl_generating_output, Cint, ()) == 0
                redirect_stdout(original_stdout)
                close(out_wr)

                redirect_stderr(original_stderr)
                close(err_wr)

                if logger.stream == stderr
                    Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), logstate))
                end
            end
        end
    end
end

@testset "aggregateTools" begin

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1, mode2], 3, [2.0, 200.0])
    mol2 = Molecule([mode1], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    aggTools = AggregateTools(aggCore)

    @test aggTools.elIndices == [[1, 1], [2, 1], [1, 2]]
    @test aggTools.vibIndices == [[[1, 1], [1]], [[1, 1], [2]], [[1, 1], [3]], [[1, 2], [1]], [[1, 2], [2]], [[1, 2], [3]], [[1, 3], [1]], [[1, 3], [2]], [[1, 3], [3]], [[2, 1], [1]], [[2, 1], [2]], [[2, 1], [3]], [[2, 2], [1]], [[2, 2], [2]], [[2, 2], [3]], [[2, 3], [1]], [[2, 3], [2]], [[2, 3], [3]], [[3, 1], [1]], [[3, 1], [2]], [[3, 1], [3]], [[3, 2], [1]], [[3, 2], [2]], [[3, 2], [3]], [[3, 3], [1]], [[3, 3], [2]], [[3, 3], [3]]]
    @test aggTools.indices == [[[1, 1], [[1, 1], [1]]], [[1, 1], [[1, 1], [2]]], [[1, 1], [[1, 1], [3]]], [[1, 1], [[1, 2], [1]]], [[1, 1], [[1, 2], [2]]], [[1, 1], [[1, 2], [3]]], [[1, 1], [[1, 3], [1]]], [[1, 1], [[1, 3], [2]]], [[1, 1], [[1, 3], [3]]], [[1, 1], [[2, 1], [1]]], [[1, 1], [[2, 1], [2]]], [[1, 1], [[2, 1], [3]]], [[1, 1], [[2, 2], [1]]], [[1, 1], [[2, 2], [2]]], [[1, 1], [[2, 2], [3]]], [[1, 1], [[2, 3], [1]]], [[1, 1], [[2, 3], [2]]], [[1, 1], [[2, 3], [3]]], [[1, 1], [[3, 1], [1]]], [[1, 1], [[3, 1], [2]]], [[1, 1], [[3, 1], [3]]], [[1, 1], [[3, 2], [1]]], [[1, 1], [[3, 2], [2]]], [[1, 1], [[3, 2], [3]]], [[1, 1], [[3, 3], [1]]], [[1, 1], [[3, 3], [2]]], [[1, 1], [[3, 3], [3]]], [[2, 1], [[1, 1], [1]]], [[2, 1], [[1, 1], [2]]], [[2, 1], [[1, 1], [3]]], [[2, 1], [[1, 2], [1]]], [[2, 1], [[1, 2], [2]]], [[2, 1], [[1, 2], [3]]], [[2, 1], [[1, 3], [1]]], [[2, 1], [[1, 3], [2]]], [[2, 1], [[1, 3], [3]]], [[2, 1], [[2, 1], [1]]], [[2, 1], [[2, 1], [2]]], [[2, 1], [[2, 1], [3]]], [[2, 1], [[2, 2], [1]]], [[2, 1], [[2, 2], [2]]], [[2, 1], [[2, 2], [3]]], [[2, 1], [[2, 3], [1]]], [[2, 1], [[2, 3], [2]]], [[2, 1], [[2, 3], [3]]], [[2, 1], [[3, 1], [1]]], [[2, 1], [[3, 1], [2]]], [[2, 1], [[3, 1], [3]]], [[2, 1], [[3, 2], [1]]], [[2, 1], [[3, 2], [2]]], [[2, 1], [[3, 2], [3]]], [[2, 1], [[3, 3], [1]]], [[2, 1], [[3, 3], [2]]], [[2, 1], [[3, 3], [3]]], [[1, 2], [[1, 1], [1]]], [[1, 2], [[1, 1], [2]]], [[1, 2], [[1, 1], [3]]], [[1, 2], [[1, 2], [1]]], [[1, 2], [[1, 2], [2]]], [[1, 2], [[1, 2], [3]]], [[1, 2], [[1, 3], [1]]], [[1, 2], [[1, 3], [2]]], [[1, 2], [[1, 3], [3]]], [[1, 2], [[2, 1], [1]]], [[1, 2], [[2, 1], [2]]], [[1, 2], [[2, 1], [3]]], [[1, 2], [[2, 2], [1]]], [[1, 2], [[2, 2], [2]]], [[1, 2], [[2, 2], [3]]], [[1, 2], [[2, 3], [1]]], [[1, 2], [[2, 3], [2]]], [[1, 2], [[2, 3], [3]]], [[1, 2], [[3, 1], [1]]], [[1, 2], [[3, 1], [2]]], [[1, 2], [[3, 1], [3]]], [[1, 2], [[3, 2], [1]]], [[1, 2], [[3, 2], [2]]], [[1, 2], [[3, 2], [3]]], [[1, 2], [[3, 3], [1]]], [[1, 2], [[3, 3], [2]]], [[1, 2], [[3, 3], [3]]]]
    @test aggTools.indicesMap == [[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27], [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54], [55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81]]

    mode1 = Mode(0.2, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 2, [2.0, 200.0])
    mol2 = Molecule([mode1], 2, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    aggTools = AggregateTools(aggCore)

    ref = [1.0 0.0 0.0 0.0 0.7788007830714049 0.0 0.5506953149031838 0.0 0.7788007830714049 0.5506953149031838 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.7788007830714049 0.0 0.5506953149031838 -0.5506953149031837 0.3894003915357024 -0.0 0.0; 0.0 0.0 1.0 0.0 -0.5506953149031837 -0.0 0.3894003915357024 0.0 0.0 0.0 0.7788007830714049 0.5506953149031838; 0.0 0.0 0.0 1.0 -0.0 -0.5506953149031837 0.0 0.3894003915357024 -0.0 0.0 -0.5506953149031837 0.3894003915357024; 0.7788007830714049 0.0 -0.5506953149031837 -0.0 1.0 0.0 0.0 0.0 0.6065306597126334 0.42888194248035344 -0.42888194248035333 -0.3032653298563167; 0.0 0.7788007830714049 -0.0 -0.5506953149031837 0.0 1.0 0.0 0.0 -0.42888194248035333 0.30326532985631666 0.30326532985631666 -0.21444097124017664; 0.5506953149031838 0.0 0.3894003915357024 0.0 0.0 0.0 1.0 0.0 0.42888194248035344 0.30326532985631677 0.30326532985631666 0.2144409712401767; 0.0 0.5506953149031838 0.0 0.3894003915357024 0.0 0.0 0.0 1.0 -0.3032653298563167 0.2144409712401767 -0.21444097124017664 0.15163266492815833; 0.7788007830714049 -0.5506953149031837 0.0 -0.0 0.6065306597126334 -0.42888194248035333 0.42888194248035344 -0.3032653298563167 1.0 0.0 0.0 0.0; 0.5506953149031838 0.3894003915357024 0.0 0.0 0.42888194248035344 0.30326532985631666 0.30326532985631677 0.2144409712401767 0.0 1.0 0.0 0.0; 0.0 -0.0 0.7788007830714049 -0.5506953149031837 -0.42888194248035333 0.30326532985631666 0.30326532985631666 -0.21444097124017664 0.0 0.0 1.0 0.0; 0.0 0.0 0.5506953149031838 0.3894003915357024 -0.3032653298563167 -0.21444097124017664 0.2144409712401767 0.15163266492815833 0.0 0.0 0.0 1.0]
    @test 1e-14 > D(aggTools.FCfactors, ref)
    ref = [1.0 0.0 0.0 0.0 0.7788007830714049 0.0 0.5506953149031838 0.0 0.7788007830714049 0.5506953149031838 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.7788007830714049 0.0 0.5506953149031838 -0.5506953149031837 0.3894003915357024 0.0 0.0; 0.0 0.0 1.0 0.0 -0.5506953149031837 0.0 0.3894003915357024 0.0 0.0 0.0 0.7788007830714049 0.5506953149031838; 0.0 0.0 0.0 1.0 0.0 -0.5506953149031837 0.0 0.3894003915357024 0.0 0.0 -0.5506953149031837 0.3894003915357024; 0.7788007830714049 0.0 -0.5506953149031837 0.0 0.90979598956895 0.0 0.2144409712401768 0.0 0.6065306597126334 0.42888194248035344 -0.42888194248035333 -0.3032653298563167; 0.0 0.7788007830714049 0.0 -0.5506953149031837 0.0 0.90979598956895 0.0 0.2144409712401768 -0.42888194248035333 0.30326532985631666 0.30326532985631666 -0.21444097124017664; 0.5506953149031838 0.0 0.3894003915357024 0.0 0.2144409712401768 0.0 0.4548979947844751 0.0 0.42888194248035344 0.30326532985631677 0.30326532985631666 0.2144409712401767; 0.0 0.5506953149031838 0.0 0.3894003915357024 0.0 0.2144409712401768 0.0 0.4548979947844751 -0.3032653298563167 0.2144409712401767 -0.21444097124017664 0.15163266492815833; 0.7788007830714049 -0.5506953149031837 0.0 0.0 0.6065306597126334 -0.42888194248035333 0.42888194248035344 -0.3032653298563167 0.90979598956895 0.2144409712401768 0.0 0.0; 0.5506953149031838 0.3894003915357024 0.0 0.0 0.42888194248035344 0.30326532985631666 0.30326532985631677 0.2144409712401767 0.2144409712401768 0.4548979947844751 0.0 0.0; 0.0 0.0 0.7788007830714049 -0.5506953149031837 -0.42888194248035333 0.30326532985631666 0.30326532985631666 -0.21444097124017664 0.0 0.0 0.90979598956895 0.2144409712401768; 0.0 0.0 0.5506953149031838 0.3894003915357024 -0.3032653298563167 -0.21444097124017664 0.2144409712401767 0.15163266492815833 0.0 0.0 0.2144409712401768 0.4548979947844751]
    @test 1e-14 > D(aggTools.FCproduct, ref)

    @test aggTools.bSystemSize == 3
    @test aggTools.bBathSize == 4
    @test aggTools.bSize == 12

    @test aggTools.basisSystem == GenericBasis([3])
    @test aggTools.basisBath == GenericBasis([4])
    @test aggTools.basis == GenericBasis([12])

    FCfactors = getFranckCondonFactors(aggCore)
    @test 1e-14 > D(aggTools.FCfactors, FCfactors)

end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
mutable struct HistorySDEIntegrator{algType,IIP,uType,tType,SolType,CacheType} <: AbstractSDEIntegrator{algType,IIP,uType,tType}
    sol::SolType
    u::uType
    t::tType
    dt::tType
    uprev::uType
    tprev::tType
    alg::algType
    dtcache::tType
    tdir::Int
    saveiter::Int
    cache::CacheType
end

function (integrator::HistorySDEIntegrator)(t, deriv::Type=Val{0}; idxs=nothing)
    StochasticDiffEq.current_interpolant(t, integrator, idxs, deriv)
end

function (integrator::HistorySDEIntegrator)(
    val::AbstractArray, t::Union{Number,AbstractArray}, deriv::Type=Val{0}; idxs=nothing
)
    StochasticDiffEq.current_interpolant!(val, t, integrator, idxs, deriv)
end

mutable struct SDDEIntegrator{algType,IIP,uType,uEltype,tType,P,eigenType,tTypeNoUnits,uEltypeNoUnits,randType,randType2,rateType,solType,cacheType,F,G,F6,OType,noiseType,EventErrorType,CallbackCacheType,H,IType} <: AbstractSDDEIntegrator{algType,IIP,uType,tType}
    f::F
    g::G
    c::F6
    noise::noiseType
    uprev::uType
    tprev::tType
    # prev_idx::Int - from DDEIntegrator TODO
    # prev2_idx::Int - from DDEIntegrator TODO
    # fpsolver::FP - from DDEIntegrator TODO
    order_discontinuity_t0::Int  #- from DDEIntegrator TODO
    tracked_discontinuities::Vector{Discontinuity{tType,Rational{Int}}}
    # discontinuity_interp_points::Int #- from DDEIntegrator TODO
    # discontinuity_abstol::dAbsType #- from DDEIntegrator TODO
    # discontinuity_reltol::dRelType #- from DDEIntegrator TODO
    t::tType
    u::uType
    p::P
    dt::tType
    dtnew::tType
    dtpropose::tType
    dtcache::tType
    T::tType
    tdir::Int
    just_hit_tstop::Bool
    do_error_check::Bool
    isout::Bool
    event_last_time::Int
    vector_event_last_time::Int
    last_event_error::EventErrorType
    accept_step::Bool
    last_stepfail::Bool
    force_stepfail::Bool
    dtchangeable::Bool
    u_modified::Bool
    saveiter::Int
    alg::algType
    sol::solType
    cache::cacheType
    callback_cache::CallbackCacheType
    sqdt::tType
    W::randType
    P::randType2
    opts::OType
    iter::Int
    success_iter::Int
    eigen_est::eigenType
    EEst::tTypeNoUnits
    q::tTypeNoUnits
    qold::tTypeNoUnits
    q11::tTypeNoUnits
    history::H
    destats::DiffEqBase.DEStats
    integrator::IType # history integrator
end

function (integrator::SDDEIntegrator)(t, deriv::Type=Val{0}; idxs=nothing)
    StochasticDiffEq.current_interpolant(t, integrator, idxs, deriv)
end

function (integrator::SDDEIntegrator)(
    val::AbstractArray, t::Union{Number,AbstractArray}, deriv::Type=Val{0}; idxs=nothing
)
    StochasticDiffEq.current_interpolant!(val, t, integrator, idxs, deriv)
end


text:
using Test
using RigidBodyTreeInspector
using RigidBodyDynamics
using DrakeVisualizer
using CoordinateTransformations
using ColorTypes
using NBInclude
using Random

@testset "chain mechanism" begin
    mechanism = rand_chain_mechanism(Float64, [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1:5]]...)
    vis = Visualizer(mechanism; show_inertias=true);
    # We can draw the mechanism at a single state:
    state = MechanismState{Float64}(mechanism)
    rand!(state)
    settransform!(vis, state)
end

@testset "attach mechanism" begin
    mechanism = rand_chain_mechanism(Float64, [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1:5]]...)
    urdf = joinpath(dirname(@__FILE__), "..", "examples", "urdf", "Acrobot.urdf")
    mechanism2 = parse_urdf(Float64, urdf)
    attach!(mechanism, root_body(mechanism), mechanism2)
    vis = Visualizer(mechanism)
end

@testset "simulation and animation" begin
    mechanism = rand_chain_mechanism(Float64, [Revolute{Float64} for i = 1:10]...)
    vis = Visualizer(mechanism, show_inertias=true);
    state = MechanismState{Float64}(mechanism)
    zero!(state)
    settransform!(vis, state)
    times, states = simulate(state, 1, Δt = 0.001);
    animate(vis, mechanism, times, states)
end

@testset "urdf mechanism" begin
    urdf = joinpath(dirname(@__FILE__), "..", "examples", "urdf", "Acrobot.urdf")
    mechanism = parse_urdf(Float64, urdf)
    vis = Visualizer(mechanism, URDFVisuals(urdf))
    remove_fixed_tree_joints!(mechanism)
    vis = Visualizer(mechanism, URDFVisuals(urdf))
    mechanism = submechanism(mechanism, bodies(mechanism)[2])
    vis = Visualizer(mechanism, URDFVisuals(urdf))
end

@testset "geometry additions" begin
    urdf = joinpath(dirname(@__FILE__), "..", "examples", "urdf", "Acrobot.urdf")
    mechanism = parse_urdf(Float64, urdf)
    vis = Visualizer(mechanism, URDFVisuals(urdf))
    addgeometry!(vis, mechanism, default_frame(bodies(mechanism)[2]))
    addgeometry!(vis, mechanism, Point3D(default_frame(bodies(mechanism)[2]), 0.1, 0.1, 0.1))
    addgeometry!(vis, mechanism, default_frame(bodies(mechanism)[3]), HyperSphere(Point(0., 0, 0), 0.05))
    addgeometry!(vis, mechanism, default_frame(bodies(mechanism)[3]), GeometryData(HyperSphere(Point(0., 0, 0), 0.05), RGBA(0, 1, 0, 0.5), Translation(0.1, 0, 0)))
end

@testset "notebooks" begin
    @nbinclude(joinpath(@__DIR__, "..", "examples", "demo.ipynb"))
    @nbinclude(joinpath(@__DIR__, "..", "examples", "urdf.ipynb"))
end


text:
"""
`AbstractMesh{N,V<:AbstractPoint{N},F<:AbstractFace}`

Supertype for `N`-dimensional surface-mesh with faces of type `F` and `N`-dimensional vertices of type `V`

Each surface-mesh is collection of vertices and faces. Each face consists of vector of offsets into vertices vector.
"""
abstract type AbstractMesh{N,V <: AbstractPoint{N}, F <: AbstractPolygon} <: AbstractVector{F} end

Base.haslength(::AbstractMesh) = true
Base.eltype(::AbstractMesh{N,V,F}) where {N,V,F} = F
Base.IndexStyle(::Type{<:AbstractMesh}) = Base.IndexLinear()