text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    extract_eigenvalues(
        model,
        f::MOI.AbstractVectorFunction,
        d::Int,
        offset::Int,
        ::Type{T},
    ) where {T}

The vector `f` contains `t` (if `offset = 1`) or `(t, u)` (if `offset = 2`)
followed by the matrix `X` of dimension `d`.

This functions extracts the eigenvalues of `X` and returns a vector containing
`t` or `(t, u)`, a vector `MOI.VariableIndex` containing the eigenvalues of `X`,
the variables created and the index of the constraint created to extract the
eigenvalues.
"""
function extract_eigenvalues(
    model,
    f::MOI.AbstractVectorFunction,
    d::Int,
    offset::Int,
    ::Type{T},
) where {T}
    f_scalars = MOI.Utilities.eachscalar(f)
    tu = [f_scalars[i] for i in 1:offset]
    n = MOI.Utilities.trimap(d, d)
    X = f_scalars[offset.+(1:n)]
    N = MOI.Utilities.trimap(2d, 2d)
    Δ = MOI.add_variables(model, n)
    Z = [zero(MOI.ScalarAffineFunction{T}) for i in 1:(N-n)]
    for j in 1:d
        for i in j:d
            Z[MOI.Utilities.trimap(i, d + j)-n] = Δ[MOI.Utilities.trimap(i, j)]
        end
        Z[MOI.Utilities.trimap(d + j, d + j)-n] = Δ[MOI.Utilities.trimap(j, j)]
    end
    Y = MOIU.operate(vcat, T, X, MOI.Utilities.vectorize(Z))
    set = MOI.PositiveSemidefiniteConeTriangle(2d)
    sdindex = MOI.add_constraint(model, Y, set)
    D = Δ[MOI.Utilities.trimap.(1:d, 1:d)]
    return tu, D, Δ, sdindex
end

"""
    LogDetBridge{T,F,G,H,I}

The `LogDetConeTriangle` is representable by a
`PositiveSemidefiniteConeTriangle` and `ExponentialCone` constraints.

Indeed, ``\\log\\det(X) = \\log(\\delta_1) + \\cdots + \\log(\\delta_n)`` where
``\\delta_1``, ..., ``\\delta_n`` are the eigenvalues of ``X``.

Adapting the method from [1, p. 149], we see that ``t \\le u \\log(\\det(X/u))``
for ``u > 0`` if and only if there exists a lower triangular matrix ``Δ`` such
that
```math
\\begin{align*}
  \\begin{pmatrix}
    X & Δ\\\\
    Δ^\\top & \\mathrm{Diag}(Δ)
  \\end{pmatrix} & \\succeq 0\\\\
  t & \\le u \\log(Δ_{11}/u) + u \\log(Δ_{22}/u) + \\cdots + u \\log(Δ_{nn}/u)
\\end{align*}
```

[1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex
    optimization: analysis, algorithms, and engineering applications*. Society
    for Industrial and Applied Mathematics, 2001.
```
"""
struct LogDetBridge{T,F,G,H,I} <: AbstractBridge
    Δ::Vector{MOI.VariableIndex}
    l::Vector{MOI.VariableIndex}
    sdindex::MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle}
    lcindex::Vector{MOI.ConstraintIndex{G,MOI.ExponentialCone}}
    tlindex::MOI.ConstraintIndex{H,MOI.LessThan{T}}
end

const LogDet{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{LogDetBridge{T},OT}

function bridge_constraint(
    ::Type{LogDetBridge{T,F,G,H,I}},
    model,
    f::I,
    s::MOI.LogDetConeTriangle,
) where {T,F,G,H,I}
    d = s.side_dimension
    tu, D, Δ, sdindex = extract_eigenvalues(model, f, d, 2, T)
    t, u = tu
    l = MOI.add_variables(model, d)
    lcindex = [sublog(model, l[i], u, D[i], T) for i in eachindex(l)]
    tlindex = subsum(model, t, l, T)
    return LogDetBridge{T,F,G,H,I}(Δ, l, sdindex, lcindex, tlindex)
end

function MOI.supports_constraint(
    ::Type{LogDetBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.LogDetConeTriangle},
) where {T}
    return true
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:LogDetBridge})
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:LogDetBridge{T,F,G,H}},
) where {T,F,G,H}
    return Tuple{Type,Type}[
        (F, MOI.PositiveSemidefiniteConeTriangle),
        (G, MOI.ExponentialCone),
        (H, MOI.LessThan{T}),
    ]
end

function concrete_bridge_type(
    ::Type{<:LogDetBridge{T}},
    I::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.LogDetConeTriangle},
) where {T}
    S = MOIU.scalar_type(I)
    G = MOIU.promote_operation(vcat, T, S, MOI.VariableIndex)
    F = MOIU.promote_operation(vcat, T, G, T)
    H = MOIU.promote_operation(+, T, S, MOI.ScalarAffineFunction{T})
    return LogDetBridge{T,F,G,H,I}
end

"""
    sublog(
        model,
        x::MOI.VariableIndex,
        y::MOI.AbstractScalarFunction,
        z::MOI.VariableIndex,
        ::Type{T},
    ) where {T}

Constrains ``x \\le y \\log(z/y)`` and returns the constraint index.
"""
function sublog(
    model,
    x::MOI.VariableIndex,
    y::MOI.AbstractScalarFunction,
    z::MOI.VariableIndex,
    ::Type{T},
) where {T}
    return MOI.add_constraint(
        model,
        MOI.Utilities.operate(vcat, T, x, y, z),
        MOI.ExponentialCone(),
    )
end

"""
    subsum(
        model,
        t::MOI.AbstractScalarFunction,
        l::Vector{MOI.VariableIndex},
        ::Type{T}
    ) where {T}

Constrains ``t \\le l_1 + \\cdots + l_n`` where `n` is the length of `l` and
returns the constraint index.
"""
function subsum(
    model,
    t::MOI.AbstractScalarFunction,
    l::Vector{MOI.VariableIndex},
    ::Type{T},
) where {T}
    rhs = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(one(T), l), zero(T))
    f = MOI.Utilities.operate!(-, T, t, rhs)
    return MOI.Utilities.normalize_and_add_constraint(
        model,
        f,
        MOI.LessThan(zero(T)),
        allow_modify_function = true,
    )
end

function MOI.get(b::LogDetBridge, ::MOI.NumberOfVariables)::Int64
    return length(b.Δ) + length(b.l)
end

MOI.get(b::LogDetBridge, ::MOI.ListOfVariableIndices) = [b.Δ; b.l]

function MOI.get(
    ::LogDetBridge{T,F},
    ::MOI.NumberOfConstraints{F,MOI.PositiveSemidefiniteConeTriangle},
)::Int64 where {T,F}
    return 1
end

function MOI.get(
    b::LogDetBridge{T,F,G},
    ::MOI.NumberOfConstraints{G,MOI.ExponentialCone},
)::Int64 where {T,F,G}
    return length(b.lcindex)
end

function MOI.get(
    ::LogDetBridge{T,F,G,H},
    ::MOI.NumberOfConstraints{H,MOI.LessThan{T}},
)::Int64 where {T,F,G,H}
    return 1
end

function MOI.get(
    b::LogDetBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,MOI.PositiveSemidefiniteConeTriangle},
) where {T,F}
    return [b.sdindex]
end

function MOI.get(
    b::LogDetBridge{T,F,G},
    ::MOI.ListOfConstraintIndices{G,MOI.ExponentialCone},
) where {T,F,G}
    return copy(b.lcindex)
end

function MOI.get(
    b::LogDetBridge{T,F,G,H},
    ::MOI.ListOfConstraintIndices{H,MOI.LessThan{T}},
) where {T,F,G,H}
    return [b.tlindex]
end

function MOI.delete(model::MOI.ModelLike, bridge::LogDetBridge)
    MOI.delete(model, bridge.tlindex)
    MOI.delete(model, bridge.lcindex)
    MOI.delete(model, bridge.sdindex)
    MOI.delete(model, bridge.l)
    MOI.delete(model, bridge.Δ)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::LogDetBridge,
)
    d = length(bridge.lcindex)
    Δ = MOI.get(model, MOI.VariablePrimal(), bridge.Δ)
    t =
        MOI.get(model, attr, bridge.tlindex) +
        sum(MOI.get(model, attr, ci)[1] for ci in bridge.lcindex)
    u = MOI.get(model, attr, first(bridge.lcindex))[2]
    x = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t, u, x)
end

# [X Δ; Δ' Diag(Δ)] in PSD
# t - sum(l) >= 0
# (l_i, u, Δ_ii) in Exp
# (t, u, x) in LogDet <=>
# exists Δ, l such that At + Bu + Cx + DΔ + El in (PSD, >=, Exp_i)
# so LogDet* = [A'; B'; C'] (PSD, >=, Exp_i)*
# and 0 = [D'; E'] (PSD, >=, Exp_i)*
# where
# A = [0, 0, 0, 1, 0, 0, 0]
# B = [0, 0, 0, 0, 0, 1, 0]
# C = [I, 0, 0, 0, 0, 0, 0]
# D = [0, I, I(i=j), 0, 0, 0, I(i=j)]
# E = [0, 0, 0, 1, I, 0, 0]
# so given dual q = (a, b, c, d, e, f, g), we get
# t = A' q = d => t = d
# u = B' q = f => u = sum(f)
# x = C' q = a => x = a
# offdiag(b) = 0
# 0 = D' q = diag(b) + diag(c) + g => g = -diag(b) - diag(c)
# let b = 0, so g_i = -c_i
# 0 = E' q = d + e => d = -e
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::LogDetBridge,
)
    t_dual = MOI.get(model, attr, bridge.tlindex)
    u_dual =
        sum(MOI.get(model, attr, lcindex_i)[2] for lcindex_i in bridge.lcindex)
    x_dual = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t_dual, u_dual, x_dual)
end

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::LogDetBridge)
    return MOI.LogDetConeTriangle(length(bridge.l))
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::LogDetBridge{T,F,G,H,I},
) where {T,F,G,H,I}
    n = length(bridge.Δ)
    Y = MOI.get(model, attr, bridge.sdindex)::F
    X = MOI.Utilities.convert_approx(I, MOI.Utilities.eachscalar(Y)[1:n])
    slog = MOI.get(model, attr, first(bridge.lcindex))::G
    u = MOI.Utilities.eachscalar(slog)[2]
    ssu = MOI.get(model, attr, bridge.tlindex)::H
    taff = MOI.Utilities.remove_variable(ssu, bridge.l)
    SI = MOI.Utilities.scalar_type(I)
    t = MOI.Utilities.convert_approx(SI, taff)
    return MOI.Utilities.operate(vcat, T, t, u, X)::I
end

"""
    RootDetBridge{T,F,G,H}

The `RootDetConeTriangle` is representable by a
`PositiveSemidefiniteConeTriangle` and an `GeometricMeanCone` constraints; see
[1, p. 149].

Indeed, ``t \\le \\det(X)^{1/n}`` if and only if there exists a lower triangular
matrix ``Δ`` such that:
```math
\\begin{align*}
  \\begin{pmatrix}
    X & Δ\\\\
    Δ^\\top & \\mathrm{Diag}(Δ)
  \\end{pmatrix} & \\succeq 0\\\\
  t & \\le (Δ_{11} Δ_{22} \\cdots Δ_{nn})^{1/n}
\\end{align*}
```

[1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex
    optimization: analysis, algorithms, and engineering applications*. Society
    for Industrial and Applied Mathematics, 2001.
"""
struct RootDetBridge{T,F,G,H} <: AbstractBridge
    Δ::Vector{MOI.VariableIndex}
    sdindex::MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle}
    gmindex::MOI.ConstraintIndex{G,MOI.GeometricMeanCone}
end

const RootDet{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{RootDetBridge{T},OT}

function bridge_constraint(
    ::Type{RootDetBridge{T,F,G,H}},
    model,
    f::G,
    s::MOI.RootDetConeTriangle,
) where {T,F,G,H}
    d = s.side_dimension
    tu, D, Δ, sdindex = extract_eigenvalues(model, f, d, 1, T)
    t = tu[1]
    DF = MOI.VectorOfVariables(D)
    gmindex = MOI.add_constraint(
        model,
        MOI.Utilities.operate(vcat, T, t, DF),
        MOI.GeometricMeanCone(d + 1),
    )
    return RootDetBridge{T,F,G,H}(Δ, sdindex, gmindex)
end

function MOI.supports_constraint(
    ::Type{<:RootDetBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.RootDetConeTriangle},
) where {T}
    return true
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:RootDetBridge})
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:RootDetBridge{T,F,G}},
) where {T,F,G}
    return Tuple{Type,Type}[
        (F, MOI.PositiveSemidefiniteConeTriangle),
        (G, MOI.GeometricMeanCone),
    ]
end

function concrete_bridge_type(
    ::Type{<:RootDetBridge{T}},
    H::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.RootDetConeTriangle},
) where {T}
    S = MOIU.scalar_type(H)
    G = MOIU.promote_operation(vcat, T, S, MOI.VariableIndex)
    F = MOIU.promote_operation(vcat, T, G, T)
    return RootDetBridge{T,F,G,H}
end

MOI.get(b::RootDetBridge, ::MOI.NumberOfVariables)::Int64 = length(b.Δ)

MOI.get(b::RootDetBridge, ::MOI.ListOfVariableIndices) = copy(b.Δ)

function MOI.get(
    ::RootDetBridge{T,F},
    ::MOI.NumberOfConstraints{F,MOI.PositiveSemidefiniteConeTriangle},
)::Int64 where {T,F}
    return 1
end

function MOI.get(
    ::RootDetBridge{T,F,G},
    ::MOI.NumberOfConstraints{G,MOI.GeometricMeanCone},
)::Int64 where {T,F,G}
    return 1
end

function MOI.get(
    b::RootDetBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,MOI.PositiveSemidefiniteConeTriangle},
) where {T,F}
    return [b.sdindex]
end

function MOI.get(
    b::RootDetBridge{T,F,G},
    ::MOI.ListOfConstraintIndices{G,MOI.GeometricMeanCone},
) where {T,F,G}
    return [b.gmindex]
end

function MOI.delete(model::MOI.ModelLike, bridge::RootDetBridge)
    MOI.delete(model, bridge.gmindex)
    MOI.delete(model, bridge.sdindex)
    MOI.delete(model, bridge.Δ)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::RootDetBridge,
)
    t = MOI.get(model, attr, bridge.gmindex)[1]
    x = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t, x)
end

# (t, x) in RootDet <=> exists Δ such that At + Bx + CΔ in (PSD, GeoMean)
# so RootDet* = [A'; B'] (PSD, GeoMean)*
# and 0 = [C'] (PSD, GeoMean)*
# where
# A = [0, 0, 0, 1, 0]
# B = [I, 0, 0, 0, 0]
# C = [0, I, I(i=j), 0, I(i=j)]
# so given dual q = (a, b, c, d, e), we get
# t = A' q = d => t = d
# x = B' q = a => x = a
# offdiag(b) = 0
# 0 = C' q = diag(b) + diag(c) + e => e = -diag(b) - diag(c)
# let b = 0, so e_i = -c_i
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual},
    bridge::RootDetBridge,
)
    t_dual = MOI.get(model, attr, bridge.gmindex)[1]
    x_dual = MOI.get(model, attr, bridge.sdindex)[1:length(bridge.Δ)]
    return vcat(t_dual, x_dual)
end

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::RootDetBridge)
    n = length(bridge.Δ)
    d = MOI.Utilities.side_dimension_for_vectorized_dimension(n)
    return MOI.RootDetConeTriangle(d)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::RootDetBridge{T,F,G,H},
) where {T,F,G,H}
    n = length(bridge.Δ)
    Y = MOI.get(model, attr, bridge.sdindex)::F
    X = MOI.Utilities.convert_approx(H, MOI.Utilities.eachscalar(Y)[1:n])
    func = MOI.get(model, attr, bridge.gmindex)::G
    t = MOI.Utilities.eachscalar(func)[1]
    return MOI.Utilities.operate(vcat, T, t, X)::H
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using Documenter, ReadableRegex

makedocs(;
    modules=[ReadableRegex],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jkrumbiegel/ReadableRegex.jl/blob/{commit}{path}#L{line}",
    sitename="ReadableRegex.jl",
    authors="Julius Krumbiegel",
    assets=String[],
)

deploydocs(;
    repo="github.com/jkrumbiegel/ReadableRegex.jl",
)

ext: jl

max_stars_repo_name: jkrumbiegel/ReadableRegex

max_issues_repo_name: jkrumbiegel/ReadableRegex

max_forks_repo_name: jkrumbiegel/ReadableRegex


text:
using Distributions
using IndependentRandomSequences
using OrderStatistics
using Combinatorics
using Base.Test

Γ = gamma

#
#Tests for Uniform Case
#
for N = 3:5
  U = Uniform()
  sequence = IIDRandomSequence(U,N)

  #scalar order statistics
  Vs = orderstatistics(sequence)
  @test map(mean,Vs) == collect(1:N)/(N+1)

  #joint order statistics
  fulljoint = sort(sequence)
  @test jointorderstatistic(fulljoint) == fulljoint
  @test_approx_eq_eps(pdf(fulljoint,collect(1:N)/(N+1)),Γ(N+1),10^3*eps())
  for n=2:N
    for orders=combinations(N,n)
      partialjoint = jointorderstatistic(sequence,orders)
      @test mean(partialjoint ) == orders/(N+1)
    end
  end

  #range
  @test mean(range(sequence)) == (N-1)/(N+1)

  #spacing
  @test map(mean,spacings(sequence)) == ones(N+1)/(N+1)

  #jointspacing
  @test mean(jointspacing(sequence)) == ones(N+1)/(N+1)

  #idempotency property

end

#
# Qualitative ordering tests
#
N = 5
for X = [Uniform(),Normal(),Exponential()]
  sequence = IIDRandomSequence(X,N)

  #scalar order statistics
  Ys = orderstatistics(sequence)
  @test issorted(map(mean,Ys))
  @test mean(Ys[1]) < mean(X)
  @test mean(Ys[end]) > mean(X)

  #joint order statistics
  jYs = jointorderstatistics(sequence)
  for jY=jYs
    @test issorted(rand(jY))
  end

end


text:
struct FEMProblem
    s # system
    coords
    connec
    neumdofs
    fixeddofs
    b # Forcing term
    Fnotzerocoord
end

_homogeneize(p::FEMProblem) = _homogeneize(p.s.M, p.s.K, p.b)

# nodes of interest:
#    nodeidx_hor = 3
#    nodeidx_ver = 4
function load_lambs(; mesh=1,
                      file_matrices=joinpath(@__DIR__, "mesh$mesh", "matrices.mat"),
                      file_mesh=joinpath(@__DIR__, "mesh$mesh", "mesh.mat"),
                      file_loads=joinpath(@__DIR__, "mesh$mesh", "loads.mat"))

    # load matrices
    d = matread(file_matrices)
    K = d["KTred"]
    M = d["massMatred"]
    @assert isdiag(M)
    M = Diagonal(M) # it is diagonal
    n = size(K, 1)
    C = spzeros(n, n)
    sys = SecondOrderLinearContinuousSystem(M, C, K)

    neumdofs = d["neumdofs"]
    neumdofs = Int.(neumdofs)[:]

    # load mesh
    mesh = matread(file_mesh)
    coords = mesh["mesh"]["nodesCoords"]
    connec = mesh["mesh"]["conecCell"]
    connec = Int.(connec)

    # fixed degrees of freedom
    fixeddofs = setdiff(1:6*size(coords, 1), neumdofs)

    # load external force definition
    df = matread(file_loads)
    F = df["factorLoadsFextCell"][1, 3][:]
    F = F[neumdofs]
    Fnotzerocoord = findall(!iszero, F)

    return FEMProblem(sys, coords, connec, neumdofs, fixeddofs, F, Fnotzerocoord)
end

# ======================
# Funciones de carga
# ======================

Heaviside(t) = t >= 0 ? 1 : 0

# funcion auxiliar que devuelve el factor de carga aplicado en la coordenada coord
function factor_de_carga(m, NSTEPS; T=1.0, coord=5, func=Heaviside)
    tiempos = range(0.0, T, length=NSTEPS)

    F = Vector{Vector{Float64}}(undef, NSTEPS)
    for i in 1:NSTEPS
        aux = zeros(m)
        aux[coord] .= func(tiempos[i])

        F[i] = aux
    end
    return F
end

function step_size_lamb(p, δtfac)
    nElem = sum(p.connec[:, 2] .== 3) # number of triangular elements
    CFL = 0.125
    cP  = 3200
    h   = sqrt(2*3200^2 / nElem) * δtfac
    δt = CFL * h / cP
end

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

include("../common.jl")
name = "oneAPI_Level_Zero_Headers"

# Bash recipe for building across all platforms
script = raw"""
cd level-zero
install_license LICENSE

mkdir -p ${includedir}/level_zero
mv include/*.h ${includedir}/level_zero
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("include/level_zero/ze_api.h", :ze_api)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

build_tarballs(ARGS, name, api_version, sources, script, platforms, products, dependencies)

# bump counter: 1

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# It is suggested that test code for MLJ.jl include files be placed in
# a file of the same name under "test/" (and included below) and that
# this test code be wrapped in a module. Any new module name will do -
# eg, `module TestDatasets` for code testing `datasets.jl`.

using MLJ
using MLJBase
using Test
using Random

@constant junk=ConstantRegressor()

@testset "utilities" begin
  @test include("utilities.jl")
end

@testset "parameters" begin
  @test include("parameters.jl")
end

@testset "Machines" begin
  @test include("machines.jl")
end

@testset "networks" begin
  @test include("networks.jl")
end

if VERSION ≥ v"1.3.0-"
  @testset "arrows" begin
    @test include("arrows.jl")
  end
end

@testset "composites" begin
  @test include("composites.jl")
end

@testset "pipelines" begin
  @test include("pipelines.jl")
end

@testset "resampling" begin
    @test include("resampling.jl")
end

@testset "tuning" begin
    @test include("tuning.jl")
end

@testset "ensembles" begin
    @test include("ensembles.jl")
end

@testset "matching models to data" begin
    @test include("model_matching.jl")
end

@testset "tasks" begin
  @test include("tasks.jl")
end

@testset "scitypes" begin
    @test include("scitypes.jl")
end

# @testset "registry environment compatibility" begin
#     @test MLJ.Registry.@update(true)
# end

## TEST THE EXAMPLES

const exdir = joinpath(MLJ.srcdir, "..", "examples")

# uncomment remaining code to test examples:
# @testset "/examples" begin
#     @test include(joinpath(exdir, "using_tasks.jl"))
#     @test include(joinpath(exdir, "random_forest.jl"))
#     @test include(joinpath(exdir, "two_parameter_tune.jl"))
# end


text:
using SimulacoesPedro
using Test

@testset "SimulacoesPedro.jl" begin
    data = Data()
    p = initial_point(data)
    @test utotal(p,data) ≈ utotal2(p,data)
end


text:
include("01_one_site_state_basics.jl")

println("
#######################################################
# Tutorial 2                                          #
#                                                     #
# 1-site states                                       #
#######################################################
")

# "S=1/2" defines a state and operator basis:
# a set of predefined states and operators.
# Additionally: "Qubit", "Qudit", "Fermion", "Electron", etc.
i = Index(2, "S=1/2")

# Same as before
Zp = ITensor([1 0], i)

# Predefined (extensible) state definitions
Zp = state("Z+", i)

# Make a "down" state (-1 eigenvalue of `Z`),
# denoted as `|Z-⟩`.
Zm = ITensor([0 1], i)

# Predefined definition
Zm = state("Z-", i)

@visualize Zm

# Some more complicated states:
# |X+⟩ = (|Z+⟩ + |Z-⟩) / √2

# Define through setting elements
Xp = ITensor(i)
Xp[i => 1] = 1 / √2
Xp[i => 2] = 1 / √2

# Define in terms of |Z+⟩, |Z-⟩
Xp = (Zp + Zm) / √2

# Use predefined definition
Xp = state("X+", i)

# |X-⟩ = (|Z+⟩ - |Z-⟩) / √2
Xm = (Zp - Zm) / √2

# Or:
Xm = state("X-", i)

# Perform some inner products:
@show dag(Zp) * Zp

# Convert from order-0 tensor to Number:
@show (dag(Zp) * Zp)[]

# Shorthand:
#
# inner(x, y) = (dag(x) * y)[]
#
# norm(x) = √(inner(x, x))
@show inner(Zp, Zp)
@show norm(Zp)
@show inner(Zp, Zp)
@show inner(Zp, Xm)
@show inner(Zm, Xm)

@visualize dag(Zp) * Xm


text:
module julia

export topplebyforce, topplebyboundary

function topplebyforce(dominoes::String)::String
    force_on_idx = fill(0, length(dominoes))

    # Calculate the 'leftward' force on each domino by iterating backwards over the
    # string, setting the current force to an arbitrarily large number when an 'L' is
    # encountered, to 0 when an 'R' is encountered (stops the momentum), and to one
    # less than the previous force when an '.' is encountered. This simulates the force
    # falling off over distance from the 'L' domino that started the chain reaction.
    current_force = 0
    for idx in length(dominoes):-1:1
        if dominoes[idx] == 'R'; current_force = 0; end
        if dominoes[idx] == 'L'; current_force = typemax(Int64); end
        force_on_idx[idx] = current_force
        current_force -= current_force > 0 ? 1 : 0
    end

    # Now calculate the 'rightward' force in much the same was as above, except this
    # time, we determine which character to place into the output string by comparing
    # the 'leftward' force to the 'rightward' force at each position. If the 'leftward'
    # force is greater, the domino falls left. If the 'rightward' force is greater, it
    # falls right. If they're equal, they cancel out and the domino remains standing.
    current_force = 0
    result_string = ""
    for idx in 1:length(dominoes)
        if dominoes[idx] == 'R'; current_force = typemax(Int64); end
        if dominoes[idx] == 'L'; current_force = 0; end
        if force_on_idx[idx] == current_force; result_string *= '.'; end
        if force_on_idx[idx] > current_force;  result_string *= 'L'; end
        if force_on_idx[idx] < current_force;  result_string *= 'R'; end
        current_force -= current_force > 0 ? 1 : 0
    end

    return result_string
end

function topplebyboundary(dominoes::String)::String
    # Build an array of tuples, where each tuple contains the index for the first and
    # last element of the current 'region', where each region is a sequence of
    # characters with either an 'R' or an 'L' at either end.
    boundaries = []
    (low, high) = (1, 1)

    for (idx, tile) in enumerate(dominoes)
        high = copy(idx)
        if tile != '.'
            push!(boundaries, (low, high))
            low = copy(high)
        end
    end

    if last(boundaries)[2] != length(dominoes)
        push!(boundaries, (low, high))
    end

    # Armed with the knowledge of where the different 'regions' lie, we can determine
    # what happens at each index based on the following criteria:
    # - If the region begins on the left and goes to an 'L', all dominoes will be 'L'
    # - If the region starts with 'R' an goes to the edge, all dominoes will be 'R'
    # - If the region starts and ends with the same value, all dominoes will be that value
    # - If the region starts with 'R' and ends with 'L', all dominoes will fall towards
    #   the center of the region
    # By replacing characters in the original string, characters that are unaffected,
    # such as those in a region that starts with 'L' and ends with 'R', are retained
    # in the final output.
    dominoes = collect(dominoes)
    for (first, last) in boundaries
        if dominoes[first] == '.' && dominoes[last] == 'L'
            dominoes[first:last] .= 'L'
        end

        if dominoes[first] == 'R' && dominoes[last] == '.'
            dominoes[first:last] .= 'R'
        end

        if dominoes[first] == dominoes[last]
            dominoes[first:last] .= dominoes[first]
        end

        if dominoes[first] == 'R' && dominoes[last] == 'L'
            while first + 1 < last - 1
                dominoes[first + 1] = dominoes[first]
                dominoes[last - 1] = dominoes[last]
                first += 1; last -= 1
            end
        end
    end

    return join(dominoes)
end

end # module

max_stars_repo_name: JenniferToops/CodingDojo

max_issues_repo_name: JenniferToops/CodingDojo

max_forks_repo_name: JenniferToops/CodingDojo


text:
using Revise
Revise.async_steal_repl_backend()
map(f-> Revise.track(Jus, f),
    ["types.jl", "protocol.jl", "vars.jl", "commands.jl", "Jus.jl", "server.jl", "example1.jl"])
includet("load.jl")


text:
using JGE
JGE.main()


text: multiply = (a, b) -> a * b

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function random(numbits)
  numbytes = fld(numbits, 8)
  if numbits % 8 != 0
    warn("The requested number of random bits is not divisible by 8. Rounding to the nearest multiple of 8.")
    numbytes += 1
  end

  ret = zeros(Uint8, numbytes)
  ccall((:RAND_bytes, "libcrypto"), Void, (Ptr{Uint8}, Int), ret, numbytes)

  return ret
end


text:
module GaPLAC

export invnormaltransform

# GP Packages
using AbstractGPs
using KernelFunctions
using GPLikelihoods
using Distances: Euclidean

# Bayesian stuff
using Turing
using Distributions
using Turing: Variational

# Data stuff
using CSV
using DataFrames
using CairoMakie
using StatsBase

# Other stuff
using ReTest

@testset "Module" begin
    @test true
end

include("gp_parts.jl")
include("liklihoods.jl")
include("abstractgp_translations.jl")
include("utils.jl")
include("interface.jl")

end # module

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

@doc raw"""
    CompressibleEulerEquations1D(gamma)

The compressible Euler equations for an ideal gas with ratio of specific heats `gamma`
in one space dimension.
"""
struct CompressibleEulerEquations1D{RealT<:Real} <: AbstractCompressibleEulerEquations{1, 3}
  gamma::RealT               # ratio of specific heats
  inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications

  function CompressibleEulerEquations1D(gamma)
    γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1))
    new{typeof(γ)}(γ, inv_gamma_minus_one)
  end
end

varnames(::typeof(cons2cons), ::CompressibleEulerEquations1D) = ("rho", "rho_v1", "rho_e")
varnames(::typeof(cons2prim), ::CompressibleEulerEquations1D) = ("rho", "v1", "p")

"""
    initial_condition_constant(x, t, equations::CompressibleEulerEquations1D)

A constant initial condition to test free-stream preservation.
"""
function initial_condition_constant(x, t, equations::CompressibleEulerEquations1D)
  rho = 1.0
  rho_v1 = 0.1
  rho_e = 10.0
  return SVector(rho, rho_v1, rho_e)
end

"""
    initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations1D)

A smooth initial condition used for convergence tests in combination with
[`source_terms_convergence_test`](@ref)
(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
"""
function initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations1D)
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  ini = c + A * sin(ω * (x[1] - t))

  rho = ini
  rho_v1 = ini
  rho_e = ini^2

  return SVector(rho, rho_v1, rho_e)
end

"""
    source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquations1D)

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref)
(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
"""
@inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquations1D)
  # Same settings as in `initial_condition`
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  γ = equations.gamma

  x1, = x

  si, co = sincos(ω * (x1 - t))
  rho = c + A * si
  rho_x = ω * A * co

  # Note that d/dt rho = -d/dx rho.
  # This yields du2 = du3 = d/dx p (derivative of pressure).
  # Other terms vanish because of v = 1.
  du1 = zero(eltype(u))
  du2 = rho_x * (2 * rho - 0.5) * (γ - 1)
  du3 = du2

  return SVector(du1, du2, du3)
end

"""
    initial_condition_density_wave(x, t, equations::CompressibleEulerEquations1D)

A sine wave in the density with constant velocity and pressure; reduces the
compressible Euler equations to the linear advection equations.
This setup is the test case for stability of EC fluxes from paper
- Gregor J. Gassner, Magnus Svärd, Florian J. Hindenlang (2020)
  Stability issues of entropy-stable and/or split-form high-order schemes
  [arXiv: 2007.09026](https://arxiv.org/abs/2007.09026)
with the following parameters
- domain [-1, 1]
- mesh = 4x4
- polydeg = 5
"""
function initial_condition_density_wave(x, t, equations::CompressibleEulerEquations1D)
  v1 = 0.1
  rho = 1 + 0.98 * sinpi(2 * (x[1] - t * v1))
  rho_v1 = rho * v1
  p = 20
  rho_e = p / (equations.gamma - 1) + 1/2 * rho * v1^2
  return SVector(rho, rho_v1, rho_e)
end

"""
    initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations1D)

A weak blast wave taken from
- Sebastian Hennemann, Gregor J. Gassner (2020)
  A provably entropy stable subcell shock capturing approach for high order split form DG
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations1D)
  # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Set up polar coordinates
  inicenter = SVector(0.0)
  x_norm = x[1] - inicenter[1]
  r = abs(x_norm)
  # The following code is equivalent to
  # phi = atan(0.0, x_norm)
  # cos_phi = cos(phi)
  # in 1D but faster
  cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)

  # Calculate primitive variables
  rho = r > 0.5 ? 1.0 : 1.1691
  v1  = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  p   = r > 0.5 ? 1.0 : 1.245

  return prim2cons(SVector(rho, v1, p), equations)
end

"""
    initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations1D)

One dimensional variant of the setup used for convergence tests of the Euler equations
with self-gravity from
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
!!! note
    There is no additional source term necessary for the manufactured solution in one
    spatial dimension. Thus, [`source_terms_eoc_test_coupled_euler_gravity`](@ref) is not
    present there.
"""
function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations1D)
  # OBS! this assumes that γ = 2 other manufactured source terms are incorrect
  if equations.gamma != 2.0
    error("adiabatic constant must be 2 for the coupling convergence test")
  end
  c = 2.0
  A = 0.1
  ini = c + A * sinpi(x[1] - t)
  G = 1.0 # gravitational constant

  rho = ini
  v1 = 1.0
  p = 2 * ini^2 * G / pi # * 2 / ndims, but ndims==1 here

  return prim2cons(SVector(rho, v1, p), equations)
end

# Calculate 1D flux for a single point
@inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations1D)
  rho, rho_v1, rho_e = u
  v1 = rho_v1 / rho
  p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)
  # Ignore orientation since it is always "1" in 1D
  f1 = rho_v1
  f2 = rho_v1 * v1 + p
  f3 = (rho_e + p) * v1
  return SVector(f1, f2, f3)
end

"""
    flux_shima_etal(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

This flux is is a modification of the original kinetic energy preserving two-point flux by
- Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018)
  Kinetic energy and entropy preserving schemes for compressible flows
  by split convective forms
  [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058)

The modification is in the energy flux to guarantee pressure equilibrium and was developed by
- Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020)
  Preventing spurious pressure oscillations in split convective form discretizations for
  compressible flows
  [DOI: 10.1016/j.jcp.2020.110060](https://doi.org/10.1016/j.jcp.2020.110060)
"""
@inline function flux_shima_etal(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
  rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg  = 1/2 * ( v1_ll +  v1_rr)
  p_avg   = 1/2 * (  p_ll +   p_rr)
  kin_avg = 1/2 * (v1_ll * v1_rr)

  # Calculate fluxes
  # Ignore orientation since it is always "1" in 1D
  pv1_avg = 1/2 * (p_ll*v1_rr + p_rr*v1_ll)
  f1 = rho_avg * v1_avg
  f2 = rho_avg * v1_avg * v1_avg + p_avg
  f3 = p_avg*v1_avg * equations.inv_gamma_minus_one + rho_avg*v1_avg*kin_avg + pv1_avg

  return SVector(f1, f2, f3)
end

"""
    flux_kennedy_gruber(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

Kinetic energy preserving two-point flux by
- Kennedy and Gruber (2008)
  Reduced aliasing formulations of the convective terms within the
  Navier-Stokes equations for a compressible fluid
  [DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020)
"""
@inline function flux_kennedy_gruber(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_e_ll = last(u_ll)
  rho_e_rr = last(u_rr)
  rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
  rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg  = 1/2 * ( v1_ll +  v1_rr)
  p_avg   = 1/2 * (  p_ll +   p_rr)
  e_avg   = 1/2 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr)

  # Ignore orientation since it is always "1" in 1D
  f1 = rho_avg * v1_avg
  f2 = rho_avg * v1_avg * v1_avg + p_avg
  f3 = (rho_avg * e_avg + p_avg) * v1_avg

  return SVector(f1, f2, f3)
end

"""
    flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

Entropy conserving two-point flux by
- Chandrashekar (2013)
  Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes
  for Compressible Euler and Navier-Stokes Equations
  [DOI: 10.4208/cicp.170712.010313a](https://doi.org/10.4208/cicp.170712.010313a)
"""
@inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
  rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)
  beta_ll = 0.5 * rho_ll / p_ll
  beta_rr = 0.5 * rho_rr / p_rr
  specific_kin_ll = 0.5 * (v1_ll^2)
  specific_kin_rr = 0.5 * (v1_rr^2)

  # Compute the necessary mean values
  rho_avg = 0.5 * (rho_ll + rho_rr)
  rho_mean  = ln_mean(rho_ll, rho_rr)
  beta_mean = ln_mean(beta_ll, beta_rr)
  beta_avg = 0.5 * (beta_ll + beta_rr)
  v1_avg = 0.5 * (v1_ll + v1_rr)
  p_mean = 0.5 * rho_avg / beta_avg
  velocity_square_avg = specific_kin_ll + specific_kin_rr

  # Calculate fluxes
  # Ignore orientation since it is always "1" in 1D
  f1 = rho_mean * v1_avg
  f2 = f1 * v1_avg + p_mean
  f3 = f1 * 0.5*(1/(equations.gamma-1)/beta_mean - velocity_square_avg)+f2*v1_avg

  return SVector(f1, f2, f3)
end

"""
    flux_ranocha(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

Entropy conserving and kinetic energy preserving two-point flux by
- Hendrik Ranocha (2018)
  Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods
  for Hyperbolic Balance Laws
  [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743)
See also
- Hendrik Ranocha (2020)
  Entropy Conserving and Kinetic Energy Preserving Numerical Methods for
  the Euler Equations Using Summation-by-Parts Operators
  [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42)
"""
@inline function flux_ranocha(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
  rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)

  # Compute the necessary mean values
  rho_mean = ln_mean(rho_ll, rho_rr)
  # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
  # in exact arithmetic since
  #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
  #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
  inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
  v1_avg = 0.5 * (v1_ll + v1_rr)
  p_avg  = 0.5 * (p_ll + p_rr)
  velocity_square_avg = 0.5 * (v1_ll*v1_rr)

  # Calculate fluxes
  # Ignore orientation since it is always "1" in 1D
  f1 = rho_mean * v1_avg
  f2 = f1 * v1_avg + p_avg
  f3 = f1 * ( velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one ) + 0.5 * (p_ll*v1_rr + p_rr*v1_ll)

  return SVector(f1, f2, f3)
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
# maximum velocity magnitude plus the maximum speed of sound
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations1D)
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  # Calculate primitive variables and speed of sound
  v1_ll = rho_v1_ll / rho_ll
  v_mag_ll = abs(v1_ll)
  p_ll = (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * v_mag_ll^2)
  c_ll = sqrt(equations.gamma * p_ll / rho_ll)
  v1_rr = rho_v1_rr / rho_rr
  v_mag_rr = abs(v1_rr)
  p_rr = (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * v_mag_rr^2)
  c_rr = sqrt(equations.gamma * p_rr / rho_rr)

  λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
end

# Calculate minimum and maximum wave speeds for HLL-type fluxes
@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations1D)
  rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
  rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)

  λ_min = v1_ll - sqrt(equations.gamma * p_ll / rho_ll)
  λ_max = v1_rr + sqrt(equations.gamma * p_rr / rho_rr)

  return λ_min, λ_max
end

"""
    flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro
[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf)
Signal speeds: [DOI: 10.1137/S1064827593260140](https://doi.org/10.1137/S1064827593260140)
"""
function flux_hllc(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations1D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  e_ll  = rho_e_ll / rho_ll
  p_ll = (equations.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * v1_ll^2)
  c_ll = sqrt(equations.gamma*p_ll/rho_ll)

  v1_rr = rho_v1_rr / rho_rr
  e_rr  = rho_e_rr / rho_rr
  p_rr = (equations.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * v1_rr^2 )
  c_rr = sqrt(equations.gamma*p_rr/rho_rr)

  # Obtain left and right fluxes
  f_ll = flux(u_ll, orientation, equations)
  f_rr = flux(u_rr, orientation, equations)

  # Compute Roe averages
  sqrt_rho_ll = sqrt(rho_ll)
  sqrt_rho_rr = sqrt(rho_rr)
  sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr
  vel_L = v1_ll
  vel_R = v1_rr
  vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho
  ekin_roe = 0.5 * vel_roe^2
  H_ll = (rho_e_ll + p_ll) / rho_ll
  H_rr = (rho_e_rr + p_rr) / rho_rr
  H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
  c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe))

  Ssl = min(vel_L - c_ll, vel_roe - c_roe)
  Ssr = max(vel_R + c_rr, vel_roe + c_roe)
  sMu_L = Ssl - vel_L
  sMu_R = Ssr - vel_R
  if Ssl >= 0.0
    f1 = f_ll[1]
    f2 = f_ll[2]
    f3 = f_ll[3]
  elseif Ssr <= 0.0
    f1 = f_rr[1]
    f2 = f_rr[2]
    f3 = f_rr[3]
  else
    SStar = (p_rr - p_ll + rho_ll*vel_L*sMu_L - rho_rr*vel_R*sMu_R) / (rho_ll*sMu_L - rho_rr*sMu_R)
    if Ssl <= 0.0 <= SStar
      densStar = rho_ll*sMu_L / (Ssl-SStar)
      enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L))
      UStar1 = densStar
      UStar2 = densStar*SStar
      UStar3 = densStar*enerStar

      f1 = f_ll[1]+Ssl*(UStar1 - rho_ll)
      f2 = f_ll[2]+Ssl*(UStar2 - rho_v1_ll)
      f3 = f_ll[3]+Ssl*(UStar3 - rho_e_ll)
    else
      densStar = rho_rr*sMu_R / (Ssr-SStar)
      enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R))
      UStar1 = densStar
      UStar2 = densStar*SStar
      UStar3 = densStar*enerStar

      #end
      f1 = f_rr[1]+Ssr*(UStar1 - rho_rr)
      f2 = f_rr[2]+Ssr*(UStar2 - rho_v1_rr)
      f3 = f_rr[3]+Ssr*(UStar3 - rho_e_rr)
    end
  end
  return SVector(f1, f2, f3)
end

@inline function max_abs_speeds(u, equations::CompressibleEulerEquations1D)
  rho, rho_v1, rho_e = u
  v1 = rho_v1 / rho
  p = (equations.gamma - 1) * (rho_e - 1/2 * rho * v1^2)
  c = sqrt(equations.gamma * p / rho)

  return (abs(v1) + c,)
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::CompressibleEulerEquations1D)
  rho, rho_v1, rho_e = u

  v1 = rho_v1 / rho
  p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)

  return SVector(rho, v1, p)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equations::CompressibleEulerEquations1D)
  rho, rho_v1, rho_e = u

  v1 = rho_v1 / rho
  v_square = v1^2
  p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square)
  s = log(p) - equations.gamma*log(rho)
  rho_p = rho / p

  w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square
  w2 = rho_p * v1
  w3 = -rho_p

  return SVector(w1, w2, w3)
end

@inline function entropy2cons(w, equations::CompressibleEulerEquations1D)
  # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD
  # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1)
  @unpack gamma = equations

  # convert to entropy `-rho * s` used by Hughes, France, Mallet (1986)
  # instead of `-rho * s / (gamma - 1)`
  V1, V2, V5 = w .* (gamma - 1)

  # specific entropy, eq. (53)
  s = gamma - V1 + 0.5 * (V2^2) / V5

  # eq. (52)
  energy_internal = ((gamma - 1) / (-V5)^gamma)^(equations.inv_gamma_minus_one) * exp(-s * equations.inv_gamma_minus_one)

  # eq. (51)
  rho    = -V5 * energy_internal
  rho_v1 = V2 * energy_internal
  rho_e  = (1 - 0.5 * (V2^2) / V5) * energy_internal
  return SVector(rho, rho_v1, rho_e)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::CompressibleEulerEquations1D)
  rho, v1, p = prim
  rho_v1 = rho * v1
  rho_e  = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1)
  return SVector(rho, rho_v1, rho_e)
end

@inline function density(u, equations::CompressibleEulerEquations1D)
 rho = u[1]
 return rho
end

@inline function pressure(u, equations::CompressibleEulerEquations1D)
 rho, rho_v1, rho_e = u
 p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2) / rho)
 return p
end

@inline function density_pressure(u, equations::CompressibleEulerEquations1D)
 rho, rho_v1, rho_e = u
 rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2))
 return rho_times_p
end

# Calculate thermodynamic entropy for a conservative state `cons`
@inline function entropy_thermodynamic(cons, equations::CompressibleEulerEquations1D)
  # Pressure
  p = (equations.gamma - 1) * (cons[3] - 1/2 * (cons[2]^2) / cons[1])

  # Thermodynamic entropy
  s = log(p) - equations.gamma*log(cons[1])

  return s
end

# Calculate mathematical entropy for a conservative state `cons`
@inline function entropy_math(cons, equations::CompressibleEulerEquations1D)
  # Mathematical entropy
  S = -entropy_thermodynamic(cons, equations) * cons[1] * equations.inv_gamma_minus_one

  return S
end

# Default entropy is the mathematical entropy
@inline entropy(cons, equations::CompressibleEulerEquations1D) = entropy_math(cons, equations)

# Calculate total energy for a conservative state `cons`
@inline energy_total(cons, ::CompressibleEulerEquations1D) = cons[3]

# Calculate kinetic energy for a conservative state `cons`
@inline function energy_kinetic(cons, equations::CompressibleEulerEquations1D)
  return 0.5 * (cons[2]^2)/cons[1]
end

# Calculate internal energy for a conservative state `cons`
@inline function energy_internal(cons, equations::CompressibleEulerEquations1D)
  return energy_total(cons, equations) - energy_kinetic(cons, equations)
end

end # @muladd


text:
module Jupyter2Pluto
using JSON
using UUIDs

const _notebook_header = "### A Pluto.jl notebook ###"
const _cell_id_delimiter = "# ╔═╡ "
const _order_delimiter = "# ╠═"
const _order_delimiter_folded = "# ╟─"
const _cell_suffix = "\n\n"

export pluto2jupyter, jupyter2pluto

abstract type JupyterCell end

struct JupyterCodeCell <: JupyterCell
    execution_count::Union{Int, Nothing}
    content::Vector{String}
end
struct JupyterMarkdownCell <: JupyterCell
    content::String
end

abstract type PlutoCell end

struct PlutoCodeCell <: PlutoCell
    cell_id::UUID
    content::String
end

struct PlutoBlockCodeCell <: PlutoCell
    cell_id::UUID
    content::String
end

struct PlutoMarkdownCell <: PlutoCell
    cell_id::UUID
    content::String
end

function Base.show(io::IO, cell::PlutoCodeCell)
    gl = """
$(_cell_id_delimiter)$(cell.cell_id)
$(cell.content)$(_cell_suffix)"""
    print(io, gl)
end

function Base.show(io::IO, cell::PlutoBlockCodeCell)
    gl = """
$(_cell_id_delimiter)$(cell.cell_id)
begin\n$(cell.content)\nend$(_cell_suffix)"""
    print(io, gl)
end

function Base.show(io::IO, cell::PlutoMarkdownCell)
    content = replace(cell.content, "\"" => "\\\"")
    r = """
$(_cell_id_delimiter)$(cell.cell_id)
md\"\"\"
$(content)\"\"\"$(_cell_suffix)"""
    print(io, r)
end

function Base.Dict(cell::JupyterMarkdownCell)
    source = String[]
    content_list = split(cell.content, "\n")
    for (i, content) in enumerate(content_list)
        i == length(cell.content) && (push!(source, content); continue)
        push!(source, content*"\n")
    end
    Dict(
        "cell_type" => "markdown",
        "metadata"=>Dict(),
        "source"=> source
    )
end

function Base.Dict(cell::JupyterCodeCell)
    source = String[]
    for (i, content) in enumerate(cell.content)
        (isempty(content) || i == length(cell.content)) && (push!(source, content); continue)
        push!(source, content*"\n")
    end
    Dict(
        "cell_type" => "code",
        "execution_count" => cell.execution_count,
        "metadata"=>Dict(),
            "outputs"=>[],
            "source"=> source
    )
end

function cell_id()
    uuid4()
end

function generate_code(pcells::Vector{PlutoCell})
    pheader = _notebook_header*"\n\n"
    codestring = "$pheader$(join(pcells))$(order_code(pcells))"
    return codestring
end

function order_code(pcells::Vector{PlutoCell})
    orderstrprefix = "$(_cell_id_delimiter)Cell order:\n"
    orderstrprefix *join(order_code.(pcells))
end

function order_code(pcell::PlutoMarkdownCell)
    return _order_delimiter_folded*""*string(pcell.cell_id)*"\n"
end

function order_code(pcell::Union{PlutoCodeCell, PlutoBlockCodeCell})
   return _order_delimiter*""*string(pcell.cell_id)*"\n"
end

function generate_plutocells(codecell::JupyterCodeCell)
    pcells = PlutoCell[]
    for content in codecell.content
        id = cell_id()
        push!(pcells, PlutoCodeCell(id, content))
    end
    return pcells
end

function PlutoCell(codecell::JupyterCodeCell)
    id = cell_id()
    if has_multiple_expressions(codecell)
        return PlutoBlockCodeCell(id, join("    ".*codecell.content))
    end
    return PlutoCodeCell(id, join(codecell.content))
end

function PlutoCell(codecell::JupyterMarkdownCell)
    id = cell_id()
    return PlutoMarkdownCell(id, codecell.content)
end

function has_multiple_expressions(codecell::JupyterCodeCell)
    expressions = Meta.parse("begin "*join(codecell.content, "\n")*" end").args
    filter!(x -> !(x isa LineNumberNode), expressions)
    length(expressions) > 1
end

function jupyter2pluto(jupyter_file)
    jupyter_cells = try
        content = JSON.parsefile(jupyter_file)
        content["cells"]
    catch ex
        error("Jupyter notebook parse error")
    end
    pluto_cells = PlutoCell[]
    for cell in jupyter_cells
        !haskey(cell, "cell_type") || !haskey(cell,"source") && continue
        if cell["cell_type"] == "markdown" && !isempty(cell["source"])
            jmark_cell = JupyterMarkdownCell(join(cell["source"]))
            pc = PlutoCell(jmark_cell)
            push!(pluto_cells, pc)
        end
        if cell["cell_type"] == "code" && !isempty(cell["source"])
            jcode_cell = JupyterCodeCell(cell["execution_count"],cell["source"])
            pc = PlutoCell(jcode_cell)
            push!(pluto_cells, pc)
        end
    end
    output = generate_code(pluto_cells)
    dest = "$(jupyter_file).jl"
    open(dest, "w") do f
        write(f, output)
    end
    dest
end

function parse_pluto_cell(rawcell::String)
    cellist = string.(split(rawcell, '\n'))
    body = join(cellist[2:end], '\n')
    multiline_mdr = r"md\"\"\"(.*)\"\"\""s
    line_mdr = r"md\"(.*)\""s
    matches = if (mat = match(multiline_mdr, body)) != nothing
        mat
    else
        match(line_mdr, body)
    end

    if matches != nothing
        PlutoMarkdownCell(UUID(cellist[1]), matches.captures[1])
    else
        PlutoCodeCell(UUID(cellist[1]), body)
    end
end

function parse_pluto_end(rawcell::String)
    main_split= string.(split(rawcell, _order_delimiter))
    orderids = String[]
    for order in main_split
        splits = string.(split(order, _order_delimiter_folded))
        for split in splits
            !isempty(split) && push!(orderids, strip(split))
        end
    end
    orderids[2:end]
end

function order(pcells::Vector{PlutoCell}, orderids::Vector{String})
    sorted_pcells = PlutoCell[]
    for orderid in orderids
        for pcell in pcells
            if occursin(string(pcell.cell_id),  orderid)
                push!(sorted_pcells, pcell)
            end
        end
    end
    return sorted_pcells
end

let
    global JupyterCell
    global parse_pluto_load
    execution_count = 1
    function parse_pluto_load(raw::AbstractString)
        jcell = JupyterCodeCell(execution_count, split(raw, "\n"))
        execution_count += 1
        jcell
    end
    function JupyterCell(pcell::PlutoCodeCell)
        jcell = JupyterCodeCell(execution_count,split(pcell.content, "\n"))
        execution_count += 1
        jcell
    end
end

function JupyterCell(pcell::PlutoMarkdownCell)
    content = pcell.content
    JupyterMarkdownCell(pcell.content)
end

function pluto2jupyter(file)
    # parser: pluto notebook has orderidlist, map(order_id => codesnippets) ::Plutocells
    plutoraw = readchomp(file)
    plutocelllist = string.(split(plutoraw, _cell_id_delimiter))
    jupyterloadcell = parse_pluto_load(plutocelllist[1])
    i=2
    pcells = PlutoCell[]
    while(i <= length(plutocelllist)-1)
        pcell = parse_pluto_cell(plutocelllist[i])
        push!(pcells, pcell)
        i+=1
    end
    plutoend = parse_pluto_end(plutocelllist[end])
    ordered_pcells = order(pcells, plutoend)
    jcells = map( pcell -> JupyterCell(pcell), ordered_pcells)
    all_main_cells = Dict.(jcells)
    pushfirst!(all_main_cells, Dict(jupyterloadcell))
    d_cells = Dict()
    d_cells["cells"] = all_main_cells
    d_cells["metadata"] = Dict(
        "kernelspec"=> Dict(
            "display_name"=> "Julia $(VERSION)",
            "language"=> "julia",
            "name"=> "julia-$(VERSION.major).$(VERSION.minor)"
        ),
        "language_info"=> Dict(
            "file_extension"=> ".jl",
            "mimetype"=> "application/julia",
            "name"=> "julia",
            "version"=> string(VERSION)
        )
    )
    d_cells["nbformat"]= 4
    d_cells["nbformat_minor"]= 2
    dest = file*".ipynb"
    open(dest, "w") do f
        JSON.print(f, d_cells , 4)
    end
    dest

end
end # module


text:
#
"""
Deform domain, compute Jacobian of transformation, and its inverse

given 2D mapping

x = x(r,s), y = y(r,s)

compute

dXdR = [xr xs], dRdX = [rx ry], J = det(dXdR), Jinv = det(dRdX)
       [yr ys]         [sx sy]

"""
struct DeformedSpace{T,D,
                     Tspace<:AbstractSpace{T,D},
                     Tgrid,Tjacmat,Tjacimat,Tjac,Tjaci} <: AbstractSpace{T,D}
    space::Tspace
    grid::Tgrid # (x1, ..., xD,)
    dXdR::Tjacmat
    dRdX::Tjacimat
    J::Tjac
    Ji::Tjaci
    # add DeformedDomain
end

function deform(space::AbstractSpace{<:Number,D},
                mapping = nothing, isseparable = false
               ) where{D}
    if mappping === nothing
        J    = IdentityOp{D}()
        Jmat = Diagonal([J for i=1:D])
        @warn "mapping === nothing"
#       return space
        return DeformedSpace(space, grid(space), Jmat, Jmat, J, J)

    elseif isseparable # x = x(r), y = y(s)
        # eliminate cross terms by making dXdR, etc diagonal
    elseif rescaling # simple rescaling
        # make jac, dXdR, etc scaling operations
    end

    R = grid(space)
    X = mapping(R...)

    gradR = gradOp(space)

    """
    dXdR = [dx1/dr1 ... dx1/drD]
           [...     ...     ...]
           [dxD/dr1 ... dxD/drD]
    """
    dXdR = begin
        dXdR = gradR.(X) |> hcat
        dXdR = DiagonalOp.(dXdR)
        dXdR = dXdR'
    end

    J = if D == 1
        dXdR[1]
    elseif D == 2
        xr = dXdR[1]
        yr = dXdR[2]
        xs = dXdR[3]
        ys = dXdR[4]

        xr * ys - xs * yr
    elseif D == 3
        xr = dXdR[1]; yr = dXdR[2]; zr = dXdR[3]
        xs = dXdR[4]; ys = dXdR[5]; zs = dXdR[6]
        xt = dXdR[7]; yt = dXdR[8]; zt = dXdR[9]

        J = xr * (ys * zt - zs * yt) -
            xs * (yr * zt - zr * yt) +
            xt * (yr * zs - zr * ys)
    else
        det(dXdR) # errors
    end

    Ji = inv(D)

    dRdX = if D == 1
        fill(Ji, 1, 1)
    elseif D == 2
        xr = dXdR[1]
        yr = dXdR[2]
        xs = dXdR[3]
        ys = dXdR[4]

        rx =  (Ji * ys)
        ry = -(Ji * xs)
        sx = -(Ji * yr)
        sy =  (Ji * xr)

        dRdX = [rx ry
                sx sy]
    elseif D == 3 # cramer's rule
        inv(dXdR)
    else
        inv(dXdR) # need even more robust AbstractOp interface
    end

    DeformedSpace(space, X, dXdR, dRdX, J, Ji)
end

get_grid(space::DeformedSpace) = space.grid

"""
[Dx] * u = [rx sx] * [Dr] * u
[Dy]     = [ry sy]   [Ds]
"""
function gradOp(space::DeformedSpace)
    gradR = gradOp(space.space)
    dRdX  = space.dRdX
    gradX = dRdX * gradR

    gradX
end

function massOp(space::DeformedSpace)
    M = massOp(space.space)
    J = space.J

    M * J
end

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

implemented as

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

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

= [Dr]' * [G11 G12]' * [Dr]
  [Ds]    [G12 G22]    [Ds]
"""
function laplaceOp(space::DeformedSpace{<:Number, D}) where{D}

    Dr   = gradOp(space.space)
    M    = massOp(space)
    dRdX = space.dRdX

    MM = Diagonal([M for i=1:D])
    GG = if dRdX isa Diagonal
        dRdX' * MM * dRdX
    else
        dRdX' * MM * dRdX |> Symmetric # TODO avoid bottom 1/2 computation
    end

    laplOp = Dr' * GG * Dr

    return first(laplOp)
end

###
# Dealiased operators
###

function laplaceOp(space1::DeformedSpace{<:Number,D},
                   space2::DeformedSpace{<:Number,D},
                   J = nothing
                  ) where{D}
    J12 = J !== nothing ? J : interpOp(space1, space2)

    Dr1   = gradOp(space1.space)
    M2    = MassOp(space2)
    dRdX2 = space2.dRdX

    JD = J12 * Dr1

    MM2 = Diagonal([M2 for i=1:D])
    GG2 = dRdX' * MM2 * dRdX |> Symmetric

    laplOp = JD' * GG2 * JD

    return first(laplOp)
end
#


text:
export HTTPError

struct HTTPError <: Exception
    code::UInt16
    message::AbstractString
    data::Dict{AbstractString, Any}
    safe::Bool
    shouldclose::Bool
end

function HTTPError(code::Integer;
    message::AbstractString="",
    data::Dict{AbstractString, Any}=Dict{AbstractString, Any}(),
    safe::Bool=code<500,
    shouldclose::Bool=false)

    HTTPError(code, message, data, safe, shouldclose)
end


text:
using Base.Test, Compat

addprocs(1)

import NBInclude; @everywhere using NBInclude

@test include(joinpath("includes", "test1.jl")) == 314159
@test f(5) == 6
@test normpath(myfile) == abspath("test.ipynb") * ":In[6]"
@test normpath(myfile2) == abspath("test.ipynb") * ":In[+7]"

@test remotecall_fetch(nbinclude, 2, "test.ipynb") == 314159


text:
export Vandermonde

struct Vandermonde{T} <: AbstractArray{T, 2}
    c :: Vector{T}
end

getindex(V::Vandermonde, i::Int, j::Int) = V.c[i]^(j-1)
isassigned(V::Vandermonde, i::Int, j::Int) = isassigned(V.c, i)

size(V::Vandermonde, r::Int) = (r==1 || r==2) ? length(V.c) :
    throw(ArgumentError("Invalid dimension $r"))
size(V::Vandermonde) = length(V.c), length(V.c)

function Matrix(V::Vandermonde{T}) where T
    n=size(V, 1)
    M=Array{T}(undef, n, n)
    M[:,1] .= 1
    for j=2:n
        for i=1:n
	    M[i,j] = M[i,j-1]*V.c[i]
        end
    end
    M
end

function Matrix(V::Adjoint{T,Vandermonde{T}}) where T
    n=size(V, 1)
    M=Array{T}(undef, n, n)
    M[1,:] .= 1
    for j=1:n
        for i=2:n
	    M[i,j] = M[i-1,j]*adjoint(V.parent.c[j])
        end
    end
    M
end

function Matrix(V::Transpose{T,Vandermonde{T}}) where T
    n=size(V, 1)
    M=Array{T}(undef, n, n)
    M[1,:] .= 1
    for j=1:n
        for i=2:n
	    M[i,j] = M[i-1,j]*V.parent.c[j]
        end
    end
    M
end

function \(V::Adjoint{T1,Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where T1 where T2
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    pvand!(adjoint(V.parent.c), x)
    return x
end

function \(V::Transpose{T1,Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where T1 where T2
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    pvand!(V.parent.c, x)
    return x
end

function \(V::Vandermonde{T1}, y::AbstractVecOrMat{T2}) where T1 where T2
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    dvand!(V.c, x)
    return x
end

function vandtype(T1::Type, T2::Type)
    # Figure out the return type of Vandermonde{T1} \ Vector{T2}
    T = promote_type(T1, T2)
    S = typeof(oneunit(T)/oneunit(T1))
    return S
end

"""
    pvand!(a, b) -> b

Solves system ``A^T*x = b`` in-place.

``A^T`` is transpose of Vandermonde matrix ``A_{ij} = a_i^{j-1}``.

Algorithm by Bjorck & Pereyra,
Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903,
https://doi.org/10.2307/2004623
"""
function pvand!(alpha, B)
    n = length(alpha);
    if n != size(B,1)
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    end
    nrhs = size(B,2)
    @inbounds begin
        for j=1:nrhs
            for k=1:n-1
                for i=n:-1:k+1
                    B[i,j] = B[i,j]-alpha[k]*B[i-1,j]
                end
            end
            for k=n-1:-1:1
                for i=k+1:n
                    B[i,j] = B[i,j]/(alpha[i]-alpha[i-k])
                end
                for i=k:n-1
                    B[i,j] = B[i,j]-B[i+1,j]
                end
            end
        end
    end
end
"""
    dvand!(a, b) -> b

Solves system ``A*x = b`` in-place.

``A`` is Vandermonde matrix ``A_{ij} = a_i^{j-1}``.

Algorithm by Bjorck & Pereyra,
Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903,
https://doi.org/10.2307/2004623
"""
function dvand!(alpha, B)
    n = length(alpha)
    if n != size(B,1)
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    end
    nrhs = size(B,2)
    @inbounds begin
        for j=1:nrhs
            for k=1:n-1
                for i=n:-1:k+1
                    B[i,j] = (B[i,j]-B[i-1,j])/(alpha[i]-alpha[i-k])
                end
            end
            for k=n-1:-1:1
                for i=k:n-1
                    B[i,j] = B[i,j]-alpha[k]*B[i+1,j]
                end
            end
        end
    end
end


text:
# ---
# title: 916. Word Subsets
# id: problem916
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/word-subsets/description/>
# hidden: true
# ---
#
# We are given two arrays `A` and `B` of words.  Each word is a string of
# lowercase letters.
#
# Now, say that word `b` is a subset of word `a` ** ** if every letter in `b`
# occurs in `a`, **including multiplicity**.  For example, `"wrr"` is a subset
# of `"warrior"`, but is not a subset of `"world"`.
#
# Now say a word `a` from `A` is _universal_ if for every `b` in `B`, `b` is a
# subset of `a`.
#
# Return a list of all universal words in `A`.  You can return the words in any
# order.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = ["amazon","apple","facebook","google","leetcode"], B = ["e","o"]
#     Output: ["facebook","google","leetcode"]
#
#
# **Example 2:**
#
#
#
#     Input: A = ["amazon","apple","facebook","google","leetcode"], B = ["l","e"]
#     Output: ["apple","google","leetcode"]
#
#
# **Example 3:**
#
#
#
#     Input: A = ["amazon","apple","facebook","google","leetcode"], B = ["e","oo"]
#     Output: ["facebook","google"]
#
#
# **Example 4:**
#
#
#
#     Input: A = ["amazon","apple","facebook","google","leetcode"], B = ["lo","eo"]
#     Output: ["google","leetcode"]
#
#
# **Example 5:**
#
#
#
#     Input: A = ["amazon","apple","facebook","google","leetcode"], B = ["ec","oc","ceo"]
#     Output: ["facebook","leetcode"]
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length, B.length <= 10000`
#   2. `1 <= A[i].length, B[i].length <= 10`
#   3. `A[i]` and `B[i]` consist only of lowercase letters.
#   4. All words in `A[i]` are unique: there isn't `i != j` with `A[i] == A[j]`.
#
#
## @lc code=start
using LeetCode

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


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetVMProfile <: SwaggerModel
    osProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetOSProfile } # spec name: osProfile
    storageProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetStorageProfile } # spec name: storageProfile
    networkProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetNetworkProfile } # spec name: networkProfile
    securityProfile::Any # spec type: Union{ Nothing, SecurityProfile } # spec name: securityProfile
    diagnosticsProfile::Any # spec type: Union{ Nothing, DiagnosticsProfile } # spec name: diagnosticsProfile
    extensionProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetExtensionProfile } # spec name: extensionProfile
    licenseType::Any # spec type: Union{ Nothing, String } # spec name: licenseType
    priority::Any # spec type: Union{ Nothing, Priority } # spec name: priority
    evictionPolicy::Any # spec type: Union{ Nothing, EvictionPolicy } # spec name: evictionPolicy
    billingProfile::Any # spec type: Union{ Nothing, BillingProfile } # spec name: billingProfile
    scheduledEventsProfile::Any # spec type: Union{ Nothing, ScheduledEventsProfile } # spec name: scheduledEventsProfile

    function VirtualMachineScaleSetVMProfile(;osProfile=nothing, storageProfile=nothing, networkProfile=nothing, securityProfile=nothing, diagnosticsProfile=nothing, extensionProfile=nothing, licenseType=nothing, priority=nothing, evictionPolicy=nothing, billingProfile=nothing, scheduledEventsProfile=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("osProfile"), osProfile)
        setfield!(o, Symbol("osProfile"), osProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("storageProfile"), storageProfile)
        setfield!(o, Symbol("storageProfile"), storageProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("networkProfile"), networkProfile)
        setfield!(o, Symbol("networkProfile"), networkProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("securityProfile"), securityProfile)
        setfield!(o, Symbol("securityProfile"), securityProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("diagnosticsProfile"), diagnosticsProfile)
        setfield!(o, Symbol("diagnosticsProfile"), diagnosticsProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("extensionProfile"), extensionProfile)
        setfield!(o, Symbol("extensionProfile"), extensionProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("licenseType"), licenseType)
        setfield!(o, Symbol("licenseType"), licenseType)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("priority"), priority)
        setfield!(o, Symbol("priority"), priority)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("evictionPolicy"), evictionPolicy)
        setfield!(o, Symbol("evictionPolicy"), evictionPolicy)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("billingProfile"), billingProfile)
        setfield!(o, Symbol("billingProfile"), billingProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("scheduledEventsProfile"), scheduledEventsProfile)
        setfield!(o, Symbol("scheduledEventsProfile"), scheduledEventsProfile)
        o
    end
end # type VirtualMachineScaleSetVMProfile

const _property_map_VirtualMachineScaleSetVMProfile = Dict{Symbol,Symbol}(Symbol("osProfile")=>Symbol("osProfile"), Symbol("storageProfile")=>Symbol("storageProfile"), Symbol("networkProfile")=>Symbol("networkProfile"), Symbol("securityProfile")=>Symbol("securityProfile"), Symbol("diagnosticsProfile")=>Symbol("diagnosticsProfile"), Symbol("extensionProfile")=>Symbol("extensionProfile"), Symbol("licenseType")=>Symbol("licenseType"), Symbol("priority")=>Symbol("priority"), Symbol("evictionPolicy")=>Symbol("evictionPolicy"), Symbol("billingProfile")=>Symbol("billingProfile"), Symbol("scheduledEventsProfile")=>Symbol("scheduledEventsProfile"))
const _property_types_VirtualMachineScaleSetVMProfile = Dict{Symbol,String}(Symbol("osProfile")=>"VirtualMachineScaleSetOSProfile", Symbol("storageProfile")=>"VirtualMachineScaleSetStorageProfile", Symbol("networkProfile")=>"VirtualMachineScaleSetNetworkProfile", Symbol("securityProfile")=>"SecurityProfile", Symbol("diagnosticsProfile")=>"DiagnosticsProfile", Symbol("extensionProfile")=>"VirtualMachineScaleSetExtensionProfile", Symbol("licenseType")=>"String", Symbol("priority")=>"Priority", Symbol("evictionPolicy")=>"EvictionPolicy", Symbol("billingProfile")=>"BillingProfile", Symbol("scheduledEventsProfile")=>"ScheduledEventsProfile")
Base.propertynames(::Type{ VirtualMachineScaleSetVMProfile }) = collect(keys(_property_map_VirtualMachineScaleSetVMProfile))
Swagger.property_type(::Type{ VirtualMachineScaleSetVMProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetVMProfile[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetVMProfile }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetVMProfile[property_name]

function check_required(o::VirtualMachineScaleSetVMProfile)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetVMProfile }, name::Symbol, val)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.8.1
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# +
using GLM
using CSV
using Random
using StatsBase
using DataFrames
using Dagitty
using Turing
using StatsPlots
using StatisticalRethinking, StructuralCausalModels
using StatisticalRethinkingPlots
using Cairo, Fontconfig
using Logging

default(labels = false)
Logging.disable_logging(Logging.Warn);
# -

# Code 6.1

# +
Random.seed!(1917)
N = 200   # grant proposals
p = 0.1   # proportion to select

# uncorrelated newsworthiness and trustworthiness
nw = rand(Normal(), N)
tw = rand(Normal(), N)

# select top 10% of combined score
s = nw .+ tw
q = quantile(s, 1 - p)
selected = s .>= q
cor(tw[selected], nw[selected])
# -

scatter(nw[.!selected], tw[.!selected]; xlab = "newsworthiness", ylab = "trustworthiness", label = "not selected")
scatter!(nw[selected], tw[selected]; label = "selected")

# # 6.1 Multicollinearity

# Code 6.2

Random.seed!(100)
N = 100
height = rand(Normal(10, 2), N)
leg_prop = rand(Uniform(0.4, 0.5), N)
leg_left = leg_prop .* height .+ rand(Normal(0, 0.02), N)
leg_right = leg_prop .* height .+ rand(Normal(0, 0.02), N)
d = DataFrame(:height => height, :leg_left => leg_left, :leg_right => leg_right);
first(d, 6)

# Code 6.3

# +
@model function model_m6_1(leg_left, leg_right, height)
    a ~ Normal(10, 100)
    bl ~ Normal(2, 10)
    br ~ Normal(2, 10)
    μ = @. a + bl * leg_left + br * leg_right
    σ ~ Exponential(1)
    height ~ MvNormal(μ, σ)
end

m6_1 = sample(model_m6_1(d.leg_left, d.leg_right, d.height), NUTS(), 1000)
m6_1_df = DataFrame(m6_1);
precis(m6_1_df)
# -

# Code 6.4

coeftab_plot(m6_1_df)

# Code 6.5

scatter(m6_1_df.br, m6_1_df.bl; alpha = 0.1)

# Code 6.6

@df m6_1_df density(:br + :bl; lw = 2, xlab = "sum of bl and br")

# Code 6.7

# +
@model function model_m6_2(leg_left, height)
    a ~ Normal(10, 100)
    bl ~ Normal(2, 10)
    μ = @. a + bl * leg_left
    σ ~ Exponential(1)
    height ~ MvNormal(μ, σ)
end

m6_2 = sample(model_m6_2(d.leg_left, d.height), NUTS(), 1000)
m6_2_df = DataFrame(m6_2);
precis(m6_2_df)
# -

std(m6_1_df.bl), std(m6_1_df.br), std(m6_1_df.bl + m6_1_df.br)

# Code 6.8

# +
d = DataFrame(CSV.File("data/milk.csv", missingstring = "NA"))

# get rid of dots in column names
rename!(n -> replace(n, "." => "_"), d)

d[!, :K] = standardize(ZScoreTransform, d.kcal_per_g)
d[!, :F] = standardize(ZScoreTransform, d.perc_fat)
d[!, :L] = standardize(ZScoreTransform, d.perc_lactose);
# -

# Code 6.9

# +
@model function model_m6_3(F, K)
    a ~ Normal(0, 0.2)
    bF ~ Normal(0, 0.5)
    μ = @. a + F * bF
    σ ~ Exponential(1)
    K ~ MvNormal(μ, σ)
end

m6_3 = sample(model_m6_3(d.F, d.K), NUTS(500, 0.65), 1000)
m6_3_df = DataFrame(m6_3);

@model function model_m6_4(L, K)
    a ~ Normal(0, 0.2)
    bL ~ Normal(0, 0.5)
    μ = @. a + L * bL
    σ ~ Exponential(1)
    K ~ MvNormal(μ, σ)
end

m6_4 = sample(model_m6_4(d.L, d.K), NUTS(), 1000)
m6_4_df = DataFrame(m6_4);

precis(m6_3_df)
precis(m6_4_df)

coeftab_plot(m6_3_df)
coeftab_plot(m6_4_df)
# -

# Code 6.10

# +
@model function model_m6_5(F, L, K)
    a ~ Normal(0, 0.2)
    bF ~ Normal(0, 0.5)
    bL ~ Normal(0, 0.5)
    μ = @. a + F * bF + L * bL
    σ ~ Exponential(1)
    K ~ MvNormal(μ, σ)
end

m6_5 = sample(model_m6_5(d.F, d.L, d.K), NUTS(), 1000)
m6_5_df = DataFrame(m6_5);
precis(m6_5_df)
# -

# Code 6.11

@df d corrplot([:kcal_per_g :perc_fat :perc_lactose]; seriestype = :scatter, bins = 10, grid = false)

# Code 6.12

# +
# get mean stderr for linear model's scale
function stderr_for_r(r)
    σ = sqrt(1 - r^2) * var(d.perc_fat)
    fat_scaled = r .* d.perc_fat
    stderr_x = [
        begin
            x = d.perc_fat .+ rand(MvNormal(fat_scaled, σ))
            # add the intercept to the model
            X = hcat(ones(length(x)), x)
            m = lm(X, d.kcal_per_g)
            stderror(m)[2]
        end
        for _ = 1:100
    ]
    s = mean(stderr_x)
end

r_seq = range(0, 0.99; step = 0.01)
s = stderr_for_r.(r_seq)
plot(r_seq, s; lw = 2, xlab = "correlation")
# -

# # 6.2 Post-treatment bias

# Code 6.13

# +
Random.seed!(70)
# number of plants
N = 100
h0 = rand(Normal(10, 2), N)
treatment = repeat(0:1, inner = div(N, 2))
fungus = [rand(Binomial(1, 0.5 - treat * 0.4)) for treat in treatment]
h1 = h0 .+ rand(MvNormal(5 .- 3 .* fungus, 1))

d = DataFrame(:h0 => h0, :h1 => h1, :treatment => treatment, :fungus => fungus)
precis(d)
# -

# Code 6.14

sim_p = rand(LogNormal(0, 0.25), 10_000)
precis(DataFrame(:sim_p => sim_p))

# Code 6.15

# +
@model function model_m6_6(h0, h1)
    p ~ LogNormal(0, 0.25)
    σ ~ Exponential(1)
    μ = h0 .* p
    h1 ~ MvNormal(μ, σ)
end

m6_6 = sample(model_m6_6(d.h0, d.h1), NUTS(), 1000)
m6_6_df = DataFrame(m6_6);
precis(m6_6_df)
# -

# Code 6.16

# +
@model function model_m6_7(h0, treatment, fungus, h1)
    a ~ LogNormal(0, 0.2)
    bt ~ Normal(0, 0.5)
    bf ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    p = @. a + bt * treatment + bf * fungus
    μ = h0 .* p
    h1 ~ MvNormal(μ, σ)
end

m6_7 = sample(model_m6_7(d.h0, d.treatment, d.fungus, d.h1), NUTS(), 1000)
m6_7_df = DataFrame(m6_7);
precis(m6_7_df)
coeftab_plot(m6_7_df)
# -

# Code 6.17

# +
@model function model_m6_8(h0, treatment, h1)
    a ~ LogNormal(0, 0.2)
    bt ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    p = @. a + bt * treatment
    μ = h0 .* p
    h1 ~ MvNormal(μ, σ)
end

m6_8 = sample(model_m6_8(d.h0, d.treatment, d.h1), NUTS(), 1000)
m6_8_df = DataFrame(m6_8);
precis(m6_8_df)
coeftab_plot(m6_8_df)
# -

# Code 6.18

plant_dag = Dagitty.DAG(:H₀ => :H₁, :F => :H₁, :T => :F)
drawdag(plant_dag, [2, 0, 1, 3], [0, 0, 0, 0])

# Code 6.19

implied_conditional_independencies_min(plant_dag)

# Code 6.20

# +
Random.seed!(70)
# number of plants
N = 1000
h0 = rand(Normal(10, 2), N)
treatment = repeat(0:1, inner = div(N, 2))
M = rand(Bernoulli(), N)
fungus = [
    rand(Binomial(1, 0.5 - treat * 0.4 + 0.4 * m))
    for (treat, m) ∈ zip(treatment, M)
]
h1 = h0 .+ rand(MvNormal(5 .+ 3 .* M, 1))

d2 = DataFrame(:h0 => h0, :h1 => h1, :treatment => treatment, :fungus => fungus);
precis(d2)
# -

m6_7 = sample(model_m6_7(d2.h0, d2.treatment, d2.fungus, d2.h1), NUTS(), 1000)
precis(DataFrame(m6_7))

m6_8 = sample(model_m6_8(d2.h0, d2.treatment, d2.h1), NUTS(), 1000)
precis(DataFrame(m6_8))

# # 6.3 Collider bias

# Code 6.21

d = sim_happiness(seed = 1977, n_years = 1000)
precis(d)

# +
d_m = d[d.married.==1, [:age, :happiness]]
d_u = d[d.married.==0, [:age, :happiness]]

scatter(d_m.age, d_m.happiness; label = "married", xlab = "age", ylab = "happiness")
scatter!(d_u.age, d_u.happiness; c = :white)
# -

# Code 6.22

d2 = d[d.age.>17, :]
d2[!, :A] = @. (d2.age - 18) / (65 - 18);

# Code 6.23

d2[!, :mid] = d2.married .+ 1;

# +
@model function model_m6_9(mid, A, happiness)
    a ~ MvNormal([0, 0], 1)
    bA ~ Normal(0, 2)
    μ = a[mid] .+ bA .* A
    σ ~ Exponential(1)
    happiness ~ MvNormal(μ, σ)
end

m6_9 = sample(model_m6_9(d2.mid, d2.A, d2.happiness), NUTS(), 1000)
m6_9_df = DataFrame(m6_9);
precis(m6_9_df)
# -

# Code 6.24

# +
@model function model_m6_10(A, happiness)
    a ~ Normal()
    bA ~ Normal(0, 2)
    μ = a .+ bA .* A
    σ ~ Exponential(1)
    happiness ~ MvNormal(μ, σ)
end

m6_10 = sample(model_m6_10(d2.A, d2.happiness), NUTS(), 1000)
m6_10_df = DataFrame(m6_10)
precis(m6_10_df)
# -

# Code 6.25

N = 200
b_GP = 1
b_GC = 0
b_PC = 1
b_U = 2

# Code 6.26

Random.seed!(42)
U = 2 .* rand(Bernoulli(), N) .- 1
G = rand(Normal(), N)
P = rand(MvNormal(@. b_GP * G + b_U * U))
C = rand(MvNormal(@. b_PC * P + b_GC * G + b_U * U))
d = DataFrame(:C => C, :P => P, :G => G, :U => U);

# Code 6.27

# +
@model function model_m6_11(P, G, C)
    a ~ Normal()
    b_PC ~ Normal()
    b_GC ~ Normal()
    μ = @. a + b_PC * P + b_GC * G
    σ ~ Exponential(1)
    C ~ MvNormal(μ, σ)
end

m6_11 = sample(model_m6_11(d.P, d.G, d.C), NUTS(), 1000)
m6_11_df = DataFrame(m6_11);
precis(m6_11_df)
# -

# Code 6.28

# +
@model function model_m6_12(P, G, U, C)
    a ~ Normal(0, 1)
    b_PC ~ Normal()
    b_GC ~ Normal()
    b_U ~ Normal()
    μ = @. a + b_PC * P + b_GC * G + b_U * U
    σ ~ Exponential(1)
    C ~ MvNormal(μ, σ)
end

m6_12 = sample(model_m6_12(d.P, d.G, d.U, d.C), NUTS(), 1000)
m6_12_df = DataFrame(m6_12);
precis(m6_12_df)

# NOTE: something not quite right here, as we should obtain the same values as the parameters used for generating the data. Does Turing scale by default?
# -

# # 6.4 Confronting confounding

# Code 6.29

# +
dag_61 = Dagitty.DAG(
    :X => :Y,
    :U => :X, :A => :U,
    :A => :C, :C => :Y,
    :U => :B, :C => :B,
)

all_backdoor_adjustment_sets(dag_61, :X, :Y)

dag_61 = StructuralCausalModels.DAG("dag_61",
    "DAG(
    U ~ A,
    C ~ A,
    B ~ U + C,
    X ~ U,
    Y ~ X + C
    )"
)
StructuralCausalModels.adjustment_sets(dag_61, :X, :Y)
# -

# Code 6.30

dag_62 = Dagitty.DAG(
    :A => :D,
    :A => :M, :M => :D,
    :S => :A, :S => :M,
    :S => :W, :W => :D,
)
all_backdoor_adjustment_sets(dag_62, :W, :D)

implied_conditional_independencies_min(dag_62)

#  StructuralCausalModels version:
dag_62b = StructuralCausalModels.DAG("dag_62b",
    "DAG(
        A ~ S,
        M ~ A+S,
        W ~ S,
        D ~ A + M + W
    )"
)
basis_set(dag_62b)
adjustment_sets(dag_62b, :W, :D)


text:
function alt_constraint(mdl::Model)

    R, Rᴸ = mdl.Prices[1:2]
    ω, τ = mdl.AddPrices

    #initialize VFalter if not done already
    size(mdl.VFalter) != size(mdl.VF) && (mdl.VFalter = zeros(Float64, size(mdl.VF)))

    #Indices for get_vars
    ndim = length.( mdl.GridDef.Grids[1:(end-1)] )
    cins = CartesianIndices( tuple( ndim[1:(length(ndim)÷2)]... ) )
    # cins = CartesianIndex.(cins,cins)

    #utility today if defaulted
    #TODO: This does not have to be calculated here, each call. It can be precalculated once after prices change.Find solution and fix it
    u_alter = zeros(Float64, size(mdl.VF) )
    vfalter = view(mdl.VFalter,:,:)
    for status = mdl.GridDef.Grids[end]

        @Threads.threads for idx = cins

            v = get_vars( (Tuple(CartesianIndex(idx,idx))... ,status), mdl_init.GridDef)
            dᴿ =  v.l₁ * δᶠ * δˡ;

            #Consumption today if defaulted
            c = max( (( v.d₁ - dᴿ) * ( ρ * R / 𝜋) + ( 1 - τ ) * ω * (v.λᵢ==1) + ϐ * (v.λᵢ==2)), ϐ)
            u_alter[LinearIndices(cins)[idx],status] = mdl.HAutility(c, true,true,false)

        end

    end

    c_if_defaulted = Χ * [(( 1 - τ ) * ω), ϐ]
    vfalter[:,:] = vfalter[:,:] - u_alter # TODO: se todo above

    # vfalter = mdl.HAutility.( c_if_defaulted, 0.0 ) +
    vfalter[:,:] = u_alter .+
        mdl.β * (
        (1 - ζ) * ( ( mdl.λ * mdl.HAutility.(c_if_defaulted, true, true, false))' .+
        mdl.β * ( vfalter[:,:] * mdl.λ' ) ) .+ ζ * ( mdl.λ * mdl.VF[1,:] )'
        )

    return(mdl.VFalter)

end


text:
using ReachabilityAnalysis, SparseArrays
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=5e-3)),
    "X0"   => Hyperrectangle(spzeros(28), sparsevec(1:8, fill(0.1, 8), 28)),
    "info" => Dict(
        "name" => "helicopter",
        "dim" => 28,
        "linear" => true,
        "continuous" => true,
        "domain" => "Automatic Control",
        "hybrid" => false,
        "description" => "An 8-dimensional helicopter model with a 20-dimensional controller. The controlled plant is a 28-dimensional continuous linear time-invariant (LTI) system. The plant is a small disturbance model of a helicopter, given as an 8-dimensional LTI system. The controller we examine is an H∞ mixed-sensitivity design for rejecting atmospheric turbulence, given as a 20-dimensional LTI system. The helicopter is representative of the Westland Lynx, a two-engined multi-purpose military helicopter of approximately 9000 lbs gross weight, with a four-blade semi-rigid main rotor.",
        "source" => "S. Skogestad and I. Postlethwaite. Multivariable Feedback Control: Analysis and Design. John Wiley & Sons, 2005. The model is taken from the Example Mdels from SpaceEx state space explorer (see http://spaceex.imag.fr/download-6).",
        "varnames" => Dict("x₁" => "Pitch Attitude (Θ)", "x₂" => "Roll Attitude (Φ)", "x₃" => "Roll Rate (p)", "x₄" => "Pitch Rate (q)",
                           "x₅" => "Yaw Rate (r)", "x₆" => "Forward Velocity (vx)", "x₇" => "Lateral Velocity (vy)", "x₈" => "Vertical Velocity (vz)"),
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
#=
DA_reflslab is the first iteration to compute the TD reflectance from a slab.
This code is very 'MATLABY'... as in you generally want to avoid loops in MATLAB and vectorize everything.
Though it is slow, I find this very composable and readable (I feel there should always be a trade off)
With that said the more optimized version is ~2.5x faster with the ability to scale much better.
The main performance increase came from completely devectorizing the code and iteratively looping through
the summation and time vector.
Some improvement also came by imposing the @inbounds macro which skips the checking bounds.
The primary advantage of creating loops and devectorizing is it allows multi-threading.
The Threads.@threads runs the for loops on multiple threads on your cpu. This ran on 4 threads but as mentioned
	this operation is easily scalable with more threads in CPU.
*** Julia's default is 1 thread where MATLAB I believe defaults to largest available. You have to initiate Julia
running multiple threads....
=#

function DA_reflslab(t, β::Array{Float64,1}, ρ::Float64,ndet::Float64, nmed::Float64, s::Float64)
	#s is slab thickness

	n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	z3m = Array{Float64}(undef, length(xs))
	z4m = Array{Float64}(undef, length(xs))

    Rt1 = Array{Float64}(undef, length(t))
    Rt2 = Array{Float64}(undef, length(t), length(z3m))

	if n == 1.0
		A= 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

    z3m = Float64[-2m.*s .- 4m.*ze .- zs for m in xs]
	z4m = Float64[-2m.*s .- (4m .- 2).*ze .+ zs for m in xs]

    Rt1 = @. -exp(-(ρ^2/(4D*ν*t)) - μa*ν*t)
   	Rt1 = @. Rt1/(2*(4π*D*ν)^(3/2)*t^(5/2))

	Rt2 = @. z3m'*exp(-(z3m'^2 / (4D*ν*t))) - z4m'*exp(-(z4m'^2 / (4D*ν*t)))

    Rt = Rt1.*sum(Rt2, dims=2)
	replace!(Rt, NaN => 0)

	return Rt
end

#=@benchmark DA_reflslab(0:0.01:10,[0.1,10.0], 1.0,1.0,1.0,10.0)
BenchmarkTools.Trial:
  memory estimate:  369.83 KiB
  allocs estimate:  18
  --------------
  minimum time:     722.578 μs (0.00% GC)
  median time:      811.696 μs (0.00% GC)
  mean time:        856.330 μs (2.28% GC)
  maximum time:     4.086 ms (73.78% GC)
  --------------
  samples:          5835
  evals/sample:     1
=#

function DA_reflslab1(t, β::Array{Float64,1}, ρ::Float64,ndet::Float64, nmed::Float64, s::Float64)
	#s is slab thickness

	  n::Float64 = nmed/ndet
      μa::Float64 = β[1]
      μsp::Float64 = β[2]
      D::Float64 = 1/3μsp
      ν::Float64 = 29.9792345/nmed
      xs::UnitRange{Int64} = -10:10

      Rt1 = Array{Float64}(undef, length(t))
	  Rt2 = zeros(Float64, length(t))

		@inbounds begin

			if n == 1.0
				A= 1.0
			elseif n > 1.0
				A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 		5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
			else
				A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
			end

			zs::Float64 = 1/μsp
			ze::Float64 = 2A*D

		    Rt1 = @. -exp(-(ρ^2/(4D*ν*t)) - μa*ν*t)/(2*(4π*D*ν)^(3/2)*t^(5/2))

				Threads.@threads for n in eachindex(t)
					for m in xs

						z3m = -2m.*s .- 4m.*ze .- zs
						z4m = -2m.*s .- (4m .- 2).*ze .+ zs

			 			Rt2[n] += z3m*exp(-(z3m^2 / (4D*ν*t[n]))) - z4m*exp(-(z4m^2 / (4D*ν*t[n])))
					end

				end
			end

	return replace!(Rt1.*Rt2, NaN => 0)

end

#=
@benchmark DA_reflslab1(0:0.01:10,[0.1,10.0], 1.0,1.0,1.0,10.0)
BenchmarkTools.Trial:
  memory estimate:  34.98 KiB
  allocs estimate:  26
  --------------
  minimum time:     289.363 μs (0.00% GC)
  median time:      318.652 μs (0.00% GC)
  mean time:        433.391 μs (0.81% GC)
  maximum time:     10.766 ms (57.23% GC)
  --------------
  samples:          10000
  evals/sample:     1

  =#

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
function creatematerialstable!(data::Dict{Symbol, Any})
  push!(data, :materials=>Dict{Symbol, Dict{Symbol, Float64}}())
  #
  # Insert default values for :steel
  #
  merge!(data[:materials], Dict(:steel=>Dict{Symbol, Float64}()))
  push!(data[:materials][:steel], :sg=>0.283)
  push!(data[:materials][:steel], :ez => 29.0e6)
  push!(data[:materials][:steel], :er => 29.0e6)
  push!(data[:materials][:steel], :g=>11.5e6)

  merge!(data[:materials], Dict(:monel=>Dict{Symbol, Float64}()))
  push!(data[:materials][:monel], :sg=>0.319)
  push!(data[:materials][:monel], :ez => 26.0e6)
  push!(data[:materials][:monel], :er => 26.0e6)
  push!(data[:materials][:monel], :g=>9.9e6)

  merge!(data[:materials], Dict(:titanium=>Dict{Symbol, Float64}()))
  push!(data[:materials][:titanium], :sg=>0.164)
  push!(data[:materials][:titanium], :ez => 15.5e6)
  push!(data[:materials][:titanium], :er => 15.5e6)
  push!(data[:materials][:titanium], :g=>5.9e6)

  merge!(data[:materials], Dict(:aluminium=>Dict{Symbol, Float64}()))
  push!(data[:materials][:aluminium], :sg=>0.1)
  push!(data[:materials][:aluminium], :ez => 10.3e6)
  push!(data[:materials][:aluminium], :er => 10.3e6)
  push!(data[:materials][:aluminium], :g=>3.9e6)
end