text:
### Common indexing functions

# Note: each array view type should implement
#
#   - uget(a, i):       retrieve the i-th element w.r.t. the underlying base
#   - uset!(a, v, i):   set the value of the i-th element w.r.t. the underlying base
#

idx(i::Integer) = convert(Int, i)::Int
idx(i::Int) = i

### uindex (for non-contiguous views)

xoffset0(s::NTuple{1,Int}, i1::Int) = (i1-1) * s[1]
xoffset0(s::NTuple{2,Int}, i1::Int, i2::Int) = (i1-1) * s[1] + (i2-1) * s[2]
xoffset0(s::NTuple{3,Int}, i1::Int, i2::Int, i3::Int) =
    (i1-1) * s[1] + (i2-1) * s[2] + (i3-1) * s[3]
xoffset0(s::NTuple{4,Int}, i1::Int, i2::Int, i3::Int, i4::Int) =
    (i1-1) * s[1] + (i2-1) * s[2] + (i3-1) * s[3] + (i4-1) * s[4]
xoffset0(s::NTuple{5,Int}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) =
    (i1-1) * s[1] + (i2-1) * s[2] + (i3-1) * s[3] + (i4-1) * s[4] + (i5-1) * s[5]

# for cases with s[1] == 1
#
xoffset(s::NTuple{1,Int}, i1::Int) = (i1-1)
xoffset(s::NTuple{2,Int}, i1::Int, i2::Int) = (i1-1) + (i2-1) * s[2]
xoffset(s::NTuple{3,Int}, i1::Int, i2::Int, i3::Int) =
    (i1-1) + (i2-1) * s[2] + (i3-1) * s[3]
xoffset(s::NTuple{4,Int}, i1::Int, i2::Int, i3::Int, i4::Int) =
    (i1-1) + (i2-1) * s[2] + (i3-1) * s[3] + (i4-1) * s[4]
xoffset(s::NTuple{5,Int}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) =
    (i1-1) + (i2-1) * s[2] + (i3-1) * s[3] + (i4-1) * s[4] + (i5-1) * s[5]

# 1D

uindex(a::NonContViews{T,1}, i1::Int) where {T} = 1 + xoffset(strides(a), i1)
uindex(a::NonContViews{T,1,0}, i1::Int) where {T} = 1 + xoffset0(strides(a), i1)
uindex(a::NonContViews{T,1}, i1::Int, i2::Int) where {T} = uindex(a, i1)
uindex(a::NonContViews{T,1}, i1::Int, i2::Int, i3::Int) where {T} = uindex(a, i1)
uindex(a::NonContViews{T,1}, i1::Int, i2::Int, i3::Int, i4::Int) where {T} = uindex(a, i1)
uindex(a::NonContViews{T,1}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} = uindex(a, i1)

# 2D

uindex(a::NonContViews{T,2}, i1::Int) where {T} = ((j1, j2) = ind2sub(size(a), i1); uindex(a, j1, j2))

uindex(a::NonContViews{T,2}, i1::Int, i2::Int) where {T} = 1 + xoffset(strides(a), i1, i2)
uindex(a::NonContViews{T,2,0}, i1::Int, i2::Int) where {T} = 1 + xoffset0(strides(a), i1, i2)

uindex(a::NonContViews{T,2}, i1::Int, i2::Int, i3::Int) where {T} = uindex(a, i1, i2)
uindex(a::NonContViews{T,2}, i1::Int, i2::Int, i3::Int, i4::Int) where {T} = uindex(a, i1, i2)
uindex(a::NonContViews{T,2}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} = uindex(a, i1, i2)

# 3D

uindex(a::NonContViews{T,3}, i1::Int) where {T} = ((j1, j2, j3) = ind2sub(size(a), i1); uindex(a, j1, j2, j3))
uindex(a::NonContViews{T,3}, i1::Int, i2::Int) where {T} = ((j2, j3) = ind2sub(size(a)[2:3], i2); uindex(a, i1, j2, j3))

uindex(a::NonContViews{T,3,0}, i1::Int, i2::Int, i3::Int) where {T} = 1 + xoffset0(strides(a), i1, i2, i3)
uindex(a::NonContViews{T,3}, i1::Int, i2::Int, i3::Int) where {T} = 1 + xoffset(strides(a), i1, i2, i3)

uindex(a::NonContViews{T,3}, i1::Int, i2::Int, i3::Int, i4::Int) where {T} = uindex(a, i1, i2, i3)
uindex(a::NonContViews{T,3}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} = uindex(a, i1, i2, i3)

# 4D

uindex(a::NonContViews{T,4}, i1::Int) where {T} = ((j1, j2, j3, j4) = ind2sub(size(a), i1); uindex(a, j1, j2, j3, j4))
uindex(a::NonContViews{T,4}, i1::Int, i2::Int) where {T} = ((j2, j3, j4) = ind2sub(size(a)[2:4], i2); uindex(a, i1, j2, j3, j4))
uindex(a::NonContViews{T,4}, i1::Int, i2::Int, i3::Int) where {T} = ((j3, j4) = ind2sub(size(a)[3:4], i3); uindex(a, i1, i2, j3, j4))

uindex(a::NonContViews{T,4,0}, i1::Int, i2::Int, i3::Int, i4::Int) where {T} = 1 + xoffset0(strides(a), i1, i2, i3, i4)
uindex(a::NonContViews{T,4}, i1::Int, i2::Int, i3::Int, i4::Int) where {T} = 1 + xoffset(strides(a), i1, i2, i3, i4)

uindex(a::NonContViews{T,4}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} = uindex(a, i1, i2, i3, i4)

# 5D

uindex(a::NonContViews{T,5}, i1::Int) where {T} =
    ((j1, j2, j3, j4, j5) = ind2sub(size(a), i1); uindex(a, j1, j2, j3, j4, j5))
uindex(a::NonContViews{T,5}, i1::Int, i2::Int) where {T} =
    ((j2, j3, j4, j5) = ind2sub(size(a)[2:5], i2); uindex(a, i1, j2, j3, j4, j5))
uindex(a::NonContViews{T,5}, i1::Int, i2::Int, i3::Int) where {T} =
    ((j3, j4, j5) = ind2sub(size(a)[3:5], i3); uindex(a, i1, i2, j3, j4, j5))
uindex(a::NonContViews{T,5}, i1::Int, i2::Int, i3::Int, i4::Int) where {T} =
    ((j4, j5) = ind2sub(size(a)[4:5], i4); uindex(a, i1, i2, i3, j4, j5))

uindex(a::NonContViews{T,5,0}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} =
    1 + xoffset0(strides(a), i1, i2, i3, i4, i5)
uindex(a::NonContViews{T,5}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} =
    1 + xoffset(strides(a), i1, i2, i3, i4, i5)

### getindex

getindex(a::ArrayView) = uget(a, 1)

getindex(a::ContViews, i::Int) = uget(a, i)
getindex(a::ContViews, i1::Int, i2::Int) = uget(a, sub2ind(size(a), i1, i2))
getindex(a::ContViews, i1::Int, i2::Int, i3::Int) = uget(a, sub2ind(size(a), i1, i2, i3))
getindex(a::ContViews, i1::Int, i2::Int, i3::Int, i4::Int) = uget(a, sub2ind(size(a), i1, i2, i3, i4))
getindex(a::ContViews, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) = uget(a, sub2ind(size(a), i1, i2, i3, i4, i5))

getindex(a::NonContViews, i::Int) = uget(a, uindex(a, i))
getindex(a::NonContViews, i1::Int, i2::Int) = uget(a, uindex(a, i1, i2))
getindex(a::NonContViews, i1::Int, i2::Int, i3::Int) = uget(a, uindex(a, i1, i2, i3))
getindex(a::NonContViews, i1::Int, i2::Int, i3::Int, i4::Int) = uget(a, uindex(a, i1, i2, i3, i4))
getindex(a::NonContViews, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) = uget(a, uindex(a, i1, i2, i3, i4, i5))

getindex(a::StridedArrayView, i::Integer) = getindex(a, idx(i))
getindex(a::StridedArrayView, i1::Integer, i2::Integer) = getindex(a, idx(i1), idx(i2))
getindex(a::StridedArrayView, i1::Integer, i2::Integer, i3::Integer) = getindex(a, idx(i1), idx(i2), idx(i3))
getindex(a::StridedArrayView, i1::Integer, i2::Integer, i3::Integer, i4::Integer) = getindex(a, idx(i1), idx(i2), idx(i3), idx(i4))
getindex(a::StridedArrayView, i1::Integer, i2::Integer, i3::Integer, i4::Integer, i5::Integer) = getindex(a, idx(i1), idx(i2), idx(i3), idx(i4), idx(i5))

### setindex!

setindex!(a::ContViews{T}, v, i::Int) where {T} = uset!(a, convert(T, v), i)
setindex!(a::ContViews{T}, v, i1::Int, i2::Int) where {T} = uset!(a, convert(T, v), sub2ind(size(a), i1, i2))
setindex!(a::ContViews{T}, v, i1::Int, i2::Int, i3::Int) where {T} = uset!(a, convert(T, v), sub2ind(size(a), i1, i2, i3))
setindex!(a::ContViews{T}, v, i1::Int, i2::Int, i3::Int, i4::Int) where {T} = uset!(a, convert(T, v), sub2ind(size(a), i1, i2, i3, i4))
setindex!(a::ContViews{T}, v, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} = uset!(a, convert(T, v), sub2ind(size(a), i1, i2, i3, i4, i5))

setindex!(a::NonContViews{T}, v, i::Int) where {T} = uset!(a, convert(T, v), uindex(a, i))
setindex!(a::NonContViews{T}, v, i1::Int, i2::Int) where {T} = uset!(a, convert(T, v), uindex(a, i1, i2))
setindex!(a::NonContViews{T}, v, i1::Int, i2::Int, i3::Int) where {T} = uset!(a, convert(T, v), uindex(a, i1, i2, i3))
setindex!(a::NonContViews{T}, v, i1::Int, i2::Int, i3::Int, i4::Int) where {T} = uset!(a, convert(T, v), uindex(a, i1, i2, i3, i4))
setindex!(a::NonContViews{T}, v, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) where {T} = uset!(a, convert(T, v), uindex(a, i1, i2, i3, i4, i5))

setindex!(a::StridedArrayView, v, i::Integer) = setindex!(a, v, idx(i))
setindex!(a::StridedArrayView, v, i1::Integer, i2::Integer) = setindex!(a, v, idx(i1), idx(i2))
setindex!(a::StridedArrayView, v, i1::Integer, i2::Integer, i3::Integer) = setindex!(a, v, idx(i1), idx(i2), idx(i3))
setindex!(a::StridedArrayView, v, i1::Integer, i2::Integer, i3::Integer, i4::Integer) = setindex!(a, v, idx(i1), idx(i2), idx(i3), idx(i4))
setindex!(a::StridedArrayView, v, i1::Integer, i2::Integer, i3::Integer, i4::Integer, i5::Integer) = setindex!(a, v, idx(i1), idx(i2), idx(i3), idx(i4), idx(i5))


text:
module MultiFieldFESpacesTests

using BlockArrays
using FillArrays
using Gridap.Arrays
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Integration
using Gridap.CellData
using Test

using Gridap.MultiField

order = 2

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

trian = get_triangulation(model)
degree = order
quad = CellQuadrature(trian,degree)
q = get_coordinates(quad)

ref_style = [:reference,:physical]

for ref_st in ref_style

  V = TestFESpace(model=model,order=order,reffe=:Lagrangian,conformity=:H1,valuetype=Float64,dof_space=ref_st)
  Q = TestFESpace(model=model,order=order-1,reffe=:Lagrangian,conformity=:L2,valuetype=Float64,dof_space=ref_st)

  U = TrialFESpace(V)
  P = TrialFESpace(Q)

  multi_field_style = ConsecutiveMultiFieldStyle()

  Y = MultiFieldFESpace([V,Q],multi_field_style)
  X = MultiFieldFESpace([U,P],multi_field_style)

  cell_axes = get_cell_axes(Y)
  @test isa(cell_axes[1][1],BlockedUnitRange)

  cell_axes = get_cell_axes(X)
  @test isa(cell_axes[1][1],BlockedUnitRange)

  @test num_free_dofs(X) == num_free_dofs(U) + num_free_dofs(P)
  @test num_free_dofs(X) == num_free_dofs(Y)

  dy = get_cell_basis(Y)
  @test is_test(dy)
  @test is_a_fe_cell_basis(dy)
  dv, dq = dy
  @test is_a_fe_cell_basis(dv)
  @test is_a_fe_cell_basis(dq)

  dx = get_cell_basis(X)
  @test is_a_fe_cell_basis(dx)
  @test is_trial(dx)
  du, dp = dx
  @test is_a_fe_cell_basis(du)
  @test is_a_fe_cell_basis(dp)

  cellmat = integrate(dv*du,trian,quad)
  cellvec = integrate(dv*2,trian,quad)
  cellids = get_cell_id(trian)
  cellmatvec = pair_arrays(cellmat,cellvec)
  @test isa(cellmat, VectorOfBlockArrayCoo)
  @test is_nonzero_block(cellmat,1,1)
  @test is_zero_block(cellmat,1,2)
  @test isa(cellvec, VectorOfBlockArrayCoo)
  @test is_nonzero_block(cellvec,1)
  @test is_zero_block(cellvec,2)

  matvecdata = (cellmatvec,cellids,cellids)
  matdata = (cellmat,cellids,cellids)
  vecdata = (cellvec,cellids)

  free_values = rand(num_free_dofs(X))
  xh = FEFunction(X,free_values)
  test_fe_function(xh)
  @test is_a_fe_function(xh)
  uh, ph = xh
  @test is_a_fe_function(uh)
  @test is_a_fe_function(ph)

  cell_isconstr = get_cell_isconstrained(X)
  @test cell_isconstr == Fill(false,num_cells(model))

  cell_constr = get_cell_constraints(X)
  @test isa(cell_constr,VectorOfBlockArrayCoo)

  cell_dof_ids = get_cell_dofs(X)
  @test isa(cell_dof_ids,VectorOfBlockArrayCoo)

  cf = CellField(X,get_cell_dofs(X))
  @test isa(cf,MultiFieldCellField)

  test_fe_space(X,matvecdata,matdata,vecdata)
  test_fe_space(Y,matvecdata,matdata,vecdata)

  #using Gridap.Visualization
  #writevtk(trian,"trian";nsubcells=30,cellfields=["uh" => uh, "ph"=> ph])

  f(x) = sin(4*pi*(x[1]-x[2]^2))+1
  fh = interpolate([f,f],X)
  fh = interpolate_everywhere([f,f],X)
  fh = interpolate_dirichlet([f,f],X)
end

end # module

ext: jl


text:
"""
    discretize(cont_sys, δ; [approx_model], [pade_expm], [lazy_expm], [lazy_sih])

Discretize a continuous system of ODEs with nondeterministic inputs.

## Input

- `cont_sys`     -- continuous system
- `δ`            -- step size
- `approx_model` -- the method to compute the approximation model for the
                    discretization, among:

    - `forward`    -- use forward-time interpolation
    - `backward`   -- use backward-time interpolation
    - `firstorder` -- use first order approximation of the ODE
    - `nobloating` -- do not bloat the initial states
                      (use for discrete-time reachability)

- `pade_expm`    -- (optional, default = `false`) if true, use Pade approximant
                    method to compute matrix exponentials of sparse matrices;
                    otherwise use Julia's buil-in `expm`
- `lazy_expm`    -- (optional, default = `false`) if true, compute the matrix
                    exponential in a lazy way (suitable for very large systems)
- `lazy_sih`     -- (optional, default = `true`) if true, compute the
                    symmetric interval hull in a lazy way (suitable if only a
                    few dimensions are of interest)

## Output

A discrete system.

## Notes

This function applies an approximation model to transform a continuous affine
system into a discrete affine system.
This transformation allows to do dense time reachability, i.e. such that the
trajectories of the given continuous system are included in the computed
flowpipe of the discretized system.
For discrete-time reachability, use `approx_model="nobloating"`.
"""
function discretize(cont_sys::InitialValueProblem{<:AbstractContinuousSystem},
                    δ::Float64;
                    approx_model::String="forward",
                    pade_expm::Bool=false,
                    lazy_expm::Bool=false,
                    lazy_sih::Bool=true)::InitialValueProblem{<:AbstractDiscreteSystem}

    if approx_model in ["forward", "backward"]
        return discr_bloat_interpolation(cont_sys, δ, approx_model, pade_expm,
                                         lazy_expm, lazy_sih)
    elseif approx_model == "firstorder"
        return discr_bloat_firstorder(cont_sys, δ)
    elseif approx_model == "nobloating"
        return discr_no_bloat(cont_sys, δ, pade_expm, lazy_expm)
    else
        error("The approximation model is invalid")
    end
end

"""
    bloat_firstorder(cont_sys, δ)

Compute bloating factors using first order approximation.

## Input

- `cont_sys` -- a continuous affine system
- `δ`        -- step size

## Notes

In this algorithm, the infinity norm is used.
See also: `discr_bloat_interpolation` for more accurate (less conservative)
bounds.

## Algorithm

This uses a first order approximation of the ODE, and matrix norm upper bounds,
see Le Guernic, C., & Girard, A., 2010, *Reachability analysis of linear systems
using support functions. Nonlinear Analysis: Hybrid Systems, 4(2), 250-262.*
"""
function discr_bloat_firstorder(cont_sys::InitialValueProblem{<:AbstractContinuousSystem},
                                δ::Float64)

    A, X0 = cont_sys.s.A, cont_sys.x0
    Anorm = norm(full(A), Inf)
    ϕ = expm(full(A))
    RX0 = norm(X0, Inf)

    if inputdim(cont_sys) == 0
        # linear case
        α = (exp(δ*Anorm) - 1. - δ*Anorm) * RX0
        Ω0 = CH(X0, ϕ * X0 + Ball2(zeros(size(ϕ, 1)), α))
        return DiscreteSystem(ϕ, Ω0)
    else
        # affine case; TODO: unify Constant and Varying input branches?
        Uset = inputset(cont_sys)
        if Uset isa ConstantInput
            U = next(Uset, 1)[1]
            RU = norm(U, Inf)
            α = (exp(δ*Anorm) - 1. - δ*Anorm)*(RX0 + RU/Anorm)
            β = (exp(δ*Anorm) - 1. - δ*Anorm)*RU/Anorm
            Ω0 = CH(X0, ϕ * X0 + δ * U + Ball2(zeros(size(ϕ, 1)), α))
            discr_U =  δ * U + Ball2(zeros(size(ϕ, 1)), β)
            return DiscreteSystem(ϕ, Ω0, discr_U)
        elseif Uset isa VaryingInput
            discr_U = Vector{LazySet}(length(Uset))
            for (i, Ui) in enumerate(Uset)
                RU = norm(Ui, Inf)
                α = (exp(δ*Anorm) - 1. - δ*Anorm)*(RX0 + RU/Anorm)
                β = (exp(δ*Anorm) - 1. - δ*Anorm)*RU/Anorm
                Ω0 = CH(X0, ϕ * X0 + δ * Ui + Ball2(zeros(size(ϕ, 1)), α))
                discr_U[i] =  δ * Ui + Ball2(zeros(size(ϕ, 1)), β)
            end
            return DiscreteSystem(ϕ, Ω0, discr_U)
        end
    end

end

"""
    discr_no_bloat(cont_sys, δ, pade_expm, lazy_expm)

Discretize a continuous system without bloating of the initial states, suitable
for discrete-time reachability.

## Input

- `cont_sys`     -- a continuous system
- `δ`            -- step size
- `pade_expm`    -- if `true`, use Pade approximant method to compute the
                    matrix exponential
- `lazy_expm`    -- if `true`, compute the matrix exponential in a lazy way
                    (suitable for very large systems)

## Output

A discrete system.

## Algorithm

The transformation implemented here is the following:

- `A -> Phi := exp(A*delta)`
- `U -> V := M*U`
- `X0 -> X0hat := X0`

where `M` corresponds to `Phi1(A, delta)` in Eq. (8) of *SpaceEx: Scalable
Verification of Hybrid Systems.*

In particular, there is no bloating, i.e. we don't bloat the initial states and
dont multiply the input by the step size δ, as required for the dense time case.
"""
function discr_no_bloat(cont_sys::InitialValueProblem{<:AbstractContinuousSystem},
                        δ::Float64,
                        pade_expm::Bool,
                        lazy_expm::Bool)

    A, X0 = cont_sys.s.A, cont_sys.x0
    n = size(A, 1)
    if lazy_expm
        ϕ = SparseMatrixExp(A * δ)
    else
        if pade_expm
            ϕ = padm(A * δ)
        else
            ϕ = expm(full(A * δ))
        end
    end

    # early return for homogeneous systems
    if cont_sys isa IVP{<:LinearContinuousSystem}
        Ω0 = X0
        return DiscreteSystem(ϕ, Ω0)
    end
    U = inputset(cont_sys)
    inputs = next_set(U, 1)

    # compute matrix to transform the inputs
    if lazy_expm
        P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n);
                             spzeros(n, 2*n) sparse(δ*I, n, n);
                             spzeros(n, 3*n)])
        Phi1Adelta = sparse(get_columns(P, (n+1):2*n)[1:n, :])
    else
        if pade_expm
            P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n);
                      spzeros(n, 2*n) sparse(δ*I, n, n);
                      spzeros(n, 3*n)])
        else
            P = expm(full([A*δ sparse(δ*I, n, n) spzeros(n, n);
                           spzeros(n, 2*n) sparse(δ*I, n, n);
                           spzeros(n, 3*n)]))
        end
        Phi1Adelta = P[1:n, (n+1):2*n]
    end

    discretized_U = Phi1Adelta * inputs

    Ω0 = X0

    if U isa ConstantInput
        return DiscreteSystem(ϕ, Ω0, discretized_U)
    else
        discretized_U = VaryingInput([Phi1Adelta * Ui for Ui in U])
        return DiscreteSystem(ϕ, Ω0, discretized_U)
    end
end

"""
    discr_bloat_interpolation(cont_sys, δ, approx_model, pade_expm, lazy_expm)

Compute bloating factors using forward or backward interpolation.

## Input

- `cs`           -- a continuous system
- `δ`            -- step size
- `approx_model` -- choose the approximation model among `"forward"` and
                    `"backward"`
- `pade_expm`    -- if true, use Pade approximant method to compute the
                    matrix exponential
- `lazy_expm`   --  if true, compute the matrix exponential in a lazy way
                    suitable for very large systems)

## Algorithm

See Frehse et al., CAV'11, *SpaceEx: Scalable Verification of Hybrid Systems*,
Lemma 3.

Note that in the unlikely case that A is invertible, the result can also
be obtained directly, as a function of the inverse of A and `e^{At} - I`.

The matrix `P` is such that: `ϕAabs = P[1:n, 1:n]`,
`Phi1Aabsdelta = P[1:n, (n+1):2*n]`, and `Phi2Aabs = P[1:n, (2*n+1):3*n]`.
"""
function discr_bloat_interpolation(cont_sys::InitialValueProblem{<:AbstractContinuousSystem},
                                   δ::Float64,
                                   approx_model::String,
                                   pade_expm::Bool,
                                   lazy_expm::Bool,
                                   lazy_sih::Bool)

    sih = lazy_sih ? SymmetricIntervalHull : symmetric_interval_hull

    A, X0 = cont_sys.s.A, cont_sys.x0
    n = size(A, 1)

    # compute matrix ϕ = exp(Aδ)
    if lazy_expm
        ϕ = SparseMatrixExp(A*δ)
    else
        if pade_expm
            ϕ = padm(A*δ)
        else
            ϕ = expm(full(A*δ))
        end
    end

    # early return for homogeneous systems
    if cont_sys isa IVP{<:LinearContinuousSystem}
         Ω0 = CH(X0, ϕ * X0)
        return DiscreteSystem(ϕ, Ω0)
    end
    U = inputset(cont_sys)
    inputs = next_set(U, 1)

    # compute the transformation matrix to bloat the initial states
    if lazy_expm
        P = SparseMatrixExp([abs.(A*δ) sparse(δ*I, n, n) spzeros(n, n);
                             spzeros(n, 2*n) sparse(δ*I, n, n);
                             spzeros(n, 3*n)])
        Phi2Aabs = sparse(get_columns(P, (2*n+1):3*n)[1:n, :])
    else
        if pade_expm
            P = padm([abs.(A*δ) sparse(δ*I, n, n) spzeros(n, n);
                      spzeros(n, 2*n) sparse(δ*I, n, n);
                      spzeros(n, 3*n)])
        else
            P = expm(full([abs.(A*δ) sparse(δ*I, n, n) spzeros(n, n);
                           spzeros(n, 2*n) sparse(δ*I, n, n);
                           spzeros(n, 3*n)]))
        end
        Phi2Aabs = P[1:n, (2*n+1):3*n]
    end

    if isa(inputs, ZeroSet)
        if approx_model == "forward" || approx_model == "backward"
            Ω0 = CH(X0, ϕ * X0 + δ * inputs)
        end
    else
        EPsi = sih(Phi2Aabs * sih(A * inputs))
        discretized_U = δ * inputs + EPsi
        if approx_model == "forward"
            EOmegaPlus = sih(Phi2Aabs * sih((A * A) * X0))
            Ω0 = CH(X0, ϕ * X0 + discretized_U + EOmegaPlus)
        elseif approx_model == "backward"
            EOmegaMinus = sih(Phi2Aabs * sih((A * A * ϕ) * X0))
            Ω0 = CH(X0, ϕ * X0 + discretized_U + EOmegaMinus)
        end
    end

    if U isa ConstantInput
        return DiscreteSystem(ϕ, Ω0, discretized_U)
    else
        discretized_U = [δ * Ui + sih(Phi2Aabs * sih(A * Ui)) for Ui in U]
        return DiscreteSystem(ϕ, Ω0, discretized_U)
    end
end


text:
using Test

using AtariAlgos

gamename = "breakout"
game = AtariEnv(gamename)
policy = RandomPolicy()
ep = Episode(game, policy)
for sars in ep
end
@info "total reward: $(ep.total_reward)"


text:
function LogCostGradient(X::Array,y::Array;λ=0, Norm=false)

   #Size of the X Parameters
    m=size(X,1)
    n=size(X,2)
    #if Normalize

    if Norm==true
       μ=mean(X,dims=1)
       σ=std(X,dims=1)
       X=(X.-μ)./σ
    else
        μ=zeros(1,n)
        σ=ones(1,n)
    end

    #Add vector of ones to X parameters
    X=[ones(m) X]

    CostFunction=function (θ)
                    J=((1/m).*((-1 .*y.*log.(Sigmoid(X*θ)).-(1 .-y).*log.(1 .-Sigmoid(X*θ)))'*ones(m))).+((0.5*λ/m).*(θ[2:end]'*θ[2:end]))
                    return J
                end
    θr=zeros(n+1)
    Gradient=function (θ)
             θr.=θ
             θr[1]=0
             g=(1/m).*(X'*(Sigmoid(X*θ).-y)) .+ ((λ./m).*θr)
        return g
            end

    return CostFunction, Gradient, μ, σ
end


text:
abstract type Growth{R,W} <: CellRule{R,W} end

# Basically we can define this around fitness
# and compute realized λ later

struct LogisticGrowth{R,W,LT,KT,AT,TT} <: CellRule{R,W}
    fitnessmap::LT
    K::KT
    α::AT
    dt::TT
end
LogisticGrowth{R,W}(;
    λ::LT = 1.5,
    K::KT = 100.0,
    α::AT = 1,
    dt::TT = 0.1,
) where {R,W,LT,KT,AT,TT} = LogisticGrowth{R,W}(λ, K, α, dt)

LogisticGrowth(
    layernames::T;
    λ::LT = 1.5,
    K::KT = 100.0,
    α::AT = 1,
    dt::TT = 0.1,
) where {T<:Vector{Symbol},LT,KT,AT,TT} = begin
    rules = Ruleset()
    for sym in layernames
        rules += LogisticGrowth{sym}(; λ = λ, K = K, α = α, dt = dt)
    end
    return rules
end

#LogisticGrowth(fitfunc::FitnessFunction)

function DynamicGrids.applyrule(data, rule::LogisticGrowth, X, index)
    w = fitness(fitnessmap(rule))

    dx = @fastmath rule.λ * X * (1.0 - (X / rule.K)^rule.α)
    return @fastmath X + dx * rule.dt
end


text:
export algebraicsolver, ReorderedSchurMultiplicationMatricesSolver

"""
    AbstractAlgebraicSolver

Solver of algebraic equations.
"""
abstract type AbstractAlgebraicSolver end

"""
    solvealgebraicequations(V::AbstractAlgebraicSet, algo::AbstractAlgebraicSolver)::Union{Nothing, Vector{<:Vector}}}

Solve the algebraic equations for which `V` is the set of solutions using the algorithm `algo`.
Returns a nullable which is `null` if `V` is not zero-dimensional and is the list of solutions otherwise.
"""
solvealgebraicequations(V::AbstractAlgebraicSet) = solvealgebraicequations(V, defaultalgebraicsolver(V))

"""
    AbstractMultiplicationMatricesAlgorithm

Algorithm computing multiplication matrices from algebraic equations.
"""
abstract type AbstractMultiplicationMatricesAlgorithm end

"""
    multiplicationmatrices(V::AbstractAlgebraicSet, algo::AbstractMultiplicationMatricesAlgorithm)::Union{Nothing, Vector{<:AbstractMatrix}}

Computing multiplication matrices from the algebraic equations for which `V` is the set of solution using the algorithm `algo`.
Returns a nullable which is `null` if `V` is not zero-dimensional and is the list of multiplication matrices otherwise.
"""
function multiplicationmatrices end

"""
    AbstractMultiplicationMatricesSolver

Solver of algebraic equations using multiplication matrices.
"""
abstract type AbstractMultiplicationMatricesSolver end

"""
    solvemultiplicationmatrices(Ms::AbstractVector{<:AbstractMatrix{T}}, algo::AbstractMultiplicationMatricesSolver)::Vector{Vector{T}} where T

Solve the algebraic equations having multiplication matrices `Ms` using the algorithm `algo`.
Returns the list of solutions.
"""
function solvemultiplicationmatrices end

struct SolverUsingMultiplicationMatrices{A<:AbstractMultiplicationMatricesAlgorithm, S<:AbstractMultiplicationMatricesSolver} <: AbstractAlgebraicSolver
    algo::A
    solver::S
end

function solvealgebraicequations(V::AbstractAlgebraicSet, solver::SolverUsingMultiplicationMatrices)
    Ms = multiplicationmatrices(V, solver.algo)
    if Ms === nothing
        nothing
    else
        solvemultiplicationmatrices(Ms, solver.solver)
    end
end

struct GröbnerBasisMultiplicationMatricesAlgorithm <: AbstractMultiplicationMatricesAlgorithm
end

function multiplicationmatrix(V::AbstractAlgebraicSet, v::AbstractVariable, B)
    M = Matrix{eltype(eltype(V))}(undef, length(B), length(B))
    for i in 1:length(B)
        p = rem(v * B[i], equalities(V))
        M[:, i] = coefficients(p, B)
    end
    M
end

function multiplicationmatrices(V::AbstractAlgebraicSet, algo::GröbnerBasisMultiplicationMatricesAlgorithm)
    vars = variables(V.I)
    iszd, B = monomialbasis(V.I, vars)
    if !iszd
        nothing
    else
        n = length(vars)
        if iszero(n)
            Matrix{eltype(eltype(T))}[]
        else
            [multiplicationmatrix(V, v, B) for v in vars]
        end
    end
end

include("schur.jl")

"""
Corless, R. M.; Gianni, P. M. & Trager, B. M. A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots Proceedings of the 1997 international symposium on Symbolic and algebraic computation, 1997, 133-140
"""
struct ReorderedSchurMultiplicationMatricesSolver{T, RNGT <: Random.AbstractRNG} <: AbstractMultiplicationMatricesSolver
    ɛ::T
    rng::RNGT
end
ReorderedSchurMultiplicationMatricesSolver(ɛ) = ReorderedSchurMultiplicationMatricesSolver(ɛ, Random.GLOBAL_RNG)
ReorderedSchurMultiplicationMatricesSolver{T}() where T = ReorderedSchurMultiplicationMatricesSolver(Base.rtoldefault(real(T)))

function solvemultiplicationmatrices(Ms::AbstractVector{<:AbstractMatrix{T}}, solver::ReorderedSchurMultiplicationMatricesSolver) where T
    λ = rand(solver.rng, length(Ms))
    λ /= sum(λ)
    _solvemultiplicationmatrices(Ms, λ, solver)
end

# Deterministic part
function _solvemultiplicationmatrices(Ms::AbstractVector{<:AbstractMatrix{T}}, λ, solver::ReorderedSchurMultiplicationMatricesSolver) where T<:Real
    @assert length(Ms) == length(λ)
    n = length(λ)
    Z, clusters = clusterordschur(sum(λ .* Ms), solver.ɛ)
    r = length(clusters)
    vals = [zeros(T, n) for k in 1:r]
    for k in 1:r
        nk = length(clusters[k])
        for j in clusters[k]
            q = Z[:, j]
            for i in 1:n
                vals[k][i] += dot(q, Ms[i] * q) / nk
            end
        end
    end
    vals
end

function algebraicsolver(algo::AbstractMultiplicationMatricesAlgorithm,
                         solver::AbstractMultiplicationMatricesSolver)
    SolverUsingMultiplicationMatrices(algo, solver)
end

defaultmultiplicationmatricesalgorithm(p) = GröbnerBasisMultiplicationMatricesAlgorithm()
defaultmultiplicationmatricessolver(::Type{T}) where T = ReorderedSchurMultiplicationMatricesSolver{T}()
defaultmultiplicationmatricessolver(::AbstractVector{PT}) where {T, PT<:APL{T}} = defaultmultiplicationmatricessolver(T)

function defaultalgebraicsolver(p)
    algebraicsolver(defaultmultiplicationmatricesalgorithm(p), defaultmultiplicationmatricessolver(p))
end
function defaultalgebraicsolver(p, algo::AbstractMultiplicationMatricesAlgorithm)
    algebraicsolver(algo, defaultmultiplicationmatricessolver(p))
end
function defaultalgebraicsolver(p, solver::AbstractMultiplicationMatricesSolver)
    algebraicsolver(defaultmultiplicationmatricesalgorithm(p), solver)
end


text:
using Celeste.SDSSIO
import Celeste.SDSSIO: SDSSBackground, RunCamcolField

function load_from_dataset(dataset)
    rcf = Celeste.RunCamcolField(4263, 5, 119)
    SDSSIO.load_field_images(dataset, rcf)
    SDSSIO.load_field_catalog(dataset, rcf)
end

@testset "sdssio" begin

@testset "sky interpolations" begin
    small_sky = [1.  2.  3.  4.;
                 5.  6.  7.  8.;
                 9. 10. 11. 12]
    sky_x = [0.1, 2.5]
    sky_y = [0.5, 2.5, 4.]
    sky = SDSSBackground(small_sky, sky_x, sky_y, ones(2))

    @test sky[1, 1] ≈ 1.0
    @test sky[2, 1] ≈ 7.0
    @test sky[1, 2] ≈ 2.5
    @test sky[2, 2] ≈ 8.5
    @test sky[1, 3] ≈ 4.0
    @test sky[2, 3] ≈ 10.0
end

@testset "test coordinates out of bounds" begin
    small_sky = [1.  2.  3.  4.;
            5.  6.  7.  8.;
            9. 10. 11. 12]
    sky_x = [-5.0, 4.0]
    sky_y = [-4.0, 5.0]
    sky = SDSSBackground(small_sky, sky_x, sky_y, ones(2))

    @test sky[1,1] == 1.0
    @test sky[1,2] == 4.0
    @test sky[2,1] == 9.0
    @test sky[2,2] == 12.0
end

@testset "read_photoobj handles missing table extensions" begin
    # get an example of a photoobj file with a missing table
    SampleData.get_sdss_catalog(6597, 4, 25)
end

@testset "SDSSDataSet variations" begin
    rcf = Celeste.RunCamcolField(4263, 5, 119)

    # #nsure relevant data is downloaded (this also tests I/O with
    # SDSSDataSet defaults)
    SampleData.get_sdss_images(4263, 5, 119)
    SampleData.get_sdss_catalog(4263, 5, 119)

    # test with slurp = true
    dataset = SDSSDataSet("data"; slurp = true)
    SDSSIO.load_field_images(dataset, rcf)
    SDSSIO.load_field_catalog(dataset, rcf)

    # test with SDSS layout
    rm("data/dr12", force = true, recursive=true)
    cd("data") do
        # Create the original SDSS directory layout for testing purposes
        mkpath("dr12/boss/photo/redux/301/4263/objcs/5")
        mkpath("dr12/boss/photoObj/frames/301/4263/5")
        mkpath("dr12/boss/photoObj/301/4263/5")
        for band in ['z','g','r','i','u']
            cp("4263/5/119/frame-$band-004263-5-0119.fits", "dr12/boss/photoObj/frames/301/4263/5/frame-$band-004263-5-0119.fits")
            cp("4263/5/119/fpM-004263-$(band)5-0119.fit", "dr12/boss/photo/redux/301/4263/objcs/5/fpM-004263-$(band)5-0119.fit")
        end
        cp("4263/5/119/psField-004263-5-0119.fit", "dr12/boss/photo/redux/301/4263/objcs/5/psField-004263-5-0119.fit")
        cp("4263/5/119/photoObj-004263-5-0119.fits", "dr12/boss/photoObj/301/4263/5/photoObj-004263-5-0119.fits")
        cp("4263/5/photoField-004263-5.fits", "dr12/boss/photoObj/301/4263/photoField-004263-5.fits")
    end

    dataset = SDSSDataSet("data/dr12"; dirlayout = :sdss)
    SDSSIO.load_field_images(dataset, rcf)
    SDSSIO.load_field_catalog(dataset, rcf)

    rm("data/dr12", recursive=true)

    worker = addprocs(1)[]
    fetch(@spawnat worker Base.require(:Celeste))

    # Test fetching over the network The basedir is relative, if the
    # worker doesn't try to fetch via the network after this, it'll
    # fail
    fetch(@spawnat worker cd("/tmp"))
    dataset = SDSSDataSet("data"; iostrategy = :masterrpc)
    remotecall_fetch(SDSSIO.load_field_images, worker, dataset, rcf)
    remotecall_fetch(SDSSIO.load_field_catalog, worker, dataset, rcf)
    rmprocs([worker])
end

end


text:
### Measurements.jl ---  Uncertainty propagation library
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file is the entry point of the package: it defines the new `Measurement'
# type, some functions to handle it within Julia and the new functions provided
# by the package and exposed to users.
#
### Code:

__precompile__()

module Measurements

# Calculus is used to calculate numerical derivatives in "@uncertain" macro.
using Calculus

using Requires

# Functions provided by this package and exposed to users
export Measurement, measurement, ±

# Define the "Derivatives" type, used inside "Measurement" type.  This should be
# a lightweight and immutable dictionary.
include("derivatives-type.jl")

##### New Type: Measurement
# Definition.  The Measurement type is composed by the following fields:
#   * val: the nominal value of the measurement
#   * err: the uncertainty, assumed to be standard deviation
#   * tag: a (hopefully) unique identifier, it is used to identify a specific
#     measurement in the list of derivatives.  NOTE: independent measurements
#     have tag > 0.  For dependent quantities the tag is 0.
#   * der: the list of derivates.  It is a lightweight dictionary, whose keys
#     are the tuples (nominal value, uncertainty, tag) of all independent
#     variables from which the object has been derived, the corresponding value
#     is the partial derivative of the object with respect to that independent
#     variable.  This dictionary is useful to trace the contribution of each
#     measurement and propagate the uncertainty in the case of functions with
#     more than one argument (in order to deal with correlation between
#     arguments).
struct Measurement{T<:AbstractFloat} <: AbstractFloat
    val::T
    err::T
    tag::UInt64
    der::Derivatives{T}
end
function Measurement(val::V, err::E, tag::UInt64,
                     der::Derivatives{D}) where {V,E,D}
    T = promote_type(V, E, D)
    return Measurement(T(val), T(err), tag, Derivatives{T}(der))
end
Measurement{T}(x::Measurement{S}) where {T,S} = convert(Measurement{T}, x)
Measurement{T}(x::S) where {T,S} = convert(Measurement{T}, x)
Measurement{T}(x::S) where {T,S<:Rational} = convert(Measurement{T}, x)

# Functions to quickly create an empty Derivatives object.
@generated empty_der1(x::Measurement{T}) where {T<:AbstractFloat} = Derivatives{T}()
@generated empty_der2(x::T) where {T<:AbstractFloat} = Derivatives{x}()

const tag_counters = UInt64[1]
function __init__()
    nthr = Base.Threads.nthreads()
    resize!(tag_counters, nthr)[:] = range(UInt64(1), step=typemax(UInt64)÷nthr, length=nthr)

    @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" include("unitful.jl")
    @require QuadGK="1fd47b50-473d-5c70-9696-f719f8f3bcdc" include("quadgk.jl")
    @require SpecialFunctions="276daf66-3868-5448-9aa4-cd146d93841b" include("special-functions.jl")
end

measurement(x::Measurement) = x
measurement(val::T) where {T<:AbstractFloat} = Measurement(val, zero(T), UInt64(0), empty_der2(val))
measurement(val::Real) = measurement(float(val))
function measurement(val::T, err::T) where {T<:AbstractFloat}
    newder = empty_der2(val)
    if iszero(err)
        Measurement{T}(val, err, UInt64(0), newder)
    else
        @inbounds tag = tag_counters[Base.Threads.threadid()] += 1
        return Measurement{T}(val, err, tag, Derivatives(newder, (val, err, tag)=>one(T)))
    end
end
measurement(val::Real, err::Real) = measurement(promote(float(val), float(err))...)
const ± = measurement

"""
    measurement(val::Real, [err::Real]) -> Measurement
    val ± err -> Measurement

Return a `Measurement` object with `val` as nominal value and `err` as
uncertainty.  `err` defaults to 0 if omitted.

The binary operator `±` is equivalent to `measurement`, so you can construct a
`Measurement` object by explicitely writing `123 ± 4`.
"""
measurement

include("conversions.jl")
include("comparisons-tests.jl")
include("utils.jl")
include("math.jl")
include("show.jl")
include("parsing.jl")
include("plot-recipes.jl")
include("misc.jl")

end # module


text:
import GPUCompiler
import GPUCompiler: FunctionSpec

struct Entry{F, TT}
    f::F
    specfunc::Ptr{Cvoid}
    func::Ptr{Cvoid}
end

# Slow ABI
function __call(entry::Entry{F, TT}, args::TT) where {F, TT}
    args = Any[args...]
    ccall(entry.func, Any, (Any, Ptr{Any}, Int32), entry.f, args, length(args))
end

(entry::Entry)(args...) = __call(entry, args)

const compiled_cache = Dict{UInt, Any}()

function jit(f::F,tt::TT=Tuple{}) where {F, TT<:Type}
    fspec = FunctionSpec(f, tt, #=kernel=# false, #=name=# nothing)
    GPUCompiler.cached_compilation(compiled_cache, _jit, _link, fspec)::Entry{F, tt}
end

function _link(@nospecialize(fspec::FunctionSpec), (llvm_mod, func_name, specfunc_name))
    # Now invoke the JIT
    jitted_mod = compile!(orc[], llvm_mod)

    specfunc_addr = addressin(orc[], jitted_mod, specfunc_name)
    specfunc_ptr  = pointer(specfunc_addr)

    func_addr = addressin(orc[], jitted_mod, func_name)
    func_ptr  = pointer(func_addr)

    if  specfunc_ptr === C_NULL || func_ptr === C_NULL
        @error "Compilation error" fspec specfunc_ptr func_ptr
    end

    return Entry{typeof(fspec.f), fspec.tt}(fspec.f, specfunc_ptr, func_ptr)
end

# actual compilation
function _jit(@nospecialize(fspec::FunctionSpec))
    llvm_specfunc, llvm_func, llvm_mod = codegen(fspec.f, fspec.tt)

    specfunc_name = LLVM.name(llvm_specfunc)
    func_name = LLVM.name(llvm_func)

    # set linkage to extern visible
    # otherwise `addressin` won't find them
    linkage!(llvm_func, LLVM.API.LLVMExternalLinkage)
    linkage!(llvm_specfunc, LLVM.API.LLVMExternalLinkage)

    run_pipeline!(llvm_mod)

    return (llvm_mod, func_name, specfunc_name)
end


text:
# Adapted from Plots.jl
mutable struct Current
    nullable_df::Union{DataFrame, Nothing}
    nullable_seed::Union{AbstractString, Nothing}
    nullable_path::Union{AbstractString, Nothing}
    nullable_file::Union{AbstractString, Nothing}
end
const CURRENT = Current(nothing, nothing, nothing, nothing)
isdfnull() = CURRENT.nullable_df === nothing
isseednull() = CURRENT.nullable_seed === nothing
ispathnull() = CURRENT.nullable_path === nothing
isfilenull() = CURRENT.nullable_file === nothing

function current_df()
    if isdfnull()
        @warn("No DataFrame loaded. Run the getdata_...(...) method first.")
    end
    CURRENT.nullable_df
end
current_df(df::DataFrame) = (CURRENT.nullable_df = df)

function current_seed()
    if isseednull()
        @warn("No seed loaded. Run the getdata_...(...[,seed]) method first.")
    end
    CURRENT.nullable_seed
end
current_seed(seed::AbstractString) = (CURRENT.nullable_seed = seed)
current_seed(seed::Int) = current_seed(string(seed))

function current_path()
    if ispathnull()
        @warn("No path loaded. Run the getdata_...(...) method first")
    end
    CURRENT.nullable_path
end
current_path(path::AbstractString) = (CURRENT.nullable_path = path)

function current_file()
    if isfilenull()
        @warn("No file loaded. Run savedata() or savetofolder() first.")
    end
    CURRENT.nullable_file
end
current_file(fullfilepath::AbstractString) = (CURRENT.nullable_file = fullfilepath)

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
module ExpansionItems

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export expansionitem

register_normal_element("q__expansion__item", context = @__MODULE__)

function expansionitem(args...;
                        wrap::Function = StippleUI.DEFAULT_WRAPPER,
                        kwargs...)

  wrap() do
    q__expansion__item(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

end


text:
using Random, CSV, Distributions, Plots; pyplot()
Random.seed!(0)

sampleData = CSV.read("../data/machine1.csv", header=false)[:,1]
n, N = length(sampleData), 10^6
alpha = 0.1

bootstrapSampleMeans = [mean(rand(sampleData, n)) for i in 1:N]
Lmean = quantile(bootstrapSampleMeans, alpha/2)
Umean = quantile(bootstrapSampleMeans, 1-alpha/2)

bootstrapSampleMedians = [median(rand(sampleData, n)) for i in 1:N]
Lmed = quantile(bootstrapSampleMedians, alpha/2)
Umed = quantile(bootstrapSampleMedians, 1-alpha/2)

println("Bootstrap confidence interval for the mean: ", (Lmean, Umean) )
println("Bootstrap confidence interval for the median: ", (Lmed, Umed) )

stephist(bootstrapSampleMeans, bins=1000, c=:blue,
    normed=true, label="Sample \nmeans")
plot!([Lmean, Lmean],[0,2], c=:black, ls=:dash, label="90% CI")
plot!([Umean, Umean],[0,2],c=:black, ls=:dash, label="",
    xlims=(52,54), ylims=(0,2), xlabel="Sample Means", ylabel="Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct AutoParallelOptions <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AutoParallelOptions(; kwargs...)
        obj = new(meta(AutoParallelOptions), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AutoParallelOptions
const __meta_AutoParallelOptions = Ref{ProtoMeta}()
function meta(::Type{AutoParallelOptions})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AutoParallelOptions)
            __meta_AutoParallelOptions[] = target = ProtoMeta(AutoParallelOptions)
            allflds = Pair{Symbol,Union{Type,String}}[:enable => Bool, :num_replicas => Int32]
            meta(target, AutoParallelOptions, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AutoParallelOptions[]
    end
end
function Base.getproperty(obj::AutoParallelOptions, name::Symbol)
    if name === :enable
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :num_replicas
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct ScopedAllocatorOptions <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ScopedAllocatorOptions(; kwargs...)
        obj = new(meta(ScopedAllocatorOptions), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ScopedAllocatorOptions
const __meta_ScopedAllocatorOptions = Ref{ProtoMeta}()
function meta(::Type{ScopedAllocatorOptions})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ScopedAllocatorOptions)
            __meta_ScopedAllocatorOptions[] = target = ProtoMeta(ScopedAllocatorOptions)
            allflds = Pair{Symbol,Union{Type,String}}[:enable_op => Base.Vector{AbstractString}]
            meta(target, ScopedAllocatorOptions, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ScopedAllocatorOptions[]
    end
end
function Base.getproperty(obj::ScopedAllocatorOptions, name::Symbol)
    if name === :enable_op
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    else
        getfield(obj, name)
    end
end

const RewriterConfig_Toggle = (;[
    Symbol("DEFAULT") => Int32(0),
    Symbol("ON") => Int32(1),
    Symbol("OFF") => Int32(2),
    Symbol("AGGRESSIVE") => Int32(3),
]...)

const RewriterConfig_CpuLayout = (;[
    Symbol("NO_CONVERSION_ON_CPU") => Int32(0),
    Symbol("NCHW_TO_NHWC") => Int32(1),
    Symbol("NHWC_TO_NCHW") => Int32(2),
]...)

const RewriterConfig_NumIterationsType = (;[
    Symbol("DEFAULT_NUM_ITERS") => Int32(0),
    Symbol("ONE") => Int32(1),
    Symbol("TWO") => Int32(2),
]...)

const RewriterConfig_MemOptType = (;[
    Symbol("DEFAULT_MEM_OPT") => Int32(0),
    Symbol("NO_MEM_OPT") => Int32(1),
    Symbol("MANUAL") => Int32(2),
    Symbol("SWAPPING_HEURISTICS") => Int32(4),
    Symbol("RECOMPUTATION_HEURISTICS") => Int32(5),
    Symbol("SCHEDULING_HEURISTICS") => Int32(6),
    Symbol("HEURISTICS") => Int32(3),
]...)

mutable struct RewriterConfig_CustomGraphOptimizer_ParameterMapEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RewriterConfig_CustomGraphOptimizer_ParameterMapEntry(; kwargs...)
        obj = new(meta(RewriterConfig_CustomGraphOptimizer_ParameterMapEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RewriterConfig_CustomGraphOptimizer_ParameterMapEntry (mapentry)
const __meta_RewriterConfig_CustomGraphOptimizer_ParameterMapEntry = Ref{ProtoMeta}()
function meta(::Type{RewriterConfig_CustomGraphOptimizer_ParameterMapEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RewriterConfig_CustomGraphOptimizer_ParameterMapEntry)
            __meta_RewriterConfig_CustomGraphOptimizer_ParameterMapEntry[] = target = ProtoMeta(RewriterConfig_CustomGraphOptimizer_ParameterMapEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => AttrValue]
            meta(target, RewriterConfig_CustomGraphOptimizer_ParameterMapEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RewriterConfig_CustomGraphOptimizer_ParameterMapEntry[]
    end
end
function Base.getproperty(obj::RewriterConfig_CustomGraphOptimizer_ParameterMapEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AttrValue
    else
        getfield(obj, name)
    end
end

mutable struct RewriterConfig_CustomGraphOptimizer <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RewriterConfig_CustomGraphOptimizer(; kwargs...)
        obj = new(meta(RewriterConfig_CustomGraphOptimizer), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RewriterConfig_CustomGraphOptimizer
const __meta_RewriterConfig_CustomGraphOptimizer = Ref{ProtoMeta}()
function meta(::Type{RewriterConfig_CustomGraphOptimizer})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RewriterConfig_CustomGraphOptimizer)
            __meta_RewriterConfig_CustomGraphOptimizer[] = target = ProtoMeta(RewriterConfig_CustomGraphOptimizer)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :parameter_map => Base.Dict{AbstractString,AttrValue}]
            meta(target, RewriterConfig_CustomGraphOptimizer, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RewriterConfig_CustomGraphOptimizer[]
    end
end
function Base.getproperty(obj::RewriterConfig_CustomGraphOptimizer, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :parameter_map
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,AttrValue}
    else
        getfield(obj, name)
    end
end

mutable struct RewriterConfig <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RewriterConfig(; kwargs...)
        obj = new(meta(RewriterConfig), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RewriterConfig
const __meta_RewriterConfig = Ref{ProtoMeta}()
function meta(::Type{RewriterConfig})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RewriterConfig)
            __meta_RewriterConfig[] = target = ProtoMeta(RewriterConfig)
            fnum = Int[50,1,3,13,14,24,7,8,9,10,11,2,15,18,22,23,25,19,12,17,26,4,6,20,5,21,16,100,200,300,301]
            allflds = Pair{Symbol,Union{Type,String}}[:cpu_layout_conversion => Int32, :layout_optimizer => Int32, :constant_folding => Int32, :shape_optimization => Int32, :remapping => Int32, :common_subgraph_elimination => Int32, :arithmetic_optimization => Int32, :dependency_optimization => Int32, :loop_optimization => Int32, :function_optimization => Int32, :debug_stripper => Int32, :disable_model_pruning => Bool, :scoped_allocator_optimization => Int32, :pin_to_host_optimization => Int32, :implementation_selector => Int32, :auto_mixed_precision => Int32, :auto_mixed_precision_mkl => Int32, :disable_meta_optimizer => Bool, :meta_optimizer_iterations => Int32, :min_graph_nodes => Int32, :experimental_disable_compressed_tensor_optimization => Bool, :memory_optimization => Int32, :memory_optimizer_target_node_name_scope => AbstractString, :meta_optimizer_timeout_ms => Int64, :auto_parallel => AutoParallelOptions, :fail_on_optimizer_errors => Bool, :scoped_allocator_opts => ScopedAllocatorOptions, :optimizers => Base.Vector{AbstractString}, :custom_optimizers => Base.Vector{RewriterConfig_CustomGraphOptimizer}, :inter_optimizer_verifier_config => VerifierConfig, :post_optimization_verifier_config => VerifierConfig]
            meta(target, RewriterConfig, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RewriterConfig[]
    end
end
function Base.getproperty(obj::RewriterConfig, name::Symbol)
    if name === :cpu_layout_conversion
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :layout_optimizer
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :constant_folding
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :shape_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :remapping
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :common_subgraph_elimination
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :arithmetic_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :dependency_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :loop_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :function_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :debug_stripper
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :disable_model_pruning
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :scoped_allocator_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :pin_to_host_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :implementation_selector
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :auto_mixed_precision
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :auto_mixed_precision_mkl
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :disable_meta_optimizer
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :meta_optimizer_iterations
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :min_graph_nodes
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :experimental_disable_compressed_tensor_optimization
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :memory_optimization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :memory_optimizer_target_node_name_scope
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :meta_optimizer_timeout_ms
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :auto_parallel
        return (obj.__protobuf_jl_internal_values[name])::AutoParallelOptions
    elseif name === :fail_on_optimizer_errors
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :scoped_allocator_opts
        return (obj.__protobuf_jl_internal_values[name])::ScopedAllocatorOptions
    elseif name === :optimizers
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :custom_optimizers
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{RewriterConfig_CustomGraphOptimizer}
    elseif name === :inter_optimizer_verifier_config
        return (obj.__protobuf_jl_internal_values[name])::VerifierConfig
    elseif name === :post_optimization_verifier_config
        return (obj.__protobuf_jl_internal_values[name])::VerifierConfig
    else
        getfield(obj, name)
    end
end

export AutoParallelOptions, ScopedAllocatorOptions, RewriterConfig_Toggle, RewriterConfig_CpuLayout, RewriterConfig_NumIterationsType, RewriterConfig_MemOptType, RewriterConfig_CustomGraphOptimizer_ParameterMapEntry, RewriterConfig_CustomGraphOptimizer, RewriterConfig
# mapentries: "RewriterConfig_CustomGraphOptimizer_ParameterMapEntry" => ("AbstractString", "AttrValue")


text:
#=
Created on 05/01/2021 18:39:18
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

A very basic implementation of Kronecker.jl
=#

struct Kronecker{T,TA,TB} <: AbstractMatrix{T}
    A::TA  # first matrix
    B::TB  # second matrix
    # constructor
    function Kronecker(A::TA, B::TB) where {TA<:AbstractMatrix,TB<:AbstractMatrix}
        T = promote_type(eltype(A), eltype(B))  # general type
        return new{T,TA,TB}(A, B)
    end
end

Base.eltype(K::Kronecker{T,TA,TB}) where {T,TA,TB} = T
Base.size(K::Kronecker) = size(K.A) .* size(K.B)

function  Base.getindex(K::Kronecker, i, j)
    A, B = K.A, K.B
    n, m = size(A)
    p, q = size(B)
    return A[(i-1)÷p+1,(j-1)÷q+1] * B[(i-1)%p+1,(j-1)%q+1]
end

# syntactic sugar
⊗ = Kronecker

A = rand(Bool, 100, 100)
B = randn(50, 50)

K = A ⊗ B
Kdense = kron(A, B)  # native julia version

Base.inv(K::Kronecker) = Kronecker(inv(K.A), inv(K.B))

# we use the @time macro to measure how long a certain execution takes

@time inv(K)
@time inv(Kdense)


text:
"""
This module provides some `FiniteStrain` types and `getstrain` methods to
calculate various strains.
"""
module FiniteStrains

export FiniteStrain,
       EulerianStrain,
       LagrangianStrain,
       NaturalStrain,
       InfinitesimalStrain,
       getstrain

abstract type FiniteStrain end

struct EulerianStrain <: FiniteStrain end
struct LagrangianStrain <: FiniteStrain end
struct NaturalStrain <: FiniteStrain end
struct InfinitesimalStrain <: FiniteStrain end

getstrain(::EulerianStrain, v0::Real, v::Real) = (cbrt(v0 / v)^2 - 1) / 2
getstrain(::LagrangianStrain, v0::Real, v::Real) = (cbrt(v / v0)^2 - 1) / 2
getstrain(::NaturalStrain, v0::Real, v::Real) = log(v / v0) / 3
getstrain(::InfinitesimalStrain, v0::Real, v::Real) = 1 - cbrt(v0 / v)

end


text:
"""
Constant mean
"""
struct ConstantMean{T}
	c::AbstractVector{T}
end
@functor ConstantMean
(cm::ConstantMean{T})(x) where {T} = zeros(T, size(x, 2)) .+ cm.c

function ConstantMean()
	ConstantMean([0.0])
end

reset(CM::ConstantMean) = ConstantMean()

function Base.show(io::IO, mean::ConstantMean)
	print(io, "ConstantMean(", mean.c, ")")
end

"""
simple neural network mean (single layer, linear)
"""
SimpleNeuralNetworkMean(n_features, n_out) = Dense(n_features, n_out)

reset(SNNM::Dense) = Dense(size(SNNM.W, 2), size(SNNM.W, 1))

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
### ARITHMETIC OPERATIONS

function mul!(output::Union{Tuple{AbstractVector{T},AbstractVector{T}},AbstractVectorOfArray{T}},sys::SaddleSystem{T,Ns,Nc},
  input::Union{Tuple{AbstractVector{T},AbstractVector{T}},AbstractVectorOfArray{T}}) where {T,Ns,Nc}

    @unpack A, B₂, B₁ᵀ, C = sys

    u,f = input
    r₁,r₂ = output
    length(u) == length(r₁) == Ns || error("Incompatible number of elements")
    length(f) == length(r₂) == Nc || error("Incompatible number of elements")

    r₁ .= A*u .+ B₁ᵀ*f
    r₂ .= B₂*u .+ C*f
    return output
end

function mul!(sol::Tuple{TU,TF},sys::SaddleSystem,rhs::Tuple{TU,TF}) where {T,Ns,Nc,TU,TF}
    u, f = sol
    r₁, r₂ = rhs
    return mul!((_unwrap_vec(u),_unwrap_vec(f)),sys,(_unwrap_vec(r₁),_unwrap_vec(r₂)))
end

function mul!(sol::ArrayPartition,sys::SaddleSystem,rhs::ArrayPartition) where {T,Ns,Nc,TU,TF}
    u, f = sol.x
    r₁, r₂ = rhs.x
    return mul!((_unwrap_vec(u),_unwrap_vec(f)),sys,(_unwrap_vec(r₁),_unwrap_vec(r₂)))
end

function (*)(sys::SaddleSystem,input::Tuple)
    u, f = input
    output = (similar(u),similar(f))
    mul!(output,sys,input)
    return output
end

# Routine for accepting vector inputs, parsing it into Ns and Nc parts
function mul!(sol::AbstractVector{T},sys::SaddleSystem{T,Ns,Nc},rhs::AbstractVector{T}) where {T,Ns,Nc}
    mul!(_split_vector(sol,Ns,Nc),sys,_split_vector(rhs,Ns,Nc))
    return sol
end

function (*)(sys::SaddleSystem{T},rhs::Union{AbstractVector{T},AbstractVectorOfArray{T},ArrayPartition{T}}) where {T}
    output = similar(rhs)
    mul!(output,sys,rhs)
    return output
end

## Multiplying tuples of saddle point systems

function mul!(sol,sys::Tuple{T1,T2},rhs) where {T1<:SaddleSystem,T2<:SaddleSystem}
   for (i,sysi) in enumerate(sys)
     mul!(sol[i],sysi,rhs[i])
   end
   sol
end

#=
function mul!(sol,sys::ArrayPartition{T},rhs) where {T<:SaddleSystem}
   for (i,sysi) in enumerate(sys.x)
     mul!(sol.x[i],sysi,rhs.x[i])
   end
   sol
end
=#

function (*)(sys::Tuple{T1,T2},rhs) where {T1<:SaddleSystem,T2<:SaddleSystem}
    sol = deepcopy.(rhs)
    mul!(sol,sys,rhs)
    return sol
end

#### Left division ####

function ldiv!(sol::Union{Tuple{AbstractVector{T},AbstractVector{T}},AbstractVectorOfArray{T}},sys::SaddleSystem{T,Ns,Nc},
              rhs::Union{Tuple{AbstractVector{T},AbstractVector{T}},AbstractVectorOfArray{T}}) where {T,Ns,Nc}
    @unpack A⁻¹, B₂, B₁ᵀ, B₂A⁻¹r₁, P, S⁻¹, _f_buf, A⁻¹B₁ᵀf = sys

    N = Ns+Nc
    u,f = sol
    r₁,r₂ = rhs
    length(u) == length(r₁) == Ns || error("Incompatible number of elements")
    length(f) == length(r₂) == Nc || error("Incompatible number of elements")

    mul!(u,A⁻¹,r₁)

    B₂A⁻¹r₁ .= B₂*u
    _f_buf .= r₂
    _f_buf .-= B₂A⁻¹r₁

    if Nc > 0
        f .= S⁻¹*_f_buf
        f .= P*f
    end
    A⁻¹B₁ᵀf .= A⁻¹*B₁ᵀ*f
    u .-= A⁻¹B₁ᵀf

    return sol
end

function ldiv!(sol::Tuple{TU,TF},sys::SaddleSystem,rhs::Tuple{TU,TF}) where {T,Ns,Nc,TU,TF}
    u, f = sol
    r₁, r₂ = rhs
    return ldiv!((_unwrap_vec(u),_unwrap_vec(f)),sys,(_unwrap_vec(r₁),_unwrap_vec(r₂)))
end

function ldiv!(sol::ArrayPartition,sys::SaddleSystem,rhs::ArrayPartition) where {T,Ns,Nc,TU,TF}
    u, f = sol.x
    r₁, r₂ = rhs.x
    return ldiv!((_unwrap_vec(u),_unwrap_vec(f)),sys,(_unwrap_vec(r₁),_unwrap_vec(r₂)))
end

function (\)(sys::SaddleSystem,rhs::Tuple) where {T,Ns,Nc}
    u, f = rhs
    sol = (similar(u),similar(f))
    ldiv!(sol,sys,rhs)
    return sol
end

# Routine for accepting vector inputs, parsing it into Ns and Nc parts
function ldiv!(sol::AbstractVector{T},sys::SaddleSystem{T,Ns,Nc},rhs::AbstractVector{T}) where {T,Ns,Nc}
    ldiv!(_split_vector(sol,Ns,Nc),sys,_split_vector(rhs,Ns,Nc))
    return sol
end

function (\)(sys::SaddleSystem{T},rhs::Union{AbstractVector{T},AbstractVectorOfArray{T},ArrayPartition{T}}) where {T}
    sol = similar(rhs)
    ldiv!(sol,sys,rhs)
    return sol
end

## Solving tuples of saddle point systems

function ldiv!(sol,sys::Tuple{T1,T2},rhs) where {T1<:SaddleSystem,T2<:SaddleSystem}
   for (i,sysi) in enumerate(sys)
     ldiv!(sol[i],sysi,rhs[i])
   end
   sol
end

#=
function ldiv!(sol,sys::ArrayPartition{T},rhs) where {M,T<:SaddleSystem}
   for (i,sysi) in enumerate(sys.x)
     ldiv!(sol.x[i],sysi,rhs.x[i])
   end
   sol
end
=#

function (\)(sys::Tuple{T1,T2},rhs) where {T1<:SaddleSystem,T2<:SaddleSystem}
    sol = deepcopy.(rhs)
    ldiv!(sol,sys,rhs)
    return sol
end

# vector -> tuple
_split_vector(x,Ns,Nc) = view(x,1:Ns), view(x,Ns+1:Ns+Nc)


text:
import StatsBase: standardize
import Distributions: scale

function standardize(x::AbstractArray)
  x = convert(Vector{Float64}, (x .- mean(x)))/std(x)
end

scale(x) = standardize(x)


text:
using NumericalMethodsforEngineers
using Test

N = 8
tol = 1.0e-10
limit = 1000

a = float(rand(1:10 * N, N, N))
x = float(rand(1:10, N))
b = a * x
x0 = ones(N)
(iters, converged, r2) = NumericalMethodsforEngineers.bicgstab(a, b, x0; tol=tol, limit=limit)

@test round.(a\b, digits=5) == round.(r2, digits=5)


text:
@testset "base.jl" begin
    @testset "PointSpace" begin
        @test @capture_out(print(PointSpace())) == @capture_err(print(stderr, PointSpace()))
        @test iszero(PointSpace()())
        @test iszero(rand(PointSpace()))
        @test all(iszero, rand(PointSpace(), 100))
        @test dim(PointSpace()) == 0
        @test codim(PointSpace()) == 1

        @test offsetframe(PointSpace())[1] == PointSpace()()
        @test isempty(offsetframe(PointSpace())[2])
        r = rand(5)
        @test offsetframe(PointSpace())[2] *ₚ r == r
        @test dimincrease(PointSpace()) == 0
    end

    @testset "Sphere" begin
        for n in 1:5
            sphere = Sphere(n)
            @test @capture_out(print(sphere)) == @capture_err(print(stderr, sphere))
            @test norm(sphere(rand(n)...)) ≈ 1
            @test norm(rand(sphere)) ≈ 1
            @test all(norm.(rand(sphere, 100)) .≈ 1)
            @test dim(sphere) == n
            @test codim(sphere) == 1
            if n == 1
                @test dimincrease(sphere) == 1
                of = offsetframe(sphere, 0.00)
                @test of[1] ≈ [1, 0]
                @test of[2] ≈ [-1 0; 0 0; 0 1]
                of = offsetframe(sphere, 0.25)
                @test of[1] ≈ [0, 1]
                @test of[2] ≈ [0 0; -1 0; 0 1]
            else
                @test dimincrease(sphere) == n+1
                _, f = offsetframe(sphere, rand(n)...)
                r = rand(5)
                @test f *ₚ r == vcat(zeros(n+1), r)
            end
        end
    end

    @testset "Ball" begin
        for n in 1:5
            ball = Ball(n)
            @test @capture_out(print(ball)) == @capture_err(print(stderr, ball))
            @test norm(ball(rand(n)...)) ≤ 1
            @test norm(rand(ball)) ≤ 1
            @test all(norm.(rand(ball, 100)) .≤ 1)
            @test dim(ball) == n
            @test codim(ball) == 0
            @test dimincrease(ball) == n
            r = rand(5)
            @test offsetframe(ball, rand(n)...)[2] *ₚ r == vcat(zeros(n), r)
        end
    end

    @testset "Cube" begin
        for n in 1:5
            cube = Cube(n)
            @test @capture_out(print(cube)) == @capture_err(print(stderr, cube))
            @test sum(cube(rand(n)...)) .≤ n
            @test sum(rand(cube)) ≤ n
            @test all(sum.(rand(cube, 100)) .≤ n)
            @test dim(cube) == n
            @test codim(cube) == 0
            @test dimincrease(cube) == n
            r = rand(5)
            @test offsetframe(cube, rand(n)...)[2] *ₚ r == vcat(zeros(n), r)
        end
    end

    @testset "ParametricManifold" begin
        circfun(x) = (cospi(2x), sinpi(2x))
        circ = ParametricManifold{1, 1}(circfun)
        @test @capture_out(print(circ)) == @capture_err(print(stderr, circ))
        @test norm(circ(rand())) ≈ 1
        @test norm(rand(circ)) ≈ 1
        @test all(norm.(rand(circ, 100)) .≈ 1)
        @test dim(circ) == 1
        @test codim(circ) == 1
        of = offsetframe(circ, 0.00)
        @test of[1] ≈ [1, 0]
        @test of[2] ≈ [-1 0; 0 0; 0 1]
        of = offsetframe(circ, 0.25)
        @test of[1] ≈ [0, 1]
        @test of[2] ≈ [0 0; -1 0; 0 1]
        @test dimincrease(circ) == 1
        r = rand(5)
        @test length(offsetframe(circ, rand())[2] *ₚ r) == 6
        @test circ ≅ Sphere(1)

        # TODO: others
    end
end


text:
module TypedTables

using NullableArrays, Compat

export DefaultKey, Cell, Column, Row, Table, DenseTable, KeyTable, DenseKeyTable

export rename, name, eltypes, permutecols, ncol, nrow, head, tail, storagetype, storagetypes

export unique!

export readtable, readtable!, writetable

export @Cell, @Row, @Column, @Table
export @col, @select, @filter, @filter!, @filter_mask

import Base.(==)

# TODO Do we need a unique table? It could keep a persistent hash. (Perhaps a multi-column key table with no value columns?)
# TODO union!, intersect! and setdiff!
# TODO dplyr funcitons (should we use same names or julia interfaces or both??)
#      - arrange (sorting)
#      - select/rename (combination of rename with indexing)
#      - distinct (currently unique and is already O(N log N))
#      - mutate (perhaps compute or computecol or something, given mutate has a meaning in julia)
#      - summarize
#      - sample_n and sample_frac (could use rand() interface I think)
#      - left_join, right_join, full_join (aka outer_join), semi_join, anti_join
#      - think about group_by (generates something similar to KeyTable of Tables?)
# TODO Load from an object with a dataframe interface
# TODO Figure out namespace issues for head and tail (maybe both are in Base for Julia 0.5?)
# TODO sub/slice for generating a sub-table (i.e. simply a table with different StorageTypes)
# TODO possibly implement subtable type (with no ability to insert/delete rows, though can setindex!) - is this only necessary for Row?
# TODO fix TableKey so it always references the length of the current table, not its parent
# TODO Somehow make sense of the key mess for joins
# TODO Unary relational operations like select, project, etc (see dplyr for inspiration)
# TODO Binary relational operations like left_join, etc.
# TODO other types of computed joins, map, do, etc
# TODO finish DenseTable
# TODO implement KeyTable ?
# TODO implement DenseKeyTable ?

# Julia 0.5 agenda:
# TODO Use generalized comprehensions in @select (user specified data structure)
#      and @filter (BitArray)
# TODO Use @pure functions where beneficial

include("Cell.jl")
include("Column.jl")
include("Row.jl")
include("Table.jl")

include("show.jl")
include("datamacros.jl")
include("setalgorithms.jl")
include("join.jl")
include("IO.jl")

end # module


text:
"Tests whether all elements of a collection are equal"
allsame(itr) = all(first(itr) .== itr)


text:
using TestPackage

TestPackage.greet()


text:
# module links
# """
#     Links: state,input -> output for each subsystem, of type DiffEqFlux.FastLayer

#     interface includes function of form h(xu, p) where xu is a concatenation of state and input, parameters p
# """
# using DiffEqFlux

# abstract type Link <: DiffEqFlux.FastLayer end

# include("affine_canonical.jl")
# export AffineFeedthrough, LinearFeedthrough, AffineCanonical, LinearCanonical, AffineResistive, LinearResistive
# export paramlength, initial_params, convert_verilogA

# end # module

using DiffEqFlux

abstract type Link <: DiffEqFlux.FastLayer end

include("affine_layer.jl")

""" Required as a FastLayer, define a paramsum field """
paramlength(f::AffineLayer) = f.paramsum

""" Required as a FastLayer """
initial_params(f::AffineLayer) = f.initial_params()

""" convert xu concatenation to x,u in 1d case """
@inline function extract_inputs(f::L, xu::AbstractArray{T,1}) where {T} where {L<:Link}
    x = view(xu, 1:f.state_dim)
    u = view(xu, f.state_dim+1:f.state_dim+f.input_dim)
    return x, u
end

""" convert xu concatenation to x,u in 2d case """
@inline function extract_inputs(f::L, xu::AbstractArray{<:T,2}) where {T} where {L<:Link}
    x = view(xu, 1:f.state_dim, :)
    u = view(xu, f.state_dim+1:f.state_dim+f.input_dim, :)
    return x, u
end

""" convert xu concatenation to x,u in 3d case """
@inline function extract_inputs(f::L, xu::AbstractArray{T,3}) where {T} where {L<:Link}
    x = view(xu, 1:f.state_dim, :, :)
    u = view(xu, f.state_dim+1:f.state_dim+f.input_dim, :, :)
    return x, u
end

""" Identity Link"""
struct IdentityLink{I} <: Link
    initial_params::I
    paramsum::Int
    function IdentityLink(initb = Flux.zeros)
        initial_params() = initb(0)
        new{typeof(initial_params)}(initial_params, 0)
    end
end

function (f::IdentityLink)(x, p)
    return x
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
@testset "windows" begin
    A = AxisArray([1.0 2.0;
                   3.0 4.0], Axis{:time}(1:2), Axis{:channel}(1:2))

    # Window application.
    # Test with triangular window; should be evaluated to give [1, 0.5] along
    # time dimension of length 2.
    testwindow(t) = 1-t
    # Can use Axis
    @test apply_window(A, Axis{:time}=>testwindow) == [1    2;
                                                       1.5  2]
    @test apply_window(A, Axis{:channel}=>testwindow) == [1    1;
                                                          3.0  2]
    # And integer dimension number
    @test apply_window!(A, 2, testwindow) == [1    1;
                                              3.0  2]
    # Test `A` is modified in place.
    @test A == [1    1;
                3.0  2]

    @test sinebell.((0:2)/3, skew=1.5, pow=2) == [0, 0.32311591129511724, 0.9807289153754679]
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
#	Princeton 3-D carbon cycle model, from Joos (1996) "An efficient and accurate representation ..."

let
	# http://en.wikipedia.org/wiki/Atmosphere_of_Earth
	# composition:		78.084% N2, 20.946% O2, 0.9340% Ar, 0.039% CO2	=>	28.957 g/mol  (28.96 if CO2 ignored)
	# dry air mass:		5.1352 * 10^18 kg  =  0.17734 * 10^18 "kmol air"
	# ideal gas => volume% = molar%, so:	1 GtC = 1/12 * 10^12 kmol C = 0.4689 * 10^-6 (kmol C / kmol air)
	global const ppm_per_GtC		= 0.4689	# ppm/GtC
	global const conc_CO2_preind	= 278.0		# ppm

	oceandata = [
		#oceanboxweight	oceanboxtime
		#0.70367		0.70177			# Joos t > 1 response
		#0.6113			0.70177			# Daniel fix year 1
		0.611253		0.812465		# Niclas optimized for timestep=0.01 & 25 year horizon
		0.24966			2.3488
		0.066485		15.281
		0.038344		65.359
		0.019439		347.55
		0.014819		999999999
	]

	dataDIC = [
		#coeffDIC		coeffDICtemp
		15.568e-1		-0.13993e-1
		7.4706e-3		-0.20207e-3
		-1.2748e-5		0.12015e-5
		2.4491e-7		-0.12639e-7
		-1.5468e-10		0.15326e-10
	]

	oceanboxweight	= oceandata[:,1]
	oceanboxtime 	= oceandata[:,2]
	coeffDIC		= dataDIC[:,1]
	coeffDICtemp	= dataDIC[:,2]

	constantDIC			= 1.722e17	# µmol*m3/(ppm*kg)
	mixedlayerdepth 	= 50.9		# m
	oceansurfacearea	= 3.55e14	# m2 		# (3.62e14)	 m2

	gasexchangerate		= 1/7.66	# 1/year
	oceantemperature 	= 17.7		# degC

	global const oceanboxlength = length(oceanboxweight)
	global const micromol_per_kg_per_GtC = constantDIC / (mixedlayerdepth*oceansurfacearea) * ppm_per_GtC # µmol/kg/GtC

	# ocean carbon, parameterizations of impulse responses
	global function oceancarbon!(s::ClimateState, p::ClimateParams)
		@unpack Concentration, NetFluxOcean, DeltaDICbox, DeltaCO2ocean, Temp = s
		@unpack timestep, oceantempfeedback = p
		# Mixed layer pulse response function, 3D model in Joos appendix A.2.4
		# The function has 5 exponential terms (+ one constant) with different weights and time constants
		# impulse response function: IRF(t) = sum(oceanboxweight .* exp.(-t./oceanboxtime))

		# [µmol/kg], The IRF convolution can be interpreted as a solution to a system of first order ODEs
		dt_DeltaDICbox = micromol_per_kg_per_GtC * NetFluxOcean * oceanboxweight - DeltaDICbox ./ oceanboxtime
		DeltaDICbox[:] += timestep * dt_DeltaDICbox

		# [µmol/kg], DIC = Dissolved Inorganic Carbon
		DeltaDIC = sum(DeltaDICbox)

		# [ppm] Nonlinear parameterization of ocean carbonate chemistry (Revelle buffer)
		DeltaCO2ocean = sum((coeffDIC[k] + coeffDICtemp[k]*oceantemperature) * DeltaDIC^k for k=1:5)
		#DeltaCO2ocean = sum((coeffDIC + coeffDICtemp*oceantemperature) .* DeltaDIC.^(1:5))

		# Air-sea flux
		NetFluxOcean = calc_NetFluxOcean(Concentration[:CO2], DeltaCO2ocean, Temp[1], oceantempfeedback)

		@pack! s = NetFluxOcean, DeltaDIC, DeltaDICbox, DeltaCO2ocean
	end

	global calc_NetFluxOcean(CO2concentration, DeltaCO2ocean, Temp_mixed, oceantempfeedback) =
		gasexchangerate / ppm_per_GtC *
			(CO2concentration - (conc_CO2_preind + DeltaCO2ocean) * exp(oceantempfeedback*0.0423*Temp_mixed))

	global function init_oceancarbon!(s::ClimateState, p::ClimateParams)
		@unpack Temp = s
		@unpack timestep, oceantempfeedback, equilibriumCO2 = p
		CO2emissions_year1 = 0.0
		DeltaDIC, DeltaCO2ocean = backgroundCO2_solve(equilibriumCO2)

		Concentration = Dict(g => 0.0 for g in GAS3)
		Concentration[:CO2] = conc_CO2_preind + ppm_per_GtC * timestep*CO2emissions_year1 + DeltaCO2ocean
		NetFluxOcean = calc_NetFluxOcean(Concentration[:CO2], DeltaCO2ocean, Temp[1], oceantempfeedback)
		DeltaDICbox = [zeros(length(oceanboxtime)-1); DeltaDIC]

		@pack! s = Concentration, NetFluxOcean, DeltaDICbox, DeltaDIC, DeltaCO2ocean
	end

	global function backgroundCO2_solve(equilibriumCO2)
		# set up background CO2 concentrations by doing a few quick iterations of Newton-Raphson
		# to solve OceanCO2(DeltaDIC) = equilibriumCO2 - conc_CO2_preind
		# (OceanCO2 and DeltaDIC will be zero if equilibriumCO2 is at default value of 278K)
		DeltaDIC = 330*(equilibriumCO2-conc_CO2_preind) / (330 + equilibriumCO2-conc_CO2_preind) # starting guess
		OceanCO2 = 0.0
		for i=1:5
			# split out the iterations to make the solve type stable (the function boundary helps)
			# both backgroundCO2_solve() and newtonraphson_iteration() need to be global for some reason
			DeltaDIC, OceanCO2 = newtonraphson_iteration(DeltaDIC, equilibriumCO2)
		end
		return DeltaDIC, OceanCO2
	end

	global function newtonraphson_iteration(DeltaDIC, equilibriumCO2)
		# DeltaDIC is x and OceanCO2 is y
		OceanCO2 = sum((coeffDIC[k] + coeffDICtemp[k]*oceantemperature) * DeltaDIC^k for k=1:5)
		derivative = sum((coeffDIC[k] + coeffDICtemp[k]*oceantemperature) * k * DeltaDIC^(k-1) for k=1:5)
		oceanCO2target = equilibriumCO2 - conc_CO2_preind
		return DeltaDIC - (OceanCO2 - oceanCO2target)/derivative,  OceanCO2
	end
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
using Test

include("heap_sort.jl")

@testset "test helpers" begin
    @testset "max_heapify!" begin
        input = [16 4 10 14 7 9 3 2 8 1]
        output = [16 14 10 8 7 9 3 2 4 1]
        heap = Heap(input)
        max_heapify!(heap, 2)
        @test heap.array == output
    end

    @testset "build_max_heap" begin
        input = [16 4 10 14 7 9 3 2 8 1]
        output = [16 14 10 8 7 9 3 2 4 1]
        heap = build_max_heap(input)
        @test heap.array == output

        input = [4 1 3 2 16 9 10 14 8 7]
        output = [16 14 10 8 7 9 3 2 4 1]
        heap = build_max_heap!(input)
        @test heap.array == output
    end
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
export ideal_class_monoid, islocally_isomorphic, isconjugate

###############################################################################
#
#  ICM / isisomorphic
#
###############################################################################

# Stefano Marseglia "Computing the ideal class monoid of an order"
@doc Markdown.doc"""
    ideal_class_monoid(R::NfAbsOrd)
      -> Vector{FacElem{NfOrdFracIdl, NfOrdFracIdlSet}}
    ideal_class_monoid(R::AlgAssAbsOrd)
      -> Vector{FacElem{AlgAssAbsOrdIdl, AlgAssAbsOrdIdlSet}}

Given an order $R$ in a number field or a finite product of number fields, this
function returns representatives of the isomorphism classes of fractional
ideals in $R$.
"""
function ideal_class_monoid(R::T) where { T <: Union{ NfAbsOrd, AlgAssAbsOrd } }
  @assert iscommutative(R)
  orders = overorders(R)
  result = Vector{FacElem{fractional_ideal_type(R)}}()
  for S in orders
    append!(result, _icm_bar(R, S))
  end
  return result
end

# Stefano Marseglia "Computing the ideal class monoid of an order", Prop. 4.1
@doc Markdown.doc"""
    islocally_isomorphic(I::NfAbsOrdIdl, J::NfAbsOrdIdl) -> Bool
    islocally_isomorphic(I::NfFracOrdIdl, J::NfFracOrdIdl) -> Bool
    islocally_isomorphic(I::AlgAssAbsOrdIdl, J::AlgAssAbsOrdIdl) -> Bool

Given two (fractional) ideals $I$ and $J$ of an order $R$ of an $Q$-étale
algebra, this function returns `true` if $I_p$ and $J_p$ are isomorphic for
all primes $p$ of $R$ and `false` otherwise.
"""
function islocally_isomorphic(I::T, J::T) where { T <: Union{ NfAbsOrdIdl, NfOrdFracIdl, AlgAssAbsOrdIdl } }
  IJ = colon(I, J)
  JI = colon(J, I)
  return one(_algebra(order(I))) in IJ*JI
end

# Stefano Marseglia "Computing the ideal class monoid of an order", Cor. 4.5
@doc Markdown.doc"""
    isisomorphic(I::NfAbsOrdIdl, J::NfAbsOrdIdl) -> Bool, nf_elem
    isisomorphic(I::NfFracOrdIdl, J::NfFracOrdIdl) -> Bool, nf_elem
    isisomorphic(I::AlgAssAbsOrdIdl, J::AlgAssAbsOrdIdl) -> Bool, AbsAlgAssElem

Given two (fractional) ideals $I$ and $J$ of an order $R$ of an $Q$-étale
algebra $A$, this function returns `true` and an element $a \in A$ such that
$I = aJ$ if such an element exists and `false` and $0$ otherwise.
"""
function isisomorphic(I::T, J::T) where { T <: Union{ NfAbsOrdIdl, NfOrdFracIdl, AlgAssAbsOrdIdl } }
  A = _algebra(order(I))
  if !islocally_isomorphic(I, J)
    return false, zero(A)
  end

  S = ring_of_multipliers(I)
  IS = extend(I, S)
  JS = extend(J, S)
  IJ = colon(IS, JS)
  IJ.order = S
  t, a = isprincipal(numerator(IJ, copy = false))
  if !t
    return false, zero(A)
  end
  return true, divexact(_elem_in_algebra(a, copy = false), A(denominator(IJ, copy = false)))
end

function ring_of_multipliers(I::NfOrdFracIdl)
  return ring_of_multipliers(numerator(I, copy = false)*denominator(I, copy = false))
end

###############################################################################
#
#  Internals
#
###############################################################################

# Computes representatives of the weak equivalence classes of fractional ideals
# of R with ring of multipliers S.
# Algorithm 2 in Marseglia: "Computing the ideal class monoid of an order"
function _wicm_bar(R::T, S::T) where { T <: Union{ NfAbsOrd, AlgAssAbsOrd } }
  K = _algebra(S)
  oneS = one(K)*S
  St = trace_dual(S)
  reps = Vector{fractional_ideal_type(R)}()
  I = St*colon(oneS, St)
  if one(K) in I
    push!(reps, _as_fractional_ideal_of_smaller_order(R, oneS))
    return reps
  end

  orders = overorders(S)
  for O in orders
    if O == S
      continue
    end

    StO = St*O
    if !(one(K) in StO*colon(one(K)*O, StO))
      continue
    end
    f = conductor(S, O)
    fO = f*O

    ideals = ideals_containing(O, fO, S)
    for I in ideals
      if colon(I, I) != oneS
        continue
      end
      I = _as_fractional_ideal_of_smaller_order(R, I)
      new_class = true
      for J in reps
        if islocally_isomorphic(I, J)
          new_class = false
          break
        end
      end
      if new_class
        push!(reps, I)
      end
    end
    break # We only need to do this with one over-order
  end
  return reps
end

# Computes the fractional ideals I in ICM(R) with (I:I) = S
# Part of Algorithm 3 in Marseglia: "Computing the ideal class monoid of an order"
function _icm_bar(R::T, S::T) where { T <: Union{ NfAbsOrd, AlgAssAbsOrd } }
  ideals = _wicm_bar(R, S)
  P, mP = picard_group(S)
  fac_elem_mon = FacElemMon(FracIdealSet(R))
  result = Vector{elem_type(fac_elem_mon)}()
  gens_of_P = [ _as_fractional_ideal_of_smaller_order(R, mP(P[i])) for i = 1:ngens(P) ]
  last_p = P()
  for (i, p) in enumerate(P)
    fac_elem = fac_elem_mon()
    for j = 1:ngens(P)
      if p.coeff[1, j] == last_p.coeff[1, j]
        # Avoid some dictionary operations
        continue
      end
      if iszero(p.coeff[1, j])
        continue
      end

      fac_elem.fac[gens_of_P[j]] = p.coeff[1, j]
    end
    last_p = p
    for I in ideals
      fac_elem2 = copy(fac_elem)
      if haskey(fac_elem2.fac, I)
        fac_elem2.fac[I] += fmpz(1)
      else
        fac_elem2.fac[I] = fmpz(1)
      end
      push!(result, fac_elem2)
    end
  end
  return result
end

# Computes all subgroups of S/a as fractional ideals of R.
function ideals_containing(S::T, a::T2, R::T) where { T <: Union{ NfAbsOrd, AlgAssAbsOrd }, T2 <: Union{ NfAbsOrdIdl, AlgAssAbsOrdIdl } }
  Q, mQ = quo(S, a, GrpAbFinGen)
  if order(Q) == 1
    return [ a ]
  end

  K = _algebra(S)
  d = degree(S)

  potential_basis = Vector{elem_type(_algebra(S))}(undef, d)
  ideals = fractional_ideal_type(R)[]

  for i = 1:mQ.offset
    potential_basis[i] = mQ.bottom_snf_basis[i]
  end

  offset = mQ.offset
  subs = subgroups(Q)

  function group_to_ideal(s::Tuple{GrpAbFinGen, GrpAbFinGenMap})
    H, HtoQ = image(s[2], false)
    for i = 1:(d - offset)
      v = HtoQ(H[i]).coeff
      if iszero(v)
        potential_basis[i + offset] = mQ.bottom_snf_basis[i + offset]
      else
        @assert ncols(v) == d - offset
        potential_basis[i + offset] = sum( v[1, j]*mQ.top_snf_basis[j + offset] for j = 1:(d - offset) )
      end
    end
    if typeof(R) <: AlgAssAbsOrd
      M = basis_matrix(potential_basis, FakeFmpqMat)
      return ideal(algebra(R), R, M)
    else
      M = basis_matrix(potential_basis, FakeFmpqMat)*basis_mat_inv(R, copy = false)
      return fractional_ideal(R, M)
    end
  end

  return ( group_to_ideal(s) for s in subs )
end

###############################################################################
#
#  Conjugacy classes of integral matrices
#
###############################################################################

function ideal_to_matrix(I::AlgAssAbsOrdIdl)
  O = order(I)
  A = algebra(O)
  a = primitive_element_via_number_fields(A)
  M = FakeFmpqMat(representation_matrix(a, :left))
  M = mul!(M, basis_matrix(I, copy = false), M)
  M = mul!(M, M, basis_mat_inv(I, copy = false))
  @assert isone(M.den)
  return M.num
end

function ideal_to_matrix(I::Union{ NfAbsOrdIdl, NfOrdFracIdl })
  O = order(I)
  K = nf(O)
  a = gen(K)
  M = FakeFmpqMat(representation_matrix(a))
  B = basis_matrix(I, copy = false)*basis_matrix(O, copy = false)
  C = inv(B)
  M = mul!(M, B, M)
  M = mul!(M, M, C)
  @assert isone(M.den)
  return M.num
end

function matrix_to_ideal(O::AlgAssAbsOrd, M::fmpz_mat)
  f = charpoly(M)
  A = algebra(O)
  @assert dim(A) == degree(f) # Actually A == Q[x]/f
  fields_and_maps = as_number_fields(A)
  result = zeros(A, dim(A))
  for (K, AtoK) in fields_and_maps
    MK = change_base_ring(K, M) - gen(K)*identity_matrix(K, dim(A))
    _, B = nullspace(MK)
    for j = 1:ncols(B)
      for i = 1:dim(A)
        result[i] += AtoK\B[i, j]
      end
    end
  end
  return ideal_from_lattice_gens(algebra(O), O, result), result
end

function matrix_to_ideal(O::NfAbsOrd, M::fmpz_mat)
  f = charpoly(M)
  K = nf(O)
  @assert K.pol == parent(K.pol)(f)
  result = zeros(K, degree(K))
  MK = change_base_ring(K, M) - gen(K)*identity_matrix(K, degree(K))
  _, B = nullspace(MK)
  for j = 1:ncols(B)
    for i = 1:degree(K)
      result[i] += B[i, j]
    end
  end
  return fractional_ideal_from_z_gens(O, result), result
end

# Stefano Marseglia "Computing the ideal class monoid of an order"
@doc Markdown.doc"""
    isconjugate(M::fmpz_mat, N::fmpz_mat) -> Bool, fmpz_mat

Returns `true` and a matrix $U$ with $M = U\cdot N\cdot U^{-1}$ if such a
matrix exists and `false` and $0$ otherwise.
The characteristic polynomial of $M$ is required to be square-free.
"""
function isconjugate(M::fmpz_mat, N::fmpz_mat)
  Zx, x = FlintZZ["x"]
  f = charpoly(Zx, M)
  if f != charpoly(Zx, N)
    return false, zero_matrix(FlintZZ, nrows(M), ncols(M))
  end

  fac = factor(f)
  fields = Vector{AnticNumberField}()
  for (g, e) in fac
    e != 1 ? error("The characteristic polynomial must be square-free") : nothing
    push!(fields, number_field(g)[1])
  end

  if length(fields) == 1
    return _isconjugate(equation_order(fields[1]), M, N)
  end

  A, AtoK = direct_product(fields)::Tuple{AlgAss{fmpq}, Vector{AbsAlgAssToNfAbsMor{AlgAss{fmpq}, elem_type(AlgAss{fmpq}), AnticNumberField, fmpq_mat}}}
  O = _equation_order(A)
  return _isconjugate(O, M, N)
end

function _isconjugate(O::Union{ NfAbsOrd, AlgAssAbsOrd }, M::fmpz_mat, N::fmpz_mat)
  I, basisI = matrix_to_ideal(O, M)
  J, basisJ = matrix_to_ideal(O, N)
  t, a = isisomorphic(J, I)
  @assert J == a*I
  if !t
    return false, zero_matrix(FlintZZ, nrows(M), ncols(M))
  end

  aBI = basis_matrix([ a*b for b in basisI ], FakeFmpqMat)
  BJ = basis_matrix(basisJ, FakeFmpqMat)
  UU = aBI*inv(BJ)
  @assert isone(UU.den)
  return true, UU.num
end


text:
abstract type Simulation end

struct ContinuousTimeSimulation{FP<:AbstractFilteringProblem, FA<:AbstractFilteringAlgorithm, DT<:Real} <: Simulation
    filt_prob::FP
    filt_algo::FA
    no_of_timesteps::Int
    dt::DT
end

function Base.show(io::IO, sim::ContinuousTimeSimulation)
    print(io, "Continuous-time simulation
    filtering problem:                      ", sim.filt_prob,"
    filtering algorithm:                    ", sim.filt_algo,"
    number of time steps:                   ", sim.no_of_timesteps,"
    size of time step:                      ", sim.dt)
end

"""
    run!(simulation)

Runs the simulation.
"""
function run!(simulation::Simulation; records = ())
    sfs = SimulationState(simulation.filt_prob, simulation.filt_algo)
    out = Array{Any, 2}(undef, length(records), simulation.no_of_timesteps)

    for t in 1:simulation.no_of_timesteps
        if t % 100 == 0
            print(".")
        end
        propagate!(sfs, simulation.filt_prob, simulation.filt_algo, simulation.dt)
        for (i,record) in enumerate(records)
            out[i,t] = deepcopy(record(sfs))
        end
    end
    return out
end

"""
    simulate!(filtering_algorithm, filtering_problem, no_of_timesteps, dt)

Runs a simulation of the hidden state, observation, and filtering algorithm for a duration of `no_of_timesteps`.
"""
function simulate!(filt_algo::AbstractFilteringAlgorithm{ContinuousTime, ContinuousTime}, filt_prob::AbstractFilteringProblem{S1, S2, ContinuousTime, ContinuousTime}, n_time, dt) where {S1, S2}
    simulation = ContinuousTimeSimulation(filt_prob, filt_algo, n_time, dt)
    run!(simulation)
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Baseline model
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, RData,  XLSX, LinearAlgebra, Statistics

dir = "X:/VIVES/1-Personal/Florian/git/IOTables/src/"
include(dir * "model/transform_data.jl") # Script with functions to import and transform raw data
include(dir * "model/price_hat.jl") # Script with function to obtain the price index
include(dir * "model/wage_hat.jl") # Script with function to obtain the wages and gross output

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

# Data specification
dir_raw = "C:/Users/u0148308/data/raw/" # location of raw data

# WIOD rev. 2013
source = "WIOD"
revision = "2013"
year = 2011 # specified year
N = 41 # number of countries
S = 35 # number of industries

# # WIOD rev. 2016
# source = "WIOD"
# revision = "2016"
# year = 2014 # specified year
# N = 44 # number of countries
# S = 56 # number of industries

# # OECD rev. 2021
# source = "OECD"
# revision = "2021"
# year = 1995 # specified year
# N = 71 # number of countries
# S = 45 # number of industries

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

Z, F, Y, F_ctry, TB_ctry, VA_ctry, VA_coeff, γ, α, π_Z, π_F = transform_data(dir_raw, source, revision, year, N, S)

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

# iteration parameters
vfactor = 0.4
tolerance = 1e-6
max_iteration = 100
iteration = 0
max_error = 1e7

# sectoral trade elasticity
θ = 5 # assumption
θ = fill(θ, N*S) # NS×1, work with country-industry elasticities

# trade costs
τ_hat_Z = ones(N*S, N*S) # NS×NS
τ_hat_F = ones(N*S, N) # NS×N

# initialize wages and price indices
w_hat = ones(N) # N×1

# # adjust trade balance, (if active => adjustments such that there is no trade deficit, i.e. counterfactual in itself)
TB_ctry_prime = copy(TB_ctry)
# TB_ctry_prime .= 0.0

while max_error > tolerance && iteration <= max_iteration

    # Price indices
    P_hat_Z, P_hat_F, π_hat_Z, π_hat_F, cost_hat = create_price_index_hat(w_hat, τ_hat_Z, τ_hat_F, VA_coeff, γ, π_Z, π_F, θ)

    # Counterfactual trade shares
    π_hat_Z = ifelse.(isinf.(π_hat_Z), 0.0, π_hat_Z) # remove Inf
    π_hat_F = ifelse.(isinf.(π_hat_F), 0.0, π_hat_F)

    global π_prime_Z = π_Z .* π_hat_Z # NS×NS
    global π_prime_F = π_F .* π_hat_F # NS×N

    # Labor market clearing
    w_hat_prev = copy(w_hat) # store last wage in case new optimization obtains negative wages

    w_hat, Y_prime, Z_prime, F_prime, ETB_ctry = create_wages_hat(w_hat, vfactor, π_prime_Z, π_prime_F, VA_ctry, TB_ctry_prime, γ, α)

    # update iteration parameters
    error = abs.(w_hat .- w_hat_prev)
    max_error = maximum(error) # update error
    iteration += 1 # update iteration count

    if minimum(w_hat) < 0.0
        w_hat[:] = copy(w_hat_prev)
        println("Outer loop: Iteration $iteration completed with error $max_error (wage negative, rerun with previous estimate)")
    else
        println("Outer loop: Iteration $iteration completed with error $max_error")
    end

end

w_hat

# Country consumer price index and real wage
P0_F_ctry = (ones(S,N) ./ α).^α # S×N
P0_F_ctry = [prod(P0_F_ctry[:,j]) for j in 1:N] # N×1, initial consumer price index

P_F_ctry = (P_hat_F ./ α).^α # S×N, counterfactual consumer price index
P_F_ctry = [prod(P_F_ctry[:,j]) for j in 1:N] # N×1

P_hat_F_ctry = P_F_ctry ./ P0_F_ctry # N×1, changes in country price index

w_hat_real = w_hat ./ P_hat_F_ctry # N×1, real wage changes

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

# Gross output
Y_ctry = [sum(Y[i:i+S-1]) for i in 1:S:N*S] # N×1
Y_prime_ctry = [sum(Y_prime[i:i+S-1]) for i in 1:S:N*S] # N×1

Y_hat_ctry = Y_prime_ctry ./ Y_ctry # N×1, nominal

# how to obtain real change since Y is made out of final and intermediate goods?
Y_hat_ctry_real = Y_hat_ctry ./ P_hat_F_ctry


text:
"""
	RandomVibrationParameters

Struct holding parameters/methods for Random Vibration Theory.

- `pf_method` is the method used for peak factor computation
	- `:DK80` (default) is Der Kiureghian (1980), building on Vanmarcke (1975)
	- `:CL56` is Cartwright Longuet-Higgins (1956)
- `dur_ex` is the model for excitation duration
	- `:BT14` (default) is the Boore & Thompson (2014) model - note that this is adpated to work with `r_ps`
- `dur_rms` is the model for rms duration
	- `:BT12` is the Boore & Thompson (2012) model
	- `:BT15` (default) is the Boore & Thompson (2015) model
- `dur_region` is the region specified for the duration model
	- `:WNA` (default) is western North America
	- `:ENA` is eastern North America
"""
struct RandomVibrationParameters
	pf_method::Symbol
	dur_ex::Symbol
	dur_rms::Symbol
	dur_region::Symbol
end

RandomVibrationParameters() = RandomVibrationParameters(:DK80, :BT14, :BT15, :WNA)
RandomVibrationParameters(pf) = RandomVibrationParameters(pf, :BT14, ((pf == :DK80) ? :BT15 : :BT12), :WNA)


text:
module CompressingSolvers
using Base: Order
include("./domains.jl")
include("./basis_functions.jl")
include("./multicolor_ordering.jl")
include("./create_problems.jl")
include("./reconstruction.jl")
include("utils.jl")

# Write your package code here.

end