text:
#=

####
#### Reads profiles needed for RRTMGP and related
#### calculations assuming a certain netCDF file layout
####

This file reads an example file containing atmospheric conditions (temperature, pressure, gas concentrations)
 and surface properties (emissivity, temperature), defined on `nlay` layers across a set of `ncol` columns subject to
 `nexp` perturbations, and returns them in data structures suitable for use in rte and rrtmpg. The input data
 are partitioned into a user-specified number of blocks.
For the moment only quantities relevant to longwave calculations are provided.

The example files comes from the Radiative Forcing MIP (https://www.earthsystemcog.org/projects/rfmip/)
 The protocol for this experiment allows for different specifications of which gases to consider:
all gases, (CO2, CH4, N2O) + {CFC11eq; CFC12eq + HFC-134eq}. Ozone is always included
The protocol does not specify the treatment of gases like CO

# Convention

Note that `ds` is used to denote an NC dataset.
=#

using RRTMGP.GasConcentrations
using RRTMGP.Utilities
using RRTMGP.Gases
using NCDatasets

"""
    read_atmos(ds, FT, I, gases_prescribed)

Read profiles for all columns
 - `p_lay` pressure (layers)
 - `t_lay` temperature (layers)
 - `p_lev` pressure (levels)
 - `t_lev` temperature (levels)
 - `col_dry` gas concentrations
"""
function read_atmos(ds, FT, I, gases_prescribed)

    ncol = ds.dim["col"]
    nlay = ds.dim["lay"]
    nlev = ds.dim["lev"]
    @assert nlev == nlay + 1

    p_lay = Array{FT}(ds["p_lay"][:])
    t_lay = Array{FT}(ds["t_lay"][:])
    p_lev = Array{FT}(ds["p_lev"][:])
    t_lev = Array{FT}(ds["t_lev"][:])

    gases_to_look_for = [
        h2o(),
        co2(),
        o3(),
        n2o(),
        co(),
        ch4(),
        o2(),
        n2(),
        ccl4(),
        cfc11(),
        cfc12(),
        cfc22(),
        hfc143a(),
        hfc125(),
        hfc23(),
        hfc32(),
        hfc134a(),
        cf4(),
        no2(),
    ]

    gases_in_database = filter(ug -> haskey(ds, "vmr_" * chem_name(ug)), gases_to_look_for)

    gas_concs = GasConcs(FT, I, gases_prescribed, ncol, nlay, length(gases_in_database))

    for eg in gases_in_database
        set_vmr!(gas_concs, eg, Array{FT}(ds["vmr_" * chem_name(eg)][:]))
    end

    # col_dry has unchanged allocation status on return if the variable isn't present in the netCDF file
    col_dry = haskey(ds, "col_dry") ? Array{FT}(ds["col_dry"][:]) : nothing

    return p_lay, t_lay, p_lev, t_lev, gas_concs, col_dry
end

"""
    is_sw(ds)

Does this file contain variables needed to do SW calculations?
"""
is_sw(ds) = haskey(ds, "solar_zenith_angle")

"""
    is_lw(ds)

Does this file contain variables needed to do LW calculations?
"""
is_lw(ds) = !is_sw(ds)

"""
    read_lw_bc(ds)

Read LW boundary conditions for all columns
"""
read_lw_bc(ds) = ds["t_sfc"][:], ds["emis_sfc"][:]

"""
    read_lw_rt(ds)

Read LW radiative transfer parameters
"""
read_lw_rt(ds) = length(size(ds["angle"]))

"""
    read_sw_bc(ds)

Read SW boundary conditions for all columns
"""
read_sw_bc(ds) = (
    ds["solar_zenith_angle"][:],
    ds["total_solar_irradiance"][:],
    ds["sfc_alb_direct"],
    ds["sfc_alb_diffuse"],
    haskey(ds, "tsi_scaling") ? ds["tsi_scaling"][:] : nothing,
)

"""
    read_spectral_disc!(ds, spectral_disc::AbstractOpticalProps)

Read spectral discretization
"""
function read_spectral_disc!(ds, spectral_disc::AbstractOpticalProps)
    @assert ds.dim["pair"] == 2
    band_lims_wvn = ds["band_lims_wvn"][:]
    band_lims_gpt = ds["band_lims_gpt"][:]
    spectral_disc.base = OpticalPropsBase("read_spectral_disc!", band_lims_wvn, band_lims_gpt)
end

"""
    read_sfc_test_file(ds)

Read surface SW albedo and LW emissivity spectra from the surface test file
"""
read_sfc_test_file(ds) = (ds["SW_albedo"][:], ds["LW_emissivity"][:])

"""
    read_direction(ds)

Which direction is up? Stored as a global attribute.
"""
read_direction(ds) = (ds.attrib["top_at_1"] == 1)

"""
    read_sources(ds)

Sources of upward and downward diffuse radiation, for each layer and at the surface
"""
read_sources(ds) = ds["source_up"][:], ds["source_dn"][:], ds["source_sfc"][:]

"""
    read_lw_Planck_sources!(ds, sources::SourceFuncLongWave{FT}) where FT

Longwave sources at layer centers; edges in two directions; surface
   Also directionality since this will be needed for solution
"""
function read_lw_Planck_sources!(ds, sources::SourceFuncLongWave{FT}) where {FT}
    ncol = ds.dim["col"]
    nlay = ds.dim["lay"]
    ngpt = ds.dim["gpt"]
    nband = ds.dim["band"]
    # Error checking
    @assert ds.dim["pair"] == 2
    @assert haskey(ds, "lay_src")

    # Spectral discretization
    band_lims_wvn = ds["band_lims_wvn"][:]
    band_lims_gpt = ds["band_lims_gpt"][:]

    sources.optical_props = OpticalPropsBase("SourceFuncLongWave", band_lims_wvn, band_lims_gpt)
    sources.τ .= Array{FT}(undef, ncol, nlay)

    sources.lay_source .= ds["lay_src"][:]
    sources.lev_source_inc .= ds["lev_src_inc"][:]
    sources.lev_source_dec .= ds["lev_src_dec"][:]
    sources.sfc_source .= ds["sfc_src"][:]
    return nothing

end

"""
    read_sw_solar_sources(ds)

Shortwave source at TOA
   Also directionality since this will be needed for solution
"""
read_sw_solar_sources(ds) = ds["toa_src"][:]

"""
    read_two_stream(ds)

Two-stream results: reflection and transmission for diffuse and direct radiation; also extinction
"""
function read_two_stream(ds)
    Rdif = ds["Rdif"][:]
    Tdif = ds["Tdif"][:]
    Rdir = haskey(ds, "Rdir") ? ds["Rdir"][:] : nothing
    Tdir = haskey(ds, "Tdir") ? ds["Tdir"][:] : nothing
    Tnoscat = haskey(ds, "Tnoscat") ? ds["Tnoscat"][:] : nothing
    return Rdif, Tdif, Rdir, Tdir, Tnoscat
end

"""
    read_gpt_fluxes(ds)

g-point fluxes
"""
read_gpt_fluxes(ds) =
    (ds["gpt_flux_up"][:], ds["gpt_flux_dn"][:], haskey(ds, "gpt_flux_dn_dir") ? ds["gpt_flux_dn_dir"][:] : nothing)

"""
    read_size(ds)

Find the size of the problem: columns, layers, perturbations (experiments)
"""
function read_size(ds)
    ncol = Int(ds.dim["site"])
    nlay = Int(ds.dim["layer"])
    nexp = Int(ds.dim["expt"])
    @assert ds.dim["level"] == nlay + 1
    return ncol, nlay, nexp
end

"""
    read_and_block_pt(ds)

Return layer and level

 - `p_lay` layer pressure dimensions (ncol, nlay, nblocks)
 - `p_lev` level pressure dimensions (ncol, nlay+1, nblocks)
 - `t_lay` layer temperature dimensions (ncol, nlay, nblocks)
 - `t_lev` level temperature dimensions (ncol, nlay+1, nblocks)
"""
function read_and_block_pt(ds)
    FT = Float64
    ncol_l, nlay_l, nexp_l = read_size(ds)
    @assert !any([ncol_l, nlay_l, nexp_l] .== 0)

    # Read p, T data; reshape to suit RRTMGP dimensions
    p_lay = Array(transpose(kron(ones(1, nexp_l), Array{FT}(ds["pres_layer"][:]))))
    p_lev = Array(transpose(kron(ones(1, nexp_l), Array{FT}(ds["pres_level"][:]))))
    t_lay = Array(transpose(reshape(Array{FT}(ds["temp_layer"][:]), nlay_l, ncol_l * nexp_l)))
    t_lev = Array(transpose(reshape(Array{FT}(ds["temp_level"][:]), nlay_l + 1, ncol_l * nexp_l)))

    return p_lay, p_lev, t_lay, t_lev
end

"""
    read_and_block_sw_bc(ds, blocksize)

Read and reshape shortwave boundary conditions

 - `surface_albedo` surface albedo
 - `total_solar_irradiance` total solar irradiance
 - `solar_zenith_angle` solar zenith angle
"""
function read_and_block_sw_bc(ds)
    FT = Float64
    ncol_l, nlay_l, nexp_l = read_size(ds)
    blocksize = ncol_l * nexp_l
    @assert !any([ncol_l, nlay_l, nexp_l] .== 0)

    # Check that output arrays are sized correctly : blocksize, nlay, (ncol * nexp)/blocksize
    surface_albedo = repeat(Array{FT}(ds["surface_albedo"][:]), nexp_l)
    total_solar_irradiance = repeat(Array{FT}(ds["total_solar_irradiance"][:]), nexp_l)
    solar_zenith_angle = repeat(Array{FT}(ds["solar_zenith_angle"][:]), nexp_l)
    return surface_albedo, total_solar_irradiance, solar_zenith_angle
end

"""
    read_and_block_lw_bc(ds, blocksize)

Read and reshape longwave boundary conditions

 - `surface_emissivity` surface emissivity
 - `surface_temperature` surface temperature
"""
function read_and_block_lw_bc(ds)
    FT = Float64
    ncol_l, nlay_l, nexp_l = read_size(ds)
    ncols = ncol_l * nexp_l
    @assert !any([ncol_l, nlay_l, nexp_l] .== 0)

    surface_emissivity = Array{FT}(reshape(repeat(ds["surface_emissivity"][:], 1, nexp_l), ncols))
    surface_temperature = Array{FT}(reshape(ds["surface_temperature"][:], ncols)) # alternate version

    return surface_emissivity, surface_temperature
end

"""
    determine_gas_names(ds, forcing_index)

Create a pair of string arrays - one containing the chemical name of each gas, used by the k-distribution, and
one containing the name as contained in the RFMIP input files - depending on the forcing scenario:

 - `forcing_index` (1 = all available greenhouse gases;
                    2 = CO2, CH4, N2O, CFC11eq
                    3 = CO2, CH4, N2O, CFC12eq, HFC-134eq
                    All scenarios use 3D values of ozone, water vapor so those aren't listed here
"""
function determine_gas_names(ds, forcing_index)

    @assert any(forcing_index .== [1, 2, 3])
    if forcing_index == 1
        names_in_kdist = convert.(AbstractGas, read_kdist_gas_names(ds))

        # elseif forcing_index == 2

        #   # Not part of the RFMIP specification, but oxygen is included because it's a major
        #   #    gas in some bands in the SW
        #   names_in_kdist = [co2(), ch4(), n2o(), o2(), cfc12(), cfc11()]
        # elseif forcing_index == 3

        #   # Not part of the RFMIP specification, but oxygen is included because it's a major
        #   #    gas in some bands in the SW
        #   names_in_kdist = [co2(), ch4(), n2o(), o2(), cfc12(), hfc134a()]
    end
    return names_in_kdist

end

"""
    read_kdist_gas_names(ds)

Read the names of the gases known to the k-distribution
"""
read_kdist_gas_names(ds) = lowercase.(strip.(String[join(ds["gas_names"][:][:, i]) for i in 1:ds.dim["absorber"]]))

"""
    read_and_block_gases_ty(ds, gas_names)

Read and reshape gas concentrations. RRTMGP requires gas concentrations to be supplied via a class
(GasConcs). Gas concentrations are set via a call to set_vmr!(gas_concs, name, values)
where `name` is nominally the chemical formula for the gas in question and `values` may be
a scalar, a 1-d profile assumed to apply to all columns, or an array of dimension (`ncol`, `nlay`).

This routine outputs a vector `nblocks` long of these types so each element of the array can be passed to
the rrtmgp gas optics calculation in turn.

This routine exploits RFMIP conventions: only water vapor and ozone vary by column within
each experiment.

Fields in the RFMIP file have a trailing _GM (global mean); some fields use a chemical formula and other
a descriptive name, so a map is provided between these.

 - `gas_names` Names used by the k-distribution/gas concentration type
 - `gas_conc_array` vector of [`GasConcs`](@ref)
"""
function read_and_block_gases_ty(ds, gas_names::Vector{AbstractGas})
    ncol_l, nlay_l, nexp_l = read_size(ds)
    ncols = ncol_l * nexp_l
    @assert !any([ncol_l, nlay_l, nexp_l] .== 0)
    FT = Float64
    I = Int
    gas_conc_array = GasConcs(FT, I, gas_names, ncols, nlay_l)

    # Experiment index for each column
    exp_num = repeat(collect(1:nexp_l)', ncol_l, 1)[:]
    # Water vapor and ozone depend on col, lay, exp: look just like other fields
    for gas in (h2o(), o3())
        gas_conc = Array{FT}(ds[rfmip_name(gas)][:])
        gas_conc_scaling = read_scaling(ds, FT, rfmip_name(gas))
        gas_conc_temp_3d = reshape(gas_conc, nlay_l, ncols) .* gas_conc_scaling
        gas_conc_temp_3d_t = convert(Array, transpose(gas_conc_temp_3d))
        set_vmr!(gas_conc_array, gas, gas_conc_temp_3d_t)
    end

    # All other gases are a function of experiment only
    for gas in gas_names

        # Skip 3D fields above, also NO2 since RFMIP doesn't have this
        if gas in [h2o(), o3(), no2()]
            continue
        end

        # Read the values as a function of experiment
        gas_conc_scaling = read_scaling(ds, FT, rfmip_name(gas) * "_GM")
        gas_conc_temp_1d = ds[rfmip_name(gas) * "_GM"][:] * gas_conc_scaling

        # Does every value in this block belong to the same experiment?
        if all(exp_num[2:end] .- exp_num[1] .== 0)
            # Provide a scalar value
            set_vmr!(gas_conc_array, gas, gas_conc_temp_1d[exp_num[1]])
        else
            # Create 2D field, ncols x nlay, with scalar values from each experiment
            set_vmr!(gas_conc_array, gas, repeat(gas_conc_temp_1d[exp_num[:]], 1, nlay_l))
        end
    end
    return gas_conc_array
end

read_scaling(ds, FT, varName) = parse(FT, ds[varName].attrib["units"])

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
struct ColumnSparseMatrix{T,Columns<:AbstractMatrix{T},Indices} <: AbstractMatrix{T}
    m::Int
    n::Int
    columns::Columns
    col_indices::Indices
    function ColumnSparseMatrix(m::Int, n::Int, columns::Columns, col_indices::Indices) where {Columns,Indices}
        if !issorted(col_indices)
            p = sortperm(col_indices)
            col_indices = col_indices[p]
            columns = copy(columns)
            Base.permutecols!!(columns, p)
        end
        new{eltype(columns), Columns, Indices}(m, n, columns, col_indices)
    end
end

function ColumnSparseMatrix(A::AbstractMatrix)
    m,n = size(A)
    col_indices = [j for j in 1:n
                   if !all(iszero, view(A, :, j))]
    ColumnSparseMatrix(m, n, A[:, col_indices], col_indices)
end

Base.size(A::ColumnSparseMatrix) = (A.m,A.n)

function Base.getindex(A::ColumnSparseMatrix{T}, i, j) where T
    k = findfirst(isequal(j), A.col_indices)
    if isnothing(k)
        zero(T)
    else
        A.columns[i,k]
    end
end

function Base.replace_in_print_matrix(A::ColumnSparseMatrix, i::Integer, j::Integer, s::AbstractString)
    j ∈ A.col_indices ? s : Base.replace_with_centered_mark(s)
end

function LinearAlgebra.dot(A::ColumnSparseMatrix, B::AbstractMatrix)
    T = promote_type(eltype(A), eltype(B))
    v = zero(T)
    for (j,k) in enumerate(A.col_indices)
        v += dot(view(A.columns, :, j), view(B, :, k))
    end
    v
end

LinearAlgebra.dot(A::AbstractMatrix, B::ColumnSparseMatrix) =
    conj(dot(B, A))

"""
    intersect_with_indices(u, v)

Intersect `u` and `v` and return the elements present in both, as well
as their indices in `u` and `v`, respectively.

Credits: @MasonProtter
"""
function intersect_with_indices(u, v)
    out = NamedTuple{(:i, :j, :x), Tuple{Int, Int, Base.promote_eltypeof(u, v)}}[]
    foreach(eachindex(u)) do i
        @inbounds ui = u[i]
        j = findfirst(==(ui), v)
        if j !== nothing
            @inbounds vj = v[j]
            push!(out, (i=i, j=j, x=ui))
        end
    end
    out
end

function LinearAlgebra.dot(A::ColumnSparseMatrix, B::ColumnSparseMatrix)
    T = promote_type(eltype(A), eltype(B))
    v = zero(T)
    for (i,j,_) in intersect_with_indices(A.col_indices, B.col_indices)
        v += dot(view(A.columns, :, i), view(B.columns, :, j))
    end
    v
end

function Base.:(*)(A::AbstractMatrix, B::ColumnSparseMatrix)
    ma,na = size(A)
    mb,nb = size(B)
    na == mb || throw(DimensionMismatch("Sizes of A $((ma,na)) and B $((mb,nb)) not compatible ($(na) != $(mb))"))

    T = promote_type(eltype(A), eltype(B))
    columns = A*B.columns

    ColumnSparseMatrix(ma, nb, columns, B.col_indices)
end

function Base.:(*)(A::Adjoint{<:Any,<:ColumnSparseMatrix}, B::ColumnSparseMatrix)
    ma,na = size(A)
    mb,nb = size(B)
    na == mb || throw(DimensionMismatch("Sizes of A $((ma,na)) and B $((mb,nb)) not compatible ($(na) != $(mb))"))

    T = promote_type(eltype(A), eltype(B))

    I = Int[]
    J = Int[]
    V = T[]

    Ap = parent(A)
    for (k,i) in enumerate(Ap.col_indices)
        u = view(Ap.columns, :, k)
        for (l,j) in enumerate(B.col_indices)
            v = view(B.columns, :, l)
            push!(I, i)
            push!(J, j)
            push!(V, u'v)
        end
    end

    sparse(I, J, V, ma, nb)
end

Base.:(*)(α::Number, A::ColumnSparseMatrix) =
    ColumnSparseMatrix(A.m, A.n, α*A.columns, A.col_indices)

Base.:(*)(A::ColumnSparseMatrix, α::Number) = α*A

Base.:(/)(A::ColumnSparseMatrix, α::Number) = inv(α)*A

export ColumnSparseMatrix

ext: jl


text:
"""
A TensorProductOperator represents the tensor product of other operators.

struct TensorProductOperator{T} <: DictionaryOperator{T}
"""
struct TensorProductOperator{T} <: DictionaryOperator{T}
    src             ::  Dictionary
    dest            ::  Dictionary
    operators
    scratch
    src_scratch
    dest_scratch
end

# Generic functions for composite types:
components(op::TensorProductOperator) = op.operators
component(op::TensorProductOperator, j::Int) = op.operators[j]

factors(op::TensorProductOperator) = components(op)

TensorProductOperator(operators::AbstractOperator...; T=promote_type(map(eltype, operators)...)) =
    TensorProductOperator{T}(operators...)

function TensorProductOperator{T}(operators::AbstractOperator...) where {T}
    tp_src = tensorproduct(map(src, operators)...)
    tp_dest = tensorproduct(map(dest, operators)...)
    TensorProductOperator{T}(tp_src, tp_dest, operators...)
end

function TensorProductOperator{T}(tp_src::Dictionary, tp_dest::Dictionary, operators::AbstractOperator...) where {T}
    L = length(operators)
    # Scratch contains matrices of sufficient size to hold intermediate results
    # in the application of the tensor product operator.
    # Example, for L=3 scratch is a length (L-1)-tuple of matrices of size:
    # - [M1,N2,N3]
    # - [M1,M2,N3]
    # where operator J maps a set of length Nj to a set of length Mj.
    scratch_array = [ zeros(T, [length(dest(operators[k])) for k=1:j-1]..., [length(src(operators[k])) for k=j:L]...) for j=2:L]
    scratch = (scratch_array...,)

    # scr_scratch and dest_scratch are tuples of length len that contain preallocated
    # storage to hold a vector for source and destination for each operator
    src_scratch_array = [zeros(T, src(operators[j])) for j=1:L]
    src_scratch = (src_scratch_array...,)
    dest_scratch_array = [zeros(T, dest(operators[j])) for j=1:L]
    dest_scratch = (dest_scratch_array...,)
    TensorProductOperator{T}(tp_src, tp_dest, operators, scratch, src_scratch, dest_scratch)
end

function tensorproduct(ops::IdentityOperator...)
    tp_src = tensorproduct(map(src, ops)...)
    tp_dest = tensorproduct(map(dest, ops)...)
    IdentityOperator(tp_src, tp_dest)
end

# Element-wise src and dest functions
src(op::TensorProductOperator, j::Int) = src(component(op, j))
dest(op::TensorProductOperator, j::Int) = dest(component(op, j))

unsafe_wrap_operator(dict1::TensorProductDict, dict2::TensorProductDict, op::TensorProductOperator{T}) where {T} =
    TensorProductOperator{T}(dict1, dict2, map(wrap_operator, components(dict1), components(dict2), components(op))...)
unsafe_wrap_operator(dict1::Dictionary, dict2::Dictionary, op::TensorProductOperator{T}) where {T} =
    TensorProductOperator{T}(dict1, dict2, op.operators, op.scratch, op.src_scratch, op.dest_scratch)

#getindex(op::TensorProductOperator, j::Int) = component(op, j)
adjoint(op::TensorProductOperator) = TensorProductOperator(map(adjoint, components(op))...)
conj(op::TensorProductOperator) = TensorProductOperator(map(conj, components(op))...)

pinv(op::TensorProductOperator) = TensorProductOperator(map(pinv, components(op))...)

inv(op::TensorProductOperator) = TensorProductOperator(map(inv, components(op))...)

isinplace(op::TensorProductOperator) = all(map(isinplace, op.operators))
isdiag(op::TensorProductOperator) = all(map(isdiag, op.operators))

isidentity(op::TensorProductOperator) = all(map(isidentity, op.operators))

# Matrix(op::TensorProductOperator) = kron((Matrix(op) for op in reverse(components(op)))...)
sparse(op::TensorProductOperator) = kron((sparse(op) for op in reverse(components(op)))...)

unsafe_wrap_operator(src, dest, op::TensorProductOperator{T}) where T =
    TensorProductOperator{T}(src, dest, op.operators, op.scratch, op.src_scratch, op.dest_scratch)

# We reshape any incoming or outgoing coefficients into an array of the right size
function apply!(op::TensorProductOperator, coef_dest, coef_src::AbstractVector)
    @warn "Reshaping input of tensor product operator from vector to tensor" maxlog=4
    apply!(op, coef_dest, reshape(coef_src, size(src(op))))
end

function apply!(op::TensorProductOperator, coef_dest::Vector, coef_src)
    @warn "Reshaping output of tensor product operator from vector to tensor" maxlog=4
    apply!(op, reshape(coef_dest, size(dest(op))), coef_src)
end

function apply!(op::TensorProductOperator, coef_dest::Vector, coef_src::AbstractVector)
    @warn "Reshaping input and output of tensor product operator from vector to tensor" maxlog=4
    apply!(op, reshape(coef_dest, size(dest(op))), reshape(coef_src, size(src(op))))
end

apply_inplace!(op::TensorProductOperator, coef_srcdest::AbstractVector) =
    apply_inplace!(op, reshape(coef_srcdest, size(src(op))))

apply!(op::TensorProductOperator, coef_dest, coef_src) =
    apply_tensor!(op, coef_dest, coef_src, op.operators, op.scratch, op.src_scratch, op.dest_scratch)

apply_inplace!(op::TensorProductOperator, coef_srcdest) =
    apply_inplace_tensor!(op, coef_srcdest, op.operators, op.src_scratch)

function apply_tensor!(op, coef_dest, coef_src, operators::Tuple{A}, scratch, src_scratch, dest_scratch) where {A}
    println("One-element TensorProductOperators should not exist!")
    apply!(operators[1], coef_dest, coef_src)
end

function apply_inplace_tensor!(op, coef_srcdest, operators::Tuple{A}, src_scratch) where {A}
    println("One-element TensorProductOperators should not exist!")
    apply!(operators[1], coef_srcdest)
end

function apply_tensor!(op, coef_dest, coef_src, operators::Tuple{A,B}, scratch, src_scratch, dest_scratch) where {A,B}
    src1 = src_scratch[1]
    src2 = src_scratch[2]
    dest1 = dest_scratch[1]
    dest2 = dest_scratch[2]
    op1 = operators[1]
    op2 = operators[2]
    intermediate = scratch[1]

    # println("coef_src: ", size(coef_src))
    # println("coef_dest: ", size(coef_dest))
    # println("src1: ", size(src1))
    # println("src2: ", size(src2))
    # println("dest1: ", size(dest1))
    # println("dest2: ", size(dest2))
    # println("op1: ", size(op1))
    # println("op2: ", size(op2))
    # println("intermediate: ", size(intermediate))
    #
    for j in eachindex(src2)
        for i in eachindex(src1)
            src1[i] = coef_src[i,j]
        end
        apply!(op1, dest1, src1)
        for i in eachindex(dest1)
            intermediate[i,j] = dest1[i]
        end
    end
    for i in eachindex(dest1)
        for j in eachindex(src2)
            src2[j] = intermediate[i,j]
        end
        apply!(op2, dest2, src2)
        for j in eachindex(dest2)
            coef_dest[i,j] = dest2[j]
        end
    end
    coef_dest
end

function apply_inplace_tensor!(op, coef_srcdest, operators::Tuple{A,B}, src_scratch) where {A,B}
    src1 = src_scratch[1]
    src2 = src_scratch[2]
    op1 = operators[1]
    op2 = operators[2]

    for j in eachindex(src2)
        for i in eachindex(src1)
            src1[i] = coef_srcdest[i,j]
        end
        apply!(op1, src1)
        for i in eachindex(src1)
            coef_srcdest[i,j] = src1[i]
        end
    end
    for i in eachindex(src1)
        for j in eachindex(src2)
            src2[j] = coef_srcdest[i,j]
        end
        apply!(op2, src2)
        for j in eachindex(src2)
            coef_srcdest[i,j] = src2[j]
        end
    end
    coef_srcdest
end

function apply_tensor!(op, coef_dest, coef_src, operators::Tuple{A,B,C}, scratch, src_scratch, dest_scratch) where {A,B,C}
    src1 = src_scratch[1]
    src2 = src_scratch[2]
    src3 = src_scratch[3]
    dest1 = dest_scratch[1]
    dest2 = dest_scratch[2]
    dest3 = dest_scratch[3]
    op1 = operators[1]
    op2 = operators[2]
    op3 = operators[3]
    intermediate1 = scratch[1]
    intermediate2 = scratch[2]

    # println("coef_src: ", size(coef_src))
    # println("coef_dest: ", size(coef_dest))
    # println("src1: ", size(src1))
    # println("src2: ", size(src2))
    # println("dest1: ", size(dest1))
    # println("dest2: ", size(dest2))
    # println("op1: ", size(op1))
    # println("op2: ", size(op2))
    # println("intermediate: ", size(intermediate))

    for k in eachindex(src3)
        for j in eachindex(src2)
            for i in eachindex(src1)
                src1[i] = coef_src[i,j,k]
            end
            apply!(op1, dest1, src1)
            for i in eachindex(dest1)
                intermediate1[i,j,k] = dest1[i]
            end
        end
    end
    for k in eachindex(src3)
        for i in eachindex(dest1)
            for j in eachindex(src2)
                src2[j] = intermediate1[i,j,k]
            end
            apply!(op2, dest2, src2)
            for j in eachindex(dest2)
                intermediate2[i,j,k] = dest2[j]
            end
        end
    end
    for i in eachindex(dest1)
        for j in eachindex(dest2)
            for k in eachindex(src3)
                src3[k] = intermediate2[i,j,k]
            end
            apply!(op3, dest3, src3)
            for k in eachindex(dest3)
                coef_dest[i,j,k] = dest3[k]
            end
        end
    end
    coef_dest
end

function apply_inplace_tensor!(op, coef_srcdest, operators::Tuple{A,B,C}, src_scratch) where {A,B,C}
    src1 = src_scratch[1]
    src2 = src_scratch[2]
    src3 = src_scratch[3]
    op1 = operators[1]
    op2 = operators[2]
    op3 = operators[3]

    for k in eachindex(src3)
        for j in eachindex(src2)
            for i in eachindex(src1)
                src1[i] = coef_srcdest[i,j,k]
            end
            apply!(op1, src1)
            for i in eachindex(src1)
                coef_srcdest[i,j,k] = src1[i]
            end
        end
    end
    for k in eachindex(src3)
        for i in eachindex(src1)
            for j in eachindex(src2)
                src2[j] = coef_srcdest[i,j,k]
            end
            apply!(op2, src2)
            for j in eachindex(src2)
                coef_srcdest[i,j,k] = src2[j]
            end
        end
    end
    for i in eachindex(src1)
        for j in eachindex(src2)
            for k in eachindex(src3)
                src3[k] = coef_srcdest[i,j,k]
            end
            apply!(op3, src3)
            for k in eachindex(src3)
                coef_srcdest[i,j,k] = src3[k]
            end
        end
    end
    coef_srcdest
end

function stencilarray(op::TensorProductOperator)
    A = Any[]
    push!(A, component(op,1))
    for i in 2:ncomponents(op)
        push!(A, " ⊗ ")
        push!(A, component(op,i))
    end
    A
end

stencil_parentheses(op::TensorProductOperator) = true
object_parentheses(op::TensorProductOperator) = true


text:
include("trapezoid_rule.jl");

 function ordertolen(order::Int64)
    #Three dimensional tensor product basis.
    expsz=3^order;

    #Using vector of three components.
    expsz=3*expsz;

    return expsz;
end

 function evaluate_tensor_basis(i::Int64,j::Int64,k::Int64,x::Float64,y::Float64,z::Float64)
    return exp(im*(i*x+j*y+k*z));
end

 function evaluate_tensor_basis_derivative(i::Int64,j::Int64,k::Int64,d::Int64,x::Float64,y::Float64,z::Float64)
    out=exp(im*(i*x+j*y+k*z));
    is=[i,j,k];
    out=im*is[d]*out;
    return out;
end
 function evaluate_tensor_basis_derivative(i::Int64,j::Int64,k::Int64,d1::Int64,d2::Int64,x::Float64,y::Float64,z::Float64)
    out=exp(im*(i*x+j*y+k*z));
    is=[i,j,k];
    out=(im*is[d2])*(im*is[d1])*out;
    return out;
end

 function evaluate_curl_basis_gradient(order::Int64,n::Int64,x::Float64,y::Float64,z::Float64)
    ncomps=3;
    m=div(n,ncomps);
    r=n%ncomps;
    (id1,id2,id3)=ind2sub((order,order,order),m);
    out=Array(Complex128,(ncomps,ncomps));
    if r==0
        dyx=evaluate_tensor_basis_derivative(id1,id2,id3,2,1,x,y,z);
        dzx=evaluate_tensor_basis_derivative(id1,id2,id3,3,1,x,y,z);
        dyy=evaluate_tensor_basis_derivative(id1,id2,id3,2,2,x,y,z);
        dzy=evaluate_tensor_basis_derivative(id1,id2,id3,3,2,x,y,z);
        dyz=evaluate_tensor_basis_derivative(id1,id2,id3,2,3,x,y,z);
        dzz=evaluate_tensor_basis_derivative(id1,id2,id3,3,3,x,y,z);
        out[:,1]=[0.0*im,dzx,-dyx];
        out[:,2]=[0.0*im,dzy,-dyy];
        out[:,3]=[0.0*im,dzy,-dyy];
    elseif r==1
        dxx=evaluate_tensor_basis_derivative(id1,id2,id3,1,1,x,y,z);
        dzx=evaluate_tensor_basis_derivative(id1,id2,id3,3,1,x,y,z);
        dxy=evaluate_tensor_basis_derivative(id1,id2,id3,1,2,x,y,z);
        dzy=evaluate_tensor_basis_derivative(id1,id2,id3,3,2,x,y,z);
        dxz=evaluate_tensor_basis_derivative(id1,id2,id3,1,3,x,y,z);
        dzz=evaluate_tensor_basis_derivative(id1,id2,id3,3,3,x,y,z);

        out[:,1]=[-dzx,0.0*im,dxx];
        out[:,2]=[-dzy,0.0*im,dxy];
        out[:,3]=[-dzz,0.0*im,dxz];
    else
        dxx=evaluate_tensor_basis_derivative(id1,id2,id3,1,1,x,y,z);
        dyx=evaluate_tensor_basis_derivative(id1,id2,id3,2,1,x,y,z);
        dxy=evaluate_tensor_basis_derivative(id1,id2,id3,1,2,x,y,z);
        dyy=evaluate_tensor_basis_derivative(id1,id2,id3,2,2,x,y,z);
        dxz=evaluate_tensor_basis_derivative(id1,id2,id3,1,3,x,y,z);
        dyz=evaluate_tensor_basis_derivative(id1,id2,id3,2,3,x,y,z);
        out[:,1]=[dyx,-dxx,0.0*im];
        out[:,2]=[dyy,-dxy,0.0*im];
        out[:,3]=[dyz,-dxz,0.0*im];
    end
    return out;
end

 function evaluate_curl_basis(order::Int64,n::Int64,x::Float64,y::Float64,z::Float64)
    ncomps=3;
    m=div(n,ncomps);
    r=n%ncomps;
    (id1,id2,id3)=ind2sub((order,order,order),m);
    out=(0.0*im,0.0*im,0.0*im);
    if r==0
        dy=evaluate_tensor_basis_derivative(id1,id2,id3,2,x,y,z);
        dz=evaluate_tensor_basis_derivative(id1,id2,id3,3,x,y,z);
        out=(0.0*im,dz,-dy);
    elseif r==1
        dx=evaluate_tensor_basis_derivative(id1,id2,id3,1,x,y,z);
        dz=evaluate_tensor_basis_derivative(id1,id2,id3,3,x,y,z);
        out=(-dz,0.0*im,dx);
    else
        dx=evaluate_tensor_basis_derivative(id1,id2,id3,1,x,y,z);
        dy=evaluate_tensor_basis_derivative(id1,id2,id3,2,x,y,z);
        out=(dy,-dx,0.0*im);
    end
    return out;
end

 function curl_basis_inner_product(order::Int64,m::Int64,n::Int64)
    ncomps=3;
    rm=m%ncomps;
    rn=n%ncomps;
    ordm=div(m,ncomps);
    ordn=div(n,ncomps);

    (id1,id2,id3)=ind2sub((order,order,order),ordm);
    idsm=[id1,id2,id3];
    (id1,id2,id3)=ind2sub((order,order,order),ordn);
    idsn=[id1,id2,id3];

    val=Float64((2*pi)*(2*pi)*(2*pi));

    rm=rm+1;
    rn=rn+1;
    if idsm==idsn
        if rm==rn
            if rm==1
                val=val*(idsn[3]*idsn[3]+idsn[2]*idsn[2])
            elseif rm==2
                val=val*(idsn[3]*idsn[3]+idsn[1]*idsn[1]);
            else
                val=val*(idsn[2]*idsn[2]+idsn[1]*idsn[1]);
            end
        else
            val=-val*idsn[rn]*idsm[rm];
        end
    else
        val=0.0;
    end

    return val;
end

 function curl_basis_mass_matrix(order::Int64,quadlen::Int64)
    n=ordertolen(order);
    M=Array(Complex128,(n,n));
    (x,y,z,W)=trapezoid_tensor_product(-pi,pi,quadlen);
    V1=copy(W)*im;
    V2=copy(V1);
    V3=copy(V1);
    U1=copy(V1);
    U2=copy(V1);
    U3=copy(V1);

    for i = 1 : n
        for j = 1 : n
            #Evaluate basis at quadrature nodes.
            for u1 = 1 : length(x)
                for u2 = 1 : length(y)
                    for u3 = 1 : length(z)
                        (vx,vy,vz)=evaluate_curl_basis(order,i,x[u1],y[u2],z[u3]);
                        (ux,uy,uz)=evaluate_curl_basis(order,j,x[u1],y[u2],z[u3]);
                        V1[u1,u2,u3]=vx;
                        V2[u1,u2,u3]=vy;
                        V3[u1,u2,u3]=vz;

                        U1[u1,u2,u3]=ux;
                        U2[u1,u2,u3]=uy;
                        U3[u1,u2,u3]=uz;

                    end
                end
            end

            S=conj(V1).*U1+conj(V2).*U2+conj(V3).*U3;

            M[i,j]=sum(S[:].*W[:]);
        end
    end

    return M;
end

 function filter_dependent_modes(order::Int64,quadlen::Int64)
    M=curl_basis_mass_matrix(order,quadlen);
    r=rank(M);
    (Q,R,p)=qr(M,Val{true});
    return p[1:r];
end

 function filter_dependent_modes(order::Int64)
    M=full(curl_basis_mass_matrix(order));
    r=rank(M);
    (Q,R,p)=qr(M,Val{true});
    return p[1:r];
end

 function filter_dependent_modes_new(order::Int64)
    M=curl_basis_mass_matrix(order);
    (m,n)=size(M);
    S=IntSet([1:m;]);
    for s in S
      for t in S
	if t!=s
	  if M.rowval[M.colptr[t]:M.colptr[t+1]-1]==M.rowval[M.colptr[s]:M.colptr[s+1]-1]
	    v=(norm(M[:,s]'*M[:,t]))/(norm(M[:,s])*norm(M[:,t]));
	    if abs(v-1.0)<1e-3
	      print("eliminated pair=($s,$t)\n");
	      print("diff=($(s-t))\n");
	      print("rel diff=($((s-t)/n))\n");
	      delete!(S,t);
	    end
	  end
	end
      end
    end
    return [Int64(p) for p in S];
end

 function filter_dependent_modes_new(order::Int64,alpha)
    M=curl_basis_mass_matrix(order);
    (m,n)=size(M);
    S=IntSet([1:m;]);
    for s in S
      for t in S
	if t!=s
	  v=(norm(M[:,s]'*M[:,t]))/(norm(M[:,s])*norm(M[:,t]));
	  if abs(v-1.0)<alpha
	    print("eliminated pair=($s,$t)\n");
	    print("diff=($(s-t))\n");
	    print("rel diff=($((s-t)/n))\n");
	    delete!(S,t);
	  end
	end
      end
    end
    return [Int64(p) for p in S];
end

 function curl_basis_mass_matrix(p::Array{Int64},order::Int64,quadlen::Int64)
    n=length(p);
    M=Array(Complex128,(n,n));
    (x,y,z,W)=trapezoid_tensor_product(-pi,pi,quadlen);
    V1=copy(W)*im;
    V2=copy(V1);
    V3=copy(V1);
    U1=copy(V1);
    U2=copy(V1);
    U3=copy(V1);

    for i = 1 : n
        for j = 1 : n
            #Evaluate basis at quadrature nodes.
            for u1 = 1 : length(x)
                for u2 = 1 : length(y)
                    for u3 = 1 : length(z)
                        (vx,vy,vz)=evaluate_curl_basis(order,p[i],x[u1],y[u2],z[u3]);
                        (ux,uy,uz)=evaluate_curl_basis(order,p[j],x[u1],y[u2],z[u3]);
                        V1[u1,u2,u3]=vx;
                        V2[u1,u2,u3]=vy;
                        V3[u1,u2,u3]=vz;

                        U1[u1,u2,u3]=ux;
                        U2[u1,u2,u3]=uy;
                        U3[u1,u2,u3]=uz;

                    end
                end
            end

            S=conj(V1).*U1+conj(V2).*U2+conj(V3).*U3;

            M[i,j]=sum(S[:].*W[:]);
        end
    end

    return M;
end

 function curl_basis_mass_matrix(p::Array{Int64},order::Int64)
    n=length(p);
    M=Array(Complex128,(n,n));
    for i = 1 : n
        for j = 1 : n
           M[i,j]=curl_basis_inner_product(order,p[i],p[j]);
        end
    end

    return M;
end

 function curl_basis_mass_matrix(order::Int64)
    n=ordertolen(order);
    M=spzeros(Complex128,n,n);
    for i = 1 : n
        for j = 1 : n
           M[i,j]=curl_basis_inner_product(order,i,j);
        end
    end
    return M;
end

 function eval_curl_soln(
    u::Array{Complex128},
    p::Array{Int64},
    order::Int64,
    x,y,z)

    outx::Complex128=0.0;
    outy::Complex128=0.0;
    outz::Complex128=0.0;
    for i = 1 : length(u)
        (ux,uy,uz)=evaluate_curl_basis(order,p[i],x,y,z);
        outx=outx+u[i]*ux;
        outy=outy+u[i]*uy;
        outz=outz+u[i]*uz;
    end
    return (outx,outy,outz);
end

function curl_basis_vandermondes(
    order::Int64,
    p::Array{Int64},
    qx::Array{Float64},
    qy::Array{Float64},
    qz::Array{Float64},
    qw::Array{Float64})

    n=length(p);
    nq=length(qw);
    V=Array(Array{Complex128},(3,));
    V[1]=Array(Complex128,(nq,n));
    V[2]=Array(Complex128,(nq,n));
    V[3]=Array(Complex128,(nq,n));

    for i = 1 : n
        for iq = 1 : nq
            (i1,i2,i3)=ind2sub(size(qw),iq);
            (ux,uy,uz)=evaluate_curl_basis(order,p[i],qx[i1],qy[i2],qz[i3]);
            V[1][iq,i]=ux;
            V[2][iq,i]=uy;
            V[3][iq,i]=uz;
        end
    end
    return V;
end

function curl_basis_gradient_vandermondes(
    order::Int64,
    p::Array{Int64},
    qx::Array{Float64},
    qy::Array{Float64},
    qz::Array{Float64},
    qw::Array{Float64})

    n=length(p);
    nq=length(qw);
    GV  = Array(Array{Complex128},(3,3));
    for i = 1 : 3
        for j = 1 : 3
            GV[i,j]=Array(Complex128,(nq,n));
        end
    end

    for i = 1 : n
        for iq = 1 : nq
            (i1,i2,i3)=ind2sub(size(qw),iq);
            U=evaluate_curl_basis_gradient(order,p[i],qx[i1],qy[i2],qz[i3]);
            for di=1:3
                for dj=1:3
                    GV[di,dj][iq,i]=U[di,dj];
                end
            end
        end
    end
    return GV;
end

max_stars_repo_name: ReidAtcheson/incns

max_issues_repo_name: ReidAtcheson/incns

max_forks_repo_name: ReidAtcheson/incns


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimizer.jl
# Defines optimizer structure used by EAGO. Namely, ObjectiveType, ProblemType
# EAGOParameters, InputProblem, ParsedProblem, and Optimizer.
#############################################################################

#=
Holds outputs in the Optimizer object, passing parameters EAGO specific parameters
to _options, the GlobalOptimizer is remade to ensure type stability of subroutines
prior to initializing a global solve.
=#
export Optimizer
mutable struct Optimizer{T<:AbstractFloat} <: MOI.AbstractOptimizer
    _model::InputModel{T}
    _options::GlobalOptimizerOptions{T}
    _solver::GlobalOptimizer
    _problem_type::ProblemType

    _termination_status_code::MOI.TerminationStatusCode
    _primal_status_code::MOI.ResultStatusCode

    _objective_value::T
    _objective_bound::T
    _relative_gap::T
    _iteration_count::Int
    _node_count::Int

    _is_silent::Bool

    _start_time::T
    _time_limit::T
    _time_left::T
    _run_time::T
    _parse_time::T

    _solution::Vector{T}
    _constraint_primal::Dict{CI,Union{Vector{T},T}}
    _input_to_solution_map
end
function Optimizer{T}(; kwargs...) where T<:Real
    model   = InputModel{T}()
    options = GlobalOptimizerOptions{T}(;kwargs...)
    solver  = GlobalOptimizer(T)
    problem_type = UNCLASSIFIED

    termination_status_code = MOI.OPTIMIZE_NOT_CALLED
    primal_status_code      = MOI.OTHER_RESULT_STATUS

    objective_value = Inf
    objective_bound = -Inf
    relative_gap = Inf
    iteration_count = 0
    node_count      = 0

    is_silent       = false

    start_time   = 0.0
    time_limit   = Inf
    time_left    = 1000.0
    run_time     = 0.0
    parse_time   = 0.0

    solution          = T[]
    constraint_primal = Dict{CI,T}()
    input_to_solution_map = nothing

    return Optimizer{T}(model,
                        options,
                        solver,
                        problem_type,
                        termination_status_code,
                        primal_status_code,
                        objective_value,
                        objective_bound,
                        relative_gap,
                        iteration_count,
                        node_count,
                        is_silent,
                        start_time,
                        time_limit,
                        time_left,
                        run_time,
                        parse_time,
                        solution,
                        constraint_primal,
                        input_to_solution_map)
end
Optimizer(; kwargs...) = Optimizer{Float64}(; kwargs...)

for f in fieldnames(GlobalOptimizerOptions)
    @eval @inline $(Symbol("_"*String(f)))(x::Optimizer) = x._options.$f
end

_input_model(m::Optimizer{T}) where T = m._model._input_model
_input_nlp_data(m::Optimizer) = m._model._nlp_data

function _set_cons_primal!(m::Optimizer{T}, ci::CI{F,S}, v::T) where {F <: MOI.AbstractFunction,
                                                                     S <: MOI.AbstractScalarSet,
                                                                     T <: AbstractFloat}
    m._constraint_primal[ci] = v
    return
end

function _set_cons_primal!(m::Optimizer{T}, ci::CI{F,S}, v::Vector{T}) where {F <: MOI.AbstractFunction,
                                                                              S <: MOI.AbstractVectorSet,
                                                                              T <: AbstractFloat}
    m._constraint_primal[ci] = v
    return
end


text:
#####################################################
### Project of IN480 - From Python 2.7 to Julia 0.6
### --- SIMPLEXN - parallel version ---
### Authors: Fabio Fatelli, Emanuele Loprevite
#####################################################

# Generation of the output model vertices in a multiple extrusion of a LAR model
@everywhere function plarExtrude1{T<:Real}(model::Tuple{Array{Array{T,1},1},Array{Array{Int64,1},1}}, pattern::Array{Int64,1})
	V, FV = model
	d, m = length(FV[1]), length(pattern)
	dd1 = d*(d+1)
	coords = cumsum(append!([0],abs.(pattern))) # built-in function cumsum
	outVertices = @spawn [vcat(v,z) for z in coords for v in V]
	offset, outcells, rangelimit = length(V), SharedArray{Int64}(m,dd1*length(FV)), d*m
	@sync @parallel for j in 1:length(FV)
		tube = [v+k*offset for k in 0:m for v in FV[j]]
		celltube = Int64[]
		celltube = @sync @parallel (append!) for k in 1:rangelimit
			tube[k:k+d]
		end
		outcells[:,(j-1)*dd1+1:(j-1)*dd1+dd1] = reshape(celltube,dd1,m)'
	end
	p = convert(SharedArray,find(x->x>0,pattern))
	cellGroups = SharedArray{Int64}(length(p),size(outcells)[2])
	@sync @parallel for k in 1:length(p)
		cellGroups[k,:] = outcells[p[k],:]
	end
	cellGroupsL = vec(cellGroups')
	outCellGroups = Array{Int64,1}[]
	outCellGroups = @parallel (append!) for k in 1:d+1:length(cellGroupsL)
			[cellGroupsL[k:k+d]]
		end
	return fetch(outVertices), outCellGroups
end

# Generation of simplicial grids of any dimension and shape
@everywhere function plarSimplexGrid1(shape::Array{Int64,1})
	model = [Int64[]],[[0]] # the empty simplicial model
	for item in shape # no parallel
		model = plarExtrude1(model,repmat([1],item))
	end
	return model
end

# Extraction of non-oriented (d-1)-facets of d-dimensional simplices
@everywhere using Combinatorics # for combinations() function

@everywhere function plarSimplexFacets(simplices::Array{Array{Int64,1},1})
	out = Array{Int64,1}[]
	d = length(simplices[1])
	out = @parallel (append!) for simplex in simplices
			collect(combinations(simplex,d-1))
		end
	return sort!(unique(out),lt=lexless)
end

# Transformation to triangles by sorting circularly the vertices of faces
@everywhere function pquads2tria{T<:Real}(model::Tuple{Array{Array{T,1},1},Array{Array{Int64,1},1}})
	V, FV = model
	if typeof(V) != Array{Array{Float64,1},1}
		V = convert(Array{Array{Float64,1},1},V)
	end
	out = Array{Int64,1}[]
	nverts = length(V)-1
	for face in FV # no parallel
		arr = [V[v+1] for v in face]
		centroid = sum(arr)/length(arr)
		append!(V,[centroid])
		nverts += 1
		v1, v2 = V[face[1]+1]-centroid, V[face[2]+1]-centroid
		v3 = cross(v1,v2)
		if norm(v3) < 1/(10^3)
			v1, v2 = V[face[1]+1]-centroid, V[face[3]+1]-centroid
			v3 = cross(v1,v2)
		end
		transf = inv(hcat(v1,v2,v3)')
		verts = [(V[v+1]'*transf)'[1:end-1] for v in face]
		tcentroid = sum(verts)/length(verts)
		tverts = pmap(x->x-tcentroid,verts)
		iterator = collect(zip(tverts,face))
		rverts = [[atan2(reverse(iterator[i][1])...),iterator[i][2]] for i in 1:length(iterator)]
		rvertsS = sort(rverts,lt=(x,y)->isless(x[1],y[1]))
		ord = [pair[2] for pair in rvertsS]
		append!(ord,ord[1])
		edges = [[i[2],ord[i[1]+1]] for i in enumerate(ord[1:end-1])]
		triangles = pmap(x->prepend!(x,nverts),edges)
		append!(out,triangles)
	end
	return V, out
end

max_stars_repo_name: EmaLoprevite/Simplexn

max_issues_repo_name: EmaLoprevite/Simplexn

max_forks_repo_name: EmaLoprevite/Simplexn


text:
using Libdl
using PackageCompiler

# this path must match the one used in `runpluto.sh`, i.e. julia should be
# started with this system image
# sysimage_path = joinpath("usr", "local", string("sys.", Libdl.dlext))
sysimage_path = "/usr/local/sys.so"

create_sysimage([:Pluto, :PlutoUI, :Distributions, :Plots];
                sysimage_path = sysimage_path,
                precompile_execution_file = "warmup.jl",
                cpu_target = PackageCompiler.default_app_cpu_target())


text:
# this file is used to resize the road-network image .svg files so that they all
# have a square aspect ratio. This is important for displaying the images in the
# experiment
include("utilities.jl")

fdir = "./imgs"
files = searchdir(fdir,".svg")
mod_imgs_fldr = "./imgs/squared"
function parse_strings(s_ary)
    ary = Array{Float64}(length(s_ary))
    for i = 1:length(s_ary)
        ary[i] = parse(Float64,s_ary[i])
    end
    return ary
end

for f in files
    txt = open(joinpath(fdir,f)) do file
        readlines(file)
    end
    loi = txt[2] # loi -> line of interest

    # need to find maximum of width and height, then cange the smaller to match it
    # then offset the viewBox to keep the image in the center

    #capture number with 'pt' at the end, this will be width and height
    r1 = r"(-?\d+(?>\.\d+)?)pt.+?(-?\d+(?>\.\d+)?)pt"
    # capture 4 number sequence, this is the 'viewBox' stuff
    r2 = r"(-?\d+(?>\.\d+)?) (-?\d+(?>\.\d+)?) (-?\d+(?>\.\d+)?) (-?\d+(?>\.\d+)?)"

    w_h = match(r1,loi)
    wh = parse_strings(w_h.captures)
    box = match(r2,loi)
    vb = parse_strings(box.captures)

    ind_max_d = indmax(wh)
    ind_min_d = indmin(wh)
    dmax = wh[ind_max_d]
    dmin = wh[ind_min_d]
    if dmax == dmin
        # dims already square
        continue
    end
    wh_diff = dmax - dmin
    wh[ind_min_d] = wh[ind_max_d]

    if ind_min_d == 1
        #we need to  modify x
        vb[1] -= wh_diff/2
        vb[3] += wh_diff/2
    else
        # we need to modify y
        vb[2] -= wh_diff/2
        vb[4] += wh_diff/2
    end

    mod_loi = loi
    mod_loi = replace(mod_loi,w_h.match,@sprintf("%0.3fpt\" height=\"%0.3fpt",wh[1],wh[2]))

    mod_loi = replace(mod_loi,box.match,@sprintf("%0.3f %0.3f %0.3f %0.3f",vb[1],vb[2],vb[3],vb[4]))

    println(mod_loi)

    txt[2] = mod_loi

    open("$mod_imgs_fldr/$f","w") do f
        writedlm(f,txt,'\n',quotes=false)
    end

end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text: Dict{String, Any}("34-11" => Dict{String, Any}("value" => "23", "type" => "integer"), "000111" => Dict{String, Any}("value" => "leading", "type" => "string"), "a-a-a" => Dict{String, Any}("_" => Dict{String, Any}("value" => "false", "type" => "bool")), "alpha" => Dict{String, Any}("value" => "a", "type" => "string"), "under_score" => Dict{String, Any}("value" => "___", "type" => "string"), "with-dash" => Dict{String, Any}("value" => "dashed", "type" => "string"), "123" => Dict{String, Any}("value" => "num", "type" => "string"), "2018_10" => Dict{String, Any}("001" => Dict{String, Any}("value" => "1", "type" => "integer")), "one1two2" => Dict{String, Any}("value" => "mixed", "type" => "string"), "10e3" => Dict{String, Any}("value" => "false float", "type" => "string"))


text:
"""
    Atom
A wrapped DexLang value.

Scalar values can be converted to julia objects using [`juliaize`](@ref), or `convert`.
Functions can be called directly, if you pass in atoms,
or they can be compiled and made usable on julia objects by using [`NativeFunction`](@ref).
"""
struct Atom
    ptr::Ptr{HsAtom}
    ctx::Ptr{HsContext}
end

Base.show(io::IO, atom::Atom) = show(io, print(atom.ptr))

juliaize(x::Atom) = juliaize(x.ptr)
Base.convert(::Type{T}, atom::Atom) where {T<:Number} = convert(T, juliaize(atom))

function (self::Atom)(args...)
    # TODO: Make those calls more hygenic
    env = self.ctx
    pieces = (self, args...)
    for (i, atom) in enumerate(pieces)
        # NB: Atoms can contain arbitrary references
        if atom.ctx !== PRELUDE && atom.ctx !== self.ctx
            throw(ArgumentError("Mixing atoms coming from different Dex modules is not supported yet!"))
        end
        old_env, env = env, insert(env, "julia_arg$i", atom.ptr)
        destroy_context(old_env)
    end
    return evaluate(join("julia_arg" .* string.(eachindex(pieces)), " "), self.ctx, env)
end

mutable struct DexModule
    # Needs to be mutable struct so can attach finalizer
    ctx::Ptr{HsContext}
end

"""
    DexModule(str)

For running 1 or more Dex expressions, and keeping the state.
You can get them back out of the module using `getproperty`.
They are returned as [`DexCall.Atom`](@ref)s.

# Example:

```julia
julia> m = DexModule(raw"""
       x = 42
       y = 2 * x
       """)
DexModule(Ptr{DexCall.HsContext} @0x0000000000000031)

julia> m.x
"42"

julia> m.y
"84"
```
"""
function DexModule(source::AbstractString)
    ctx = dex_eval(PRELUDE, source)
    ctx == C_NULL && throw_from_dex()
    m =  DexModule(ctx)
    finalizer(m) do _m
        destroy_context(getfield(_m, :ctx))
    end
    return m
end

function Base.getproperty(m::DexModule, name::Symbol)
    ctx = getfield(m, :ctx)
    ret = lookup(ctx, string(name))
    ret == C_NULL && throw_from_dex()
    return Atom(ret, ctx)
end

@doc raw"""
    evaluate(str)

A friendly function for running Dex code.
The string `str` must contain a single Dex expression.
Return a [`DexCall.Atom`](@ref)
# Example:
```julia
julia> evaluate(raw"sum $ for i. exp [log 2.0, log 4.0].i")
"6."
```
"""
function evaluate(str, _module=PRELUDE, env=_module)
    result = eval_expr(env, str)
    result == C_NULL && throw_from_dex()
    return Atom(result, _module)
end

max_stars_repo_name: AlexDikelsky/dex-lang

max_issues_repo_name: AlexDikelsky/dex-lang

max_forks_repo_name: AlexDikelsky/dex-lang


text:
"""
$(SIGNATURES)

GetNumberOfAperiodicities provides the number of dimensions for aperiodicity
coding. It is determined by only fs.

**Parameters**

- `fs` : Sampleing frequency

**Returns**

- `n` : Number of aperiodicities
"""
function get_number_of_aperiodicities(fs)
    ccall((:GetNumberOfAperiodicities, libworld), Cint, (Cint,), fs)
end

"""
$(SIGNATURES)

CodeAperiodicity codes the aperiodicity. The number of dimensions is
determined by fs.

**Parameters**

- `aperiodicity` : Aperiodicity before coding
- `fs` : Sampling frequency
- `fftsize` : FFT size (default : `get_fftsize_for_cheaptrick(fs)`)

**Returns**

- `coded_aperiodicity` : Coded aperiodicity
"""
function code_aperiodicity(aperiodicity, fs,
                           fftsize=get_fftsize_for_cheaptrick(fs))
    n = get_number_of_aperiodicities(fs)
    f0_length = size(aperiodicity, 2)
    coded_aperiodicity = zeros(Cdouble, n, f0_length)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    ccoded_aperiodicity = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(ccoded_aperiodicity, coded_aperiodicity)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    caperiodicity = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(caperiodicity, aperiodicity)

    ccall((:CodeAperiodicity, libworld), Cvoid,
          (Ptr{Ptr{Cdouble}}, Cint, Cint, Cint, Ptr{Ptr{Cdouble}}),
          caperiodicity, f0_length, fs, fftsize, ccoded_aperiodicity)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:f0_length, j=1:n
        coded_aperiodicity[j,i] = unsafe_load(ccoded_aperiodicity[i], j)
    end

    coded_aperiodicity
end

"""
$(SIGNATURES)

DecodeAperiodicity decoes the coded aperiodicity.

**Parameters**

- `coded_aperiodicity` : Coded aperiodicity
- `fs` : Sampling frequency
- `fftsize` : FFT size (default : `get_fftsize_for_cheaptrick(fs)`)

**Returns**

- `aperiodicity` : Decoded aperiodicity
"""
function decode_aperiodicity(coded_aperiodicity, fs,
                             fftsize=get_fftsize_for_cheaptrick(fs, CheapTrickOption(fs)))
    freqbins = fftsize>>1 + 1
    f0_length = size(coded_aperiodicity, 2)
    aperiodicity = zeros(Cdouble, freqbins, f0_length)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    ccoded_aperiodicity = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(ccoded_aperiodicity, coded_aperiodicity)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    caperiodicity = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(caperiodicity, aperiodicity)

    ccall((:DecodeAperiodicity, libworld), Cvoid,
          (Ptr{Ptr{Cdouble}}, Cint, Cint, Cint, Ptr{Ptr{Cdouble}}),
          ccoded_aperiodicity, f0_length, fs, fftsize, caperiodicity)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:f0_length, j=1:freqbins
        aperiodicity[j,i] = unsafe_load(caperiodicity[i], j)
    end

    aperiodicity
end

"""
$(SIGNATURES)

CodeSpectralEnvelope codes the spectral envelope.

**Parameters**

- `spectrogram` : spectrogram (time sequence of spectral envelope)
- `fs` : Sampling frequency
- `fftsize` : FFT size
- `number_of_dimentions` : Number of dimentions for coded spectral envelope

**Returns**

- `coded_spectral_envelope` : Coded spectral envelope
"""
function code_spectral_envelope(spectrogram, fs,
                                fftsize, number_of_dimentions)
    nd = number_of_dimentions
    f0_length = size(spectrogram, 2)
    coded_spectral_envelope = zeros(Cdouble, nd, f0_length)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    ccoded_spectral_envelope = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(ccoded_spectral_envelope, coded_spectral_envelope)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(cspectrogram, spectrogram)

    ccall((:CodeSpectralEnvelope, libworld), Cvoid,
          (Ptr{Ptr{Cdouble}}, Cint, Cint, Cint, Cint, Ptr{Ptr{Cdouble}}),
          cspectrogram, f0_length, fs, fftsize, nd, ccoded_spectral_envelope)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:f0_length, j=1:nd
        coded_spectral_envelope[j,i] = unsafe_load(ccoded_spectral_envelope[i], j)
    end

    coded_spectral_envelope
end

"""
$(SIGNATURES)

DecodeSpectralEnvelope decodes the spectral envelope.

**Parameters**

- `coded_spectral_envelope` : Coded spectral envelope
- `fs` : Sampling frequency
- `fftsize` : FFT size

**Returns**

- `spectrogram` : decoded spectral envelope
"""
function decode_spectral_envelope(coded_spectral_envelope, fs, fftsize)
    freqbins = fftsize>>1 + 1
    nd = size(coded_spectral_envelope, 1)
    f0_length = size(coded_spectral_envelope, 2)
    spectrogram = zeros(Cdouble, freqbins, f0_length)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    ccoded_spectral_envelope = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(ccoded_spectral_envelope, coded_spectral_envelope)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array{Ptr{Cdouble}}(undef, f0_length)
    ptrarray2d!(cspectrogram, spectrogram)

    ccall((:DecodeSpectralEnvelope, libworld), Cvoid,
          (Ptr{Ptr{Cdouble}}, Cint, Cint, Cint, Cint, Ptr{Ptr{Cdouble}}),
          ccoded_spectral_envelope, f0_length, fs, fftsize, nd, cspectrogram)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:f0_length, j=1:freqbins
        spectrogram[j,i] = unsafe_load(cspectrogram[i], j)
    end

    spectrogram
end


text:
using FaultTolerantControl
const FTC = FaultTolerantControl

function test()
    allocator = AdaptiveAllocator()
    x0 = State(allocator)()
    tf = 10.0
    prob, df = sim(
                   x0,
                   apply_inputs(Dynamics!(allocator); Θ̂̇=zeros(6, 4));
                   tf=tf
                  )
    df
end


text:
# Example adapted from Daniel Shiffman - The Coding Train
# https://youtu.be/0jjeOYMjmDU?list=PLRqwX-V7Uu6ZiZxtDDRCi6uhfTH4FilpH

using Alexya

@init "Fractal Tree" 800 600 # create a Canvas
@layout aside(:v, 200)

slider = @create Slider(0:π/12:2π; init=π/4, @margin(20))
φ = π/4

@use function update()
    global φ = value(slider)

    background("#515151")
    origin((@width)/2, @height)
    sethue("white")
    branch(150)
end

function branch(len)
    line(O, Point(0, -len), :stroke)
    translate(0, -len)

    if len > 4
        gsave()
        rotate(φ)
        branch(.67len)
        grestore()

        gsave()
        rotate(-φ)
        branch(.67len)
        grestore()
    end
end

start()


text:
include("utilities.jl")
global const S_box_AES = [
0x63, 0x7c, 0x77, 0x7b, 0xf2, 0x6b, 0x6f, 0xc5, 0x30, 0x01, 0x67, 0x2b, 0xfe, 0xd7, 0xab, 0x76,
0xca, 0x82, 0xc9, 0x7d, 0xfa, 0x59, 0x47, 0xf0, 0xad, 0xd4, 0xa2, 0xaf, 0x9c, 0xa4, 0x72, 0xc0,
0xb7, 0xfd, 0x93, 0x26, 0x36, 0x3f, 0xf7, 0xcc, 0x34, 0xa5, 0xe5, 0xf1, 0x71, 0xd8, 0x31, 0x15,
0x04, 0xc7, 0x23, 0xc3, 0x18, 0x96, 0x05, 0x9a, 0x07, 0x12, 0x80, 0xe2, 0xeb, 0x27, 0xb2, 0x75,
0x09, 0x83, 0x2c, 0x1a, 0x1b, 0x6e, 0x5a, 0xa0, 0x52, 0x3b, 0xd6, 0xb3, 0x29, 0xe3, 0x2f, 0x84,
0x53, 0xd1, 0x00, 0xed, 0x20, 0xfc, 0xb1, 0x5b, 0x6a, 0xcb, 0xbe, 0x39, 0x4a, 0x4c, 0x58, 0xcf,
0xd0, 0xef, 0xaa, 0xfb, 0x43, 0x4d, 0x33, 0x85, 0x45, 0xf9, 0x02, 0x7f, 0x50, 0x3c, 0x9f, 0xa8,
0x51, 0xa3, 0x40, 0x8f, 0x92, 0x9d, 0x38, 0xf5, 0xbc, 0xb6, 0xda, 0x21, 0x10, 0xff, 0xf3, 0xd2,
0xcd, 0x0c, 0x13, 0xec, 0x5f, 0x97, 0x44, 0x17, 0xc4, 0xa7, 0x7e, 0x3d, 0x64, 0x5d, 0x19, 0x73,
0x60, 0x81, 0x4f, 0xdc, 0x22, 0x2a, 0x90, 0x88, 0x46, 0xee, 0xb8, 0x14, 0xde, 0x5e, 0x0b, 0xdb,
0xe0, 0x32, 0x3a, 0x0a, 0x49, 0x06, 0x24, 0x5c, 0xc2, 0xd3, 0xac, 0x62, 0x91, 0x95, 0xe4, 0x79,
0xe7, 0xc8, 0x37, 0x6d, 0x8d, 0xd5, 0x4e, 0xa9, 0x6c, 0x56, 0xf4, 0xea, 0x65, 0x7a, 0xae, 0x08,
0xba, 0x78, 0x25, 0x2e, 0x1c, 0xa6, 0xb4, 0xc6, 0xe8, 0xdd, 0x74, 0x1f, 0x4b, 0xbd, 0x8b, 0x8a,
0x70, 0x3e, 0xb5, 0x66, 0x48, 0x03, 0xf6, 0x0e, 0x61, 0x35, 0x57, 0xb9, 0x86, 0xc1, 0x1d, 0x9e,
0xe1, 0xf8, 0x98, 0x11, 0x69, 0xd9, 0x8e, 0x94, 0x9b, 0x1e, 0x87, 0xe9, 0xce, 0x55, 0x28, 0xdf,
0x8c, 0xa1, 0x89, 0x0d, 0xbf, 0xe6, 0x42, 0x68, 0x41, 0x99, 0x2d, 0x0f, 0xb0, 0x54, 0xbb, 0x16
];

global const invS_box_AES = [
0x52, 0x09, 0x6a, 0xd5, 0x30, 0x36, 0xa5, 0x38, 0xbf, 0x40, 0xa3, 0x9e, 0x81, 0xf3, 0xd7, 0xfb,
0x7c, 0xe3, 0x39, 0x82, 0x9b, 0x2f, 0xff, 0x87, 0x34, 0x8e, 0x43, 0x44, 0xc4, 0xde, 0xe9, 0xcb,
0x54, 0x7b, 0x94, 0x32, 0xa6, 0xc2, 0x23, 0x3d, 0xee, 0x4c, 0x95, 0x0b, 0x42, 0xfa, 0xc3, 0x4e,
0x08, 0x2e, 0xa1, 0x66, 0x28, 0xd9, 0x24, 0xb2, 0x76, 0x5b, 0xa2, 0x49, 0x6d, 0x8b, 0xd1, 0x25,
0x72, 0xf8, 0xf6, 0x64, 0x86, 0x68, 0x98, 0x16, 0xd4, 0xa4, 0x5c, 0xcc, 0x5d, 0x65, 0xb6, 0x92,
0x6c, 0x70, 0x48, 0x50, 0xfd, 0xed, 0xb9, 0xda, 0x5e, 0x15, 0x46, 0x57, 0xa7, 0x8d, 0x9d, 0x84,
0x90, 0xd8, 0xab, 0x00, 0x8c, 0xbc, 0xd3, 0x0a, 0xf7, 0xe4, 0x58, 0x05, 0xb8, 0xb3, 0x45, 0x06,
0xd0, 0x2c, 0x1e, 0x8f, 0xca, 0x3f, 0x0f, 0x02, 0xc1, 0xaf, 0xbd, 0x03, 0x01, 0x13, 0x8a, 0x6b,
0x3a, 0x91, 0x11, 0x41, 0x4f, 0x67, 0xdc, 0xea, 0x97, 0xf2, 0xcf, 0xce, 0xf0, 0xb4, 0xe6, 0x73,
0x96, 0xac, 0x74, 0x22, 0xe7, 0xad, 0x35, 0x85, 0xe2, 0xf9, 0x37, 0xe8, 0x1c, 0x75, 0xdf, 0x6e,
0x47, 0xf1, 0x1a, 0x71, 0x1d, 0x29, 0xc5, 0x89, 0x6f, 0xb7, 0x62, 0x0e, 0xaa, 0x18, 0xbe, 0x1b,
0xfc, 0x56, 0x3e, 0x4b, 0xc6, 0xd2, 0x79, 0x20, 0x9a, 0xdb, 0xc0, 0xfe, 0x78, 0xcd, 0x5a, 0xf4,
0x1f, 0xdd, 0xa8, 0x33, 0x88, 0x07, 0xc7, 0x31, 0xb1, 0x12, 0x10, 0x59, 0x27, 0x80, 0xec, 0x5f,
0x60, 0x51, 0x7f, 0xa9, 0x19, 0xb5, 0x4a, 0x0d, 0x2d, 0xe5, 0x7a, 0x9f, 0x93, 0xc9, 0x9c, 0xef,
0xa0, 0xe0, 0x3b, 0x4d, 0xae, 0x2a, 0xf5, 0xb0, 0xc8, 0xeb, 0xbb, 0x3c, 0x83, 0x53, 0x99, 0x61,
0x17, 0x2b, 0x04, 0x7e, 0xba, 0x77, 0xd6, 0x26, 0xe1, 0x69, 0x14, 0x63, 0x55, 0x21, 0x0c, 0x7d
];

global const S_box_PRESENT = [0xC, 0x5, 0x6, 0xB, 0x9, 0x0, 0xA, 0xD, 0x3, 0xE, 0xF, 0x8, 0x4, 0x7, 0x1, 0x2];
global const invS_box_PRESENT = [
0, 16, 32, 48, 1, 17, 33, 49, 2, 18, 34, 50, 3, 19, 35, 51, 4, 20, 36, 52, 5, 21, 37, 53, 6, 22, 38, 54, 7, 23, 39, 55, 8, 24,
40, 56, 9, 25, 41, 57, 10, 26, 42, 58, 11, 27, 43, 59, 12, 28, 44, 60, 13, 29, 45, 61, 14, 30, 46, 62, 15, 31, 47, 63];
#add more S-P components here if needed
global HW = Dict()

global S_boxes = Dict()
S_boxes["AES"] = S_box_AES
S_boxes["PRESENT"] = S_box_PRESENT
#add more as more components are added

global invS_boxes = Dict()
invS_boxes["AES"] = invS_box_AES
invS_boxes["PRESENT"] = invS_box_PRESENT
#add more as more components are added

max_stars_repo_name: udion/zypher

max_issues_repo_name: udion/zypher

max_forks_repo_name: udion/zypher


text:
using ParallelRandomFields
using ParallelRandomFields.grf2D_Threads

using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@init_parallel_stencil(Threads, Float64, 2)

using MAT, Plots

@views function generate_grf2D()
    cov_typ = "expon"
    do_viz  = true
    do_save = false
    # Physics
    lx, ly  = 100.0, 100.0  # domain size
    sf      = 1.0           # standard deviation
    # -- exponential setup
    cl_e    = (10.0, 8.0)   # correlation lengths in [x,y]
    # -- gaussian setup
    cl_g    = 5.0           # correlation length isotropic
    k_m     = 100.0         # maximum value of the wave number
    # Numerics
    nh      = 10000         # inner parameter, number of harmonics
    nx, ny  = 64, 64        # numerical grid resolution
    # Derived numerics
    dx, dy  = lx/nx, ly/ny  # numerical grid step size
    # Array allocation
    Yf      = @zeros(nx, ny)
    # Visu init
    if do_viz
        X, Y = -lx/2:dx:lx/2, -ly/2:dy:ly/2
    end
    if cov_typ=="expon"
        # Generate the 2D exponential covariance function
        grf2D_expon!(Yf, sf, cl_e, nh, nx, ny, dx, dy; do_reset=true)
    elseif cov_typ=="gauss"
        # Generate the 2D Gaussian covariance function
        grf2D_gauss!(Yf, sf, cl_g, nh, k_m, nx, ny, dx, dy; do_reset=true)
    else
        error("trying to run with undefined covariance function")
    end
    # Visualisation
    if do_viz
        display(heatmap(X, Y, Array(Yf)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:hot, title="2D RandomField"))
    end
    if do_save  file = matopen("grf2D_$(cov_typ).mat", "w"); write(file, "grf2D", Array(Yf)); close(file)  end
    return
end

generate_grf2D()


text:
function ridge(y, x, k)
    # standardize regressors (except the constant)
    s = std(x,dims=1)
    A = eye(size(x,2))
    for i = 1:size(x,2)
        if s[i]==0
            s[i] = 1.0
            A[i,i] = 0.0 # don't shrink constant
        end
    end
    x = x ./ s
    βhat = (inv(x'x + k*A)*x'y) ./ s'
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
"""
    Triangle(a, b, c)

Create a `Triangle` that can be used for spatial lookup.

`a`, `b`, and `c` are the points of the triangle on the x/y plane.

See also `AbstractRegion`
"""
struct Triangle{T <: Real} <: AbstractRegion{T}
    a::SVector{2,T}
    b::SVector{2,T}
    c::SVector{2,T}
end

@inline function Triangle(a::AbstractVector{<:Real}, b::AbstractVector{<:Real}, c::AbstractVector{<:Real})
    T = promote_type(eltype(a), eltype(b), eltype(c))
    Triangle(convert(SVector{2,T}, a), convert(SVector{2,T}, b), convert(SVector{2,T}, c))
end

"""
    area(triangle::Triangle)

Returns the area of the `triangle`.
"""
@inline function area(t::Triangle{T}) where {T}
    # Should possibly use same convention as LineString and Polygon, and return negative area for counter-clockwise triangles
    convert(T, 0.5) * abs((t.b - t.a) × (t.c - t.a))
end

function boundingbox(t::Triangle{T}) where T
    points = [t.a, t.b, t.c]

    xmin, xmax = extrema(map(x -> x[1], points))
    ymin, ymax = extrema(map(x -> x[2], points))
    zmin, zmax = typemin(T), typemax(T)
    return BoundingBox{T}(xmin, ymin, zmin, xmax, ymax, zmax)
end

function in(p::StaticVector{3,T}, t::Triangle) where T <: Real
    a, b, c = t.a, t.b, t.c
    p2d = p[1:2]

    areaABC = cross(b - a, c - a)

    areaPBC = cross(b - p2d, c - p2d)
    areaPCA = cross(c - p2d, a - p2d)
    areaPAB = cross(a - p2d, b - p2d)

    u = areaPBC / areaABC
    v = areaPCA / areaABC
    w = areaPAB / areaABC

    return (u + v + w) ≈ 1.0 && u > 0 && v > 0 && w > 0 && u < 1 && v < 1 && w < 1
end


text:
using MVDA
using LinearAlgebra, Random, StableRNGs, Parameters, Statistics
using Test

@testset "Projections" begin
    @testset "L0" begin
        for _ in 1:10
            number_components = 1000
            x = 10 * randn(number_components)
            k = 50
            P = MVDA.L0Projection(number_components)
            xproj = P(copy(x), k)

            nzidx = findall(xi -> xi != 0, xproj)
            @test length(nzidx) ≤ k
            @test sort!(xproj[nzidx], by=abs, rev=true) == sort(x, by=abs, rev=true)[1:k]
        end
    end
end

function test_on_dataset(prob, target, X, k)
    n, p, c = length(target), size(X, 2), length(unique(target))
    has_intercept = prob.intercept

    @testset "MVDAProblem" begin
        # coefficient shape
        for field in (:coeff, :coeff_prev, :proj, :grad)
            arr = getfield(prob, field)
            @test size(arr.all, 1) == p + has_intercept && size(arr.all, 2) == c-1
            @test all(j -> length(arr.dim[j]) == p + has_intercept, eachindex(arr.dim))
        end

        # residual shape
        r1 = prob.res.main.all
        r2 = prob.res.dist.all
        r3 = prob.res.weighted.all
        @test size(r1, 1) == n && size(r1, 2) == c-1
        @test size(r2, 1) == p + has_intercept && size(r2, 2) == c-1
        @test size(r3, 1) == n && size(r3, 2) == c-1

        # verify problem dimensions
        dims = MVDA.probdims(prob)
        @test n == dims[1] && p == dims[2] && c == dims[3]

        # check vertex definitions
        v = prob.vertex
        d = [norm(v[j] - v[i]) for i in 1:c for j in i+1:c]
        dmin, dmax = extrema(d)
        @test all(j -> norm(v[j]) ≈ 1, eachindex(v))    # vertices lie on unit circle
        @test dmin ≈ dmax                               # vertices are equidistant
    end

    @testset "Evaluation" begin
        @unpack Y, X, coeff, res, proj, grad = prob
        B = coeff.all
        P = proj.all
        T = Float64

        # Initialize coefficients and parameters
        rng = StableRNG(1903)
        randn!(rng, B)
        ϵ = 0.5 * sqrt(2*c/(c-1))
        ρ = 4.7

        # Initialize other data structures.
        extras = MVDA.__mm_init__(MMSVD(), prob, nothing)
        MVDA.__mm_update_sparsity__(MMSVD(), prob, ϵ, ρ, k, extras)
        MVDA.__mm_update_rho__(MMSVD(), prob, ϵ, ρ, k, extras)
        operator = extras.projection

        # Set projection.
        MVDA.apply_projection(operator, prob, k)

        # Compute residuals and gradient.
        XB = X*B
        R1 = Y - XB
        R2 = P - B

        normr = [norm(ri) for ri in eachrow(R1)]
        W = Diagonal( [ifelse(normr[i] ≤ ϵ, 0.0, (normr[i]-ϵ) / normr[i]) for i in eachindex(normr)] )
        R1_scaled = 1/sqrt(n) * W * R1
        Z = similar(Y)
        @views for i in axes(Z, 1)
            Z[i, :] .= ifelse(normr[i] ≤ ϵ, XB[i, :], W[i,i] * Y[i, :] + (1-W[i,i]) * XB[i, :])
        end
        G = -1/n * X' * (Z - XB) - ρ * (P - B)

        @testset "Residuals" begin
            MVDA.__evaluate_residuals__(prob, ϵ, extras, true, true, true)
            @test all(j -> res.main.dim[j] ≈ 1/sqrt(n) * view(R1, :, j), 1:c-1)
            @test all(j -> res.dist.dim[j] ≈ view(R2, :, j), 1:c-1)
            @test all(j -> res.weighted.dim[j] ≈ view(R1_scaled, :, j), 1:c-1)
            @test all(i -> view(extras.Z, i, :) ≈ view(Z, i, :), 1:n)
        end

        @testset "Gradient" begin
            MVDA.__evaluate_gradient__(prob, ρ, extras)
            @test all(j -> grad.dim[j] ≈ view(G, :, j), 1:c-1)
        end

        @testset "Objective" begin
            loss, obj, dist, gradsq = MVDA.__evaluate_objective__(prob, ϵ, ρ, extras)
            @test loss ≈ dot(R1_scaled, R1_scaled)
            @test obj ≈ 1//2 * (dot(R1_scaled, R1_scaled) + ρ * dot(R2, R2))
            @test dist ≈ dot(R2, R2)
            @test gradsq ≈ dot(G, G)
        end
    end

    @testset "Descent Property" begin
        ϵ = 0.5 * sqrt(2*c/(c-1))
        ρ = 1.234

        @testset "$(algorithm)" for algorithm in (MMSVD(), SD(),)
            function test_descent_property(s, threshold)
                B = prob.coeff.all
                rng = StableRNG(1903)
                randn!(rng, B)
                copyto!(prob.coeff_prev.all, B)

                _, _, obj0A, _, _ = MVDA.anneal(algorithm, prob, ϵ, ρ, s, ninner=0, gtol=1e-6, nesterov_threshold=threshold)
                _, _, obj0B, _, _ = MVDA.anneal(algorithm, prob, ϵ, ρ, s, ninner=0, gtol=1e-6, nesterov_threshold=threshold)
                _, _, obj1, _, _ = MVDA.anneal(algorithm, prob, ϵ, ρ, s, ninner=1, gtol=1e-6, nesterov_threshold=threshold)
                _, _, obj2, _, _ = MVDA.anneal(algorithm, prob, ϵ, ρ, s, ninner=1, gtol=1e-6, nesterov_threshold=threshold)
                _, _, obj100, _, _ = MVDA.anneal(algorithm, prob, ϵ, ρ, s, ninner=98, gtol=1e-6, nesterov_threshold=threshold)
                _, _, objfinal, _, gradsq = MVDA.anneal(algorithm, prob, ϵ, ρ, s, ninner=10^4, gtol=1e-8, nesterov_threshold=threshold)

                @test obj0A == obj0B # no iterations
                @test obj0A > obj1   # decrease after 1 iteration
                @test obj1 > obj2    # decrease after 1 iteration
                @test obj0A > obj100 # decrease after 100 iterations
                @test obj0A > objfinal # decrease at final estimate
                @test gradsq < 1e-8 # convergence
                @test all(!isnan, B) # no instability
            end

            # check for different model sizes
            for s in (0.0, 0.25, 0.5, 0.75)
                # w/o Nesterov acceleration
                test_descent_property(s, 100)

                # w/ Nesterov acceleration
                test_descent_property(s, 10)
            end
        end
    end
end

# tests on example datasets
df = MVDA.dataset("iris")
target, X = Vector(df.target), Matrix(df[!,2:end])
k = 2

@testset "w/ Intercept" begin
    prob = MVDAProblem(target, X, intercept=true)
    test_on_dataset(prob, target, X, k)
end

@testset "w/o Intercept" begin
    prob = MVDAProblem(target, X, intercept=false)
    test_on_dataset(prob, target, X, k)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

# Sparse utils to make assembly of local and global matrices easier.
# Unoptimized but should do all necessary stuff for at start.

mutable struct SparseMatrixCOO{T<:Real}
    I :: Vector{Int}
    J :: Vector{Int}
    V :: Vector{T}
end

const SparseVectorCOO = SparseMatrixCOO

function SparseMatrixCOO()
    return SparseMatrixCOO{Float64}([], [], [])
end

function SparseVectorCOO(I::Vector, V::Vector)
    return SparseVectorCOO(I, ones(I), V)
end

function convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC)
    return SparseMatrixCOO(findnz(A)...)
end

function convert(::Type{SparseVectorCOO}, A::SparseVector)
    return SparseVectorCOO(findnz(A)...)
end

function convert(::Type{SparseMatrixCOO}, A::Matrix)
    return SparseMatrixCOO(findnz(A)...)
end

""" Convert from COO format to CSC.

Parameters
----------
tol
    used to drop near zero values less than tol.
"""
function sparse(A::SparseMatrixCOO; tol=1.0e-12)
    B = sparse(A.I, A.J, A.V)
    SparseArrays.droptol!(B, tol)
    return B
end

function sparse(A::SparseMatrixCOO, n::Int, m::Int; tol=1.0e-12)
    B = sparse(A.I, A.J, A.V, n, m)
    SparseArrays.droptol!(B, tol)
    return B
end

function sparse(A::SparseMatrixCOO, n::Int, m::Int, f::Function; tol=1.0e-12)
    B = sparse(A.I, A.J, A.V, n, m, f)
    SparseArrays.droptol!(B, tol)
    return B
end

# For backward compatibility, will be removed
function push!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64)
    push!(A.I, I)
    push!(A.J, J)
    push!(A.V, V)
end

function add!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64)
    push!(A.I, I)
    push!(A.J, J)
    push!(A.V, V)
end

function add!(A::SparseMatrixCOO, I::Int, V::Float64)
    push!(A.I, I)
    push!(A.J, 1)
    push!(A.V, V)
end

function empty!(A::SparseMatrixCOO)
    empty!(A.I)
    empty!(A.J)
    empty!(A.V)
end

function append!(A::SparseMatrixCOO, B::SparseMatrixCOO)
    append!(A.I, B.I)
    append!(A.J, B.J)
    append!(A.V, B.V)
end

function isempty(A::SparseMatrixCOO)
    return isempty(A.I) && isempty(A.J) && isempty(A.V)
end

function full(A::SparseMatrixCOO, args...)
    return full(sparse(A.I, A.J, A.V, args...))
end

""" Add local element matrix to sparse matrix. This basically does:

>>> A[dofs1, dofs2] = A[dofs1, dofs2] + data

Example
-------

>>> S = [3, 4]
>>> M = [6, 7, 8]
>>> data = Float64[5 6 7; 8 9 10]
>>> A = SparseMatrixCOO()
>>> add!(A, S, M, data)
>>> full(A)
4x8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
 0.0  0.0  0.0  0.0  0.0  5.0  6.0   7.0
 0.0  0.0  0.0  0.0  0.0  8.0  9.0  10.0

"""
function add!(A::SparseMatrixCOO, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix)
    n, m = size(data)
    for j=1:m
        for i=1:n
            add!(A, dofs1[i], dofs2[j], data[i,j])
        end
    end
    return nothing
end

""" Add sparse matrix of CSC to COO. """
function add!(A::SparseMatrixCOO, B::SparseMatrixCSC)
    i, j, v = findnz(B)
    C = SparseMatrixCOO(i, j, v)
    append!(A, C)
end

""" Add new data to COO Sparse vector. """
function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::Int=1)
    if length(dofs) != length(data)
        info("dofs = $dofs")
        info("data = $(vec(data))")
        error("when adding to sparse vector dimension mismatch!")
    end
    append!(A.I, dofs)
    append!(A.J, dim*ones(Int, length(dofs)))
    append!(A.V, vec(data))
end

""" Add SparseVector to SparseVectorCOO. """
function add!(a::SparseVectorCOO, b::SparseVector)
    i, v = findnz(b)
    c = SparseVectorCOO(i, v)
    append!(a, c)
    return
end

""" Combine (I,J,V) values if possible to reduce memory usage. """
function optimize!(A::SparseMatrixCOO)
    i, j, v = findnz(sparse(A))
    A.I = i
    A.J = j
    A.V = v
    return
end

""" Find all nonzero rows from sparse matrix.

Returns
-------

Ordered list of row indices.
"""
function get_nonzero_rows(A::SparseMatrixCSC)
    return sort(unique(rowvals(A)))
end

function get_nonzero_columns(A::SparseMatrixCSC)
    return get_nonzero_rows(transpose(A))
end

function size(A::SparseMatrixCOO)
    isempty(A) && return (0, 0)
    return maximum(A.I), maximum(A.J)
end

function size(A::SparseMatrixCOO, idx::Int)
    return size(A)[idx]
end

""" Resize sparse matrix A to (higher) dimension n x m. """
function resize_sparse(A, n, m)
    return sparse(findnz(A)..., n, m)
end

""" Resize sparse vector b to (higher) dimension n. """
function resize_sparsevec(b, n)
    return sparsevec(findnz(b)..., n)
end

""" Approximative comparison of two matricse A and B. """
function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO)
    A2 = sparse(A)
    B2 = sparse(B, size(A2)...)
    return isapprox(A2, B2)
end


text:
@testset "utils" begin

@testset "cloning" begin
    Context() do ctx
        # set-up
        mod = LLVM.Module("my_module"; ctx)

        param_types = [LLVM.Int32Type(ctx), LLVM.Int32Type(ctx)]
        ret_type = LLVM.Int32Type(ctx)
        fun_type = LLVM.FunctionType(ret_type, param_types)
        f = LLVM.Function(mod, "f", fun_type)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(f, "entry"; ctx)
            position!(builder, entry)

            tmp = add!(builder, parameters(f)[1], parameters(f)[2], "tmp")
            ret!(builder, tmp)

            verify(mod)
        end

        # basic clone
        let new_f = clone(f)
            @test new_f != f
            @test llvmtype(new_f) == llvmtype(f)
            for (bb1, bb2) in zip(blocks(f), blocks(new_f))
                for (inst1, inst2) in zip(instructions(bb2), instructions(bb2))
                    @test inst1 == inst2
                end
            end
        end

        # clone into, testing the value mapper (adding an argument)
        let
            new_param_types = [LLVM.Int32Type(ctx), LLVM.Int32Type(ctx), LLVM.Int32Type(ctx)]
            new_fun_type = LLVM.FunctionType(ret_type, new_param_types)
            new_f = LLVM.Function(mod, "new", new_fun_type)

            value_map = Dict{LLVM.Value, LLVM.Value}(
                parameters(f)[1] => parameters(new_f)[2],
                parameters(f)[2] => parameters(new_f)[3]
            )
            clone_into!(new_f, f; value_map)

            # operands of the add instruction should have been remapped
            add = first(instructions(first(blocks(new_f))))
            @test operands(add)[1] == parameters(new_f)[2]
            @test operands(add)[2] == parameters(new_f)[3]
        end

        # clone into, testing the type remapper (changing precision)
        let
            new_param_types = [LLVM.Int64Type(ctx), LLVM.Int64Type(ctx)]
            new_ret_type = LLVM.Int64Type(ctx)
            new_fun_type = LLVM.FunctionType(new_ret_type, new_param_types)
            new_f = LLVM.Function(mod, "new", new_fun_type)

            # we always need to map all arguments
            value_map = Dict{LLVM.Value, LLVM.Value}(
                old_param => new_param for (old_param, new_param) in
                                            zip(parameters(f), parameters(new_f)))

            function type_mapper(typ)
                if typ == LLVM.Int32Type(ctx)
                    LLVM.Int64Type(ctx)
                else
                    typ
                end
            end
            function materializer(val)
                # not needed here
                error()
            end
            clone_into!(new_f, f; value_map, type_mapper, materializer)

            # the add should now be a 64-bit addition
            add = first(instructions(first(blocks(new_f))))
            @test llvmtype(operands(add)[1]) == LLVM.Int64Type(ctx)
            @test llvmtype(operands(add)[2]) == LLVM.Int64Type(ctx)
            @test llvmtype(add) == LLVM.Int64Type(ctx)
        end
    end
end

end


text:
# Some useful wrappers
function get_function(model::MOI.ModelLike, ci::CI)
    return MOI.get(model, MOI.ConstraintFunction(), ci)
end

function get_set(model::MOI.ModelLike, ci::CI)
    return MOI.get(model, MOI.ConstraintSet(), ci)
end

function get_ci_row_dimension(model::MOI.ModelLike, ci::CI)
    return MOI.output_dimension(get_function(model, ci))
end

function get_scalar_term(
    model::MOI.ModelLike,
    ci::CI{SVF,S},
) where {S<:MOI.AbstractScalarSet}
    return [-MOI.constant(get_set(model, ci))]
end

function get_scalar_term(
    model::MOI.ModelLike,
    ci::CI{F,S},
) where {F<:MOI.AbstractScalarFunction,S<:MOI.AbstractScalarSet}
    return [
        MOI.constant(get_function(model, ci)) -
        MOI.constant(get_set(model, ci)),
    ]
end

# This is used to fill the dual objective dictionary
function get_scalar_term(
    model::MOI.ModelLike,
    i::Int,
    ci::CI{F,S},
) where {F<:MOI.AbstractVectorFunction,S<:MOI.AbstractVectorSet}
    return MOI.constant(get_function(model, ci))[i]
end

# This is used to fill the dual objective dictionary
function get_scalar_term(
    model::MOI.ModelLike,
    i::Int,
    ci::CI{SVF,S},
) where {S<:MOI.AbstractScalarSet}
    return -MOI.constant(get_set(model, ci))
end

# This is used to fill the dual objective dictionary
function get_scalar_term(
    model::MOI.ModelLike,
    i::Int,
    ci::CI{F,S},
) where {F<:MOI.AbstractScalarFunction,S<:MOI.AbstractScalarSet}

    # In this case there i only one constant in the function and one in the set
    return MOI.constant(get_function(model, ci))[1] -
           MOI.constant(get_set(model, ci))
end


text:
function _coefnames(x; show_intercept=true)
    cn = coefnames(x)
    return show_intercept ? cn : filter!(!=("(Intercept)"), cn)
end

function _coefnames(x::MixedModelBootstrap; show_intercept=true)
    cn = collect(string.(propertynames(first(x.fits).β)))
    return show_intercept ? cn : filter!(!=("(Intercept)"), cn)
end

"""
    confint_table(x::StatsBase.StatisticalModel, level=0.95)
    confint_table(x::MixedModelBootstrap, level=0.95)

Return a DataFrame of coefficient names, point estimates and confidence intervals.

`level` specifies the confidence level, e.g. `0.95` for 95% confidence intervals.

For `MixedModels`, the intervals are computed using the standard errors and the asymptotic
Wald approximation (e.g. est±1.96*se for 95% intervals). For `MixedModelBootstrap`, the intervals
are computed using `shortestcovint`, and the point estimate re-computed as the mean of the bootstrap values.

The returned table has the following columns:
- `coefname`: the names of the coefficients
- `estimate`: the point estimates
- `lower`: the lower edge of the confidence interval
- `upper`: the upper edge of the confidence interval

"""
function confint_table(x::StatsBase.StatisticalModel, level=0.95)
    # taking from the lower tail
    semultiple = zquantile((1 - level) / 2)
    se = stderror(x)

    return DataFrame(;
        coefname=coefnames(x),
        estimate=coef(x),
        lower=coef(x) + semultiple * se,
        upper=coef(x) - semultiple * se,
    )
end

function confint_table(x::MixedModelBootstrap, level=0.95)
    df = transform!(
        select!(DataFrame(x.β), Not(:iter)), :coefname => ByRow(string) => :coefname
    )
    ci(x) = shortestcovint(x, level)
    return combine(
        groupby(df, :coefname),
        :β => mean => :estimate,
        :β => NamedTuple{(:lower, :upper)} ∘ ci => [:lower, :upper],
    )
end

_npreds(x; show_intercept=true) = length(_coefnames(x; show_intercept))

"""
    ppoints(n::Integer)

Return a sequence of `n` equally-spaced points in the interval (0, 1) - so-called "probability points"
"""
ppoints(n::Integer) = inv(2n):inv(n):1

"""
    zquantile(x::AbstractFloat)

Evaluate `quantile(Normal(), x)` using only the `SpecialFunctions` package (i.e. not requiring `Distributions`).
"""
zquantile(x::T) where {T<:AbstractFloat} = -erfcinv(2x) * sqrt(T(2))


text:
function puzzle03(input::Int=325489)
    L = iseven(Int(ceil(√input))) ? ceil(√input) + 1 : ceil(√input)
    corners = L^2 .- (0:3) .* (L-1)
    width = floor(L/2)
    height = minimum(abs.((input + width) .- corners))
    partone = Int(height + width)

    # initialize matrix
    # (x, y, v)
    coords = [(0, 0, 1)]
    direction = ["right", "up", "left", "down"]
    steps = 1
    notfound = true
    while notfound
        for s in 1:steps
            x, y, v = last(coords)
            if first(direction) == "right"
                x += 1
            elseif first(direction) == "up"
                y += 1
            elseif first(direction) == "left"
                x -= 1
            elseif first(direction) == "down"
                y -= 1
            end
            neighbors = reshape([(x+xx,y+yy) for xx=-1:1, yy=-1:1], 9)
            v = sum([x[3] for x in coords if x[1:2] in neighbors])
            append!(coords, [(x, y, v)])
            if v > input
                notfound = false
                break
            end
        end
        direction = circshift(direction, -1)
        if first(direction) in ["left", "right"]
            steps += 1
        end
    end
    parttwo = last(coords)[3]

    return [partone, parttwo]
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion between Euler angles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angle_to_angle.jl
# =========================================

# Functions: angle_to_angle
# -------------------------

@testset "Euler angles => Euler angles" begin
    for T in (Float32, Float64)
        # We do not need comprehensive test here because `angle_to_angle` first
        # converts the Euler angles to DCM and then back to Euler angles. Those
        # two operations are already heavily tested.

        ea1 = EulerAngles(deg2rad(T(45)), 0, 0, :XYZ)
        ea2 = angle_to_angle(ea1, :ZYX)
        @test eltype(ea2) === T
        @test ea2.a1 ≈ T(0)
        @test ea2.a2 ≈ T(0)
        @test ea2.a3 ≈ deg2rad(T(45))
        @test ea2.rot_seq == :ZYX

        ea1 = EulerAngles(deg2rad(T(45)), 0, 0, :YXZ)
        ea2 = angle_to_angle(ea1, :ZYX)
        @test eltype(ea2) === T
        @test ea2.a1 ≈ T(0)
        @test ea2.a2 ≈ deg2rad(T(45))
        @test ea2.a3 ≈ T(0)
        @test ea2.rot_seq == :ZYX

        ea1 = EulerAngles(deg2rad(T(45)), 0, 0, :ZXY)
        ea2 = angle_to_angle(ea1, :ZYX)
        @test eltype(ea2) === T
        @test ea2.a1 ≈ deg2rad(T(45))
        @test ea2.a2 ≈ T(0)
        @test ea2.a3 ≈ T(0)
        @test ea2.rot_seq == :ZYX

        ea1 = EulerAngles(T(π / 2), T(π / 6), T(2 / 3 * π), :ZYX)
        ea2 = angle_to_angle(ea1, :ZYX)
        @test eltype(ea2) === T
        @test ea2.a1 ≈ ea1.a1
        @test ea2.a2 ≈ ea1.a2
        @test ea2.a3 ≈ ea1.a3
        @test ea2.rot_seq == :ZYX
    end
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Implementation of the geometry required by the renderer (manipulation of 3D scenes, cameras, rays, ...)

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

"""
    Vec <: StaticArrays.FieldVector{3, Float32}

A vector in 3D space with 3 fields `x`, `y`, and `z` of type `Float32`.

For inherited properties and constructors see `StaticArrays.FieldVector`.
"""
struct Vec <: FieldVector{3, Float32}
    x::Float32
    y::Float32
    z::Float32
end

@doc """
    Vec(x::Float32, y::Float32, z::Float32)

Constructor for a [`Vec`](@ref) instance.
""" Vec(::Float32, ::Float32, ::Float32)

similar_type(::Type{<:Vec}, ::Type{Float32}, s::Size{(3,)}) = Vec

################
# Miscellaneous

function show(io::IO, a::Vec)
    print(io, typeof(a), "(", join((string(el) for el ∈ a), ", "), ")")
end

function show(io::IO, ::MIME"text/plain", a::Vec)
    print(io, "Vec with eltype $(eltype(a))\n", join(("$label = $el" for (label, el) ∈ zip((:x, :y, :z), a)), ", "))
end

"""
    norm²(v::Vec)

Compute the squared norm of a [`Vec`](@ref).

# Examples

```jldoctest
julia> norm²(Vec(1, 2, 3))
14.0f0
```
"""
norm²(v::Vec) = sum(el -> el^2, v)

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

"""
    Normal{V} <: StaticArrays.FieldVector{3, Float32}

A pseudo-vector in 3D space with 3 fields `x`, `y`, and `z` of type `Float32`.
The parameter `V` tells if the normal is normalized or not.

For inherited properties and constructors see `StaticArrays.FieldVector`.
"""
struct Normal{V} <: FieldVector{3, Float32}
    x::Float32
    y::Float32
    z::Float32
end

@doc """
    Normal{V}(x::Float32, y::Float32, z::Float32)

Constructor for a [`Normal`](@ref) instance.
""" Normal{V}(::Float32, ::Float32, ::Float32)

"""
    Normal(x, y, z)

Construct a non-normalized [`Normal{false}`](@ref) with given coordinates. All values are converted in `Float32`.

# Examples

```jldoctest
julia> Normal(1.2, 3.3, 5)
Normal with eltype Float32, not normalized
x = 1.2, y = 3.3, z = 5.0
```
"""
Normal(x, y, z) = Normal{false}(map(x -> convert(Float32, x), (x, y, z)))

similar_type(::Type{<:Normal}, ::Type{Float32}, s::Size{(3,)}) = Normal

################
# Miscellaneous

function show(io::IO, a::Normal)
    print(io, typeof(a), "(", join((string(el) for el ∈ a), ", "), ")")
end

function show(io::IO, ::MIME"text/plain", n::Normal{V}) where {V}
    print(io, "Normal with eltype $(eltype(n)), ", V ? "normalized" : "not normalized", "\n",
          join(("$label = $el" for (label, el) ∈ zip((:x, :y, :z), n)), ", "))
end

"""
    normalize(n::Normal)

Normalize `n` and return a [`Normal{true}`](@ref). If `n` is already a `Normal{true}` instance, no normalization is computed and `n` is returned.

# Examples

```jldoctest
julia> n = normalize(Normal(1,2,4))
Normal with eltype Float32, normalized
x = 0.21821788, y = 0.43643576, z = 0.8728715

julia> normalize(n)
Normal with eltype Float32, normalized
x = 0.21821788, y = 0.43643576, z = 0.8728715
```
"""
normalize(n::Normal{false}) = Normal{true}(normalize(SVector{3}(n)))
normalize(n::Normal{true}) = n

"""
    norm(n::Normal{true})

Compute the squared norm of a [`Normal{true}`](@ref). Since `n` is already normalized, `1f0` is returned.
"""
norm(::Normal{true}) = 1f0

"""
    norm²(n::Normal)

Compute the squared norm of a [`Normal`](@ref). If `n` is a `Normal{true}` instance then `1f0` is returned.

# Examples

```jldoctest
julia> n = Normal(1, 2, 3);

julia> norm²(n)
14.0f0

julia> norm²(normalize(n))
1.0f0
```
"""
norm²(n::Normal{false}) = sum(el -> el^2, n)
norm²(::Normal{true}) = 1f0

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

"""
    Point

A point in a 3D space. Implemented as a wrapper struct around a `SVector{3, Float32}`.
"""
struct Point
    v::SVector{3, Float32}

    function Point(p::AbstractVector)
        size(p) == (3,) || throw(ArgumentError("argument 'p' has size = $(size(p)) but 'Point' requires an argument of size = (3,)"))
        new(SVector{3}(p))
    end
end

@doc """
    Point(p::AbstractVector)

Constructor for a [`Point`](@ref) instance.
""" Point(::AbstractVector)

"""
    Point(x, y, z)

Construct a [`Point`](@ref) with given coordinates. All values are converted in `Float32`.

# Examples

```jldoctest
julia> Point(1.2, 3.3, 5)
Point with eltype Float32
x = 1.2, y = 3.3, z = 5.0
```
"""
Point(x, y, z) = Point(SVector(map(x -> convert(Float32, x), (x, y, z))))

################
# Miscellaneous

eltype(::Point) = Float32
eltype(::Type{Point}) = Float32

length(::Point) = 3
firstindex(::Point) = 1
lastindex(p::Point) = 3
getindex(p::Point, i::Integer) = p.v[i]
iterate(p::Point, state = 1) = state > 3 ? nothing : (p[state], state +1)

function show(io::IO, p::Point)
    print(io, typeof(p), "(", join((string(el) for el ∈ p), ", "), ")")
end

function show(io::IO, ::MIME"text/plain", p::Point)
    print(io, Point, " with eltype $(eltype(p))\n", join(("$label = $el" for (label, el) ∈ zip((:x, :y, :z), p)), ", "))
end

"""
    ≈(p1::Point, p2::Point)

Check if two points are close.

# Examples

```jldoctest
julia> p = Point(1, 2, 3);

julia> p ≈ Point(1, 2, 3)
true

julia> p ≈ Point(0, 0, 0)
false
```
"""
(≈)(p1::Point, p2::Point) = p1.v ≈ p2.v

"""
    -(p1::Point, p2::Point)

Return the elementwise difference of two [`Point`](@ref) as an instance of [`Vec`](@ref).

# Examples

```jldoctest
julia> Point(1, 2, 3) - Point(4, 5, 6)
Vec with eltype Float32
x = -3.0, y = -3.0, z = -3.0
```
"""
(-)(p1::Point, p2::Point) = Vec(p1.v - p2.v)

"""
    -(p::Point, v::Vec)

Return the elementwise difference between a [`Point`](@ref) and a [`Vec`](@ref) as an instance of `Point`.

# Examples

```jldoctest
julia> Point(1, 2, 3) - Vec(4, 5, 6)
Point with eltype Float32
x = -3.0, y = -3.0, z = -3.0
```
"""
(-)(p::Point, v::Vec) = Point(p.v - v)

"""
    +(p::Point, v::Vec)

Return the elementwise sum between a [`Point`](@ref) and a [`Vec`](@ref) as an instance of `Point`.

# Examples

```jldoctest
julia> Point(1, 2, 3) + Vec(4, 5, 6)
Point with eltype Float32
x = 5.0, y = 7.0, z = 9.0
```
"""
(+)(p::Point, v::Vec) = Point(p.v + v)

"""
    *(p::Point, s...)

Multiplication operator. `x * y * z * ...`` calls this function with all arguments, i.e. `*(x, y, z, ...)`.

Return a [`Point`](@ref).

# Examples

```jldoctest
julia> Point(1, 2, 3) * 2 * 3
Point with eltype Float32
x = 6.0, y = 12.0, z = 18.0
```
"""
(*)(p::Point, s...) = (*)(p.v, s...) |> Point

"""
    convert(::Type{Vec}, p::Point)
    convert(::Type{Normal}, p::Point)

Convert a [`Point`](@ref) into the specified type ([`Vec`](@ref) or [`Normal{false}`](@ref)).
"""
convert(::Type{Vec}, p::Point) = Vec(p.v)
convert(::Type{Normal}, p::Point) = Normal{false}(p.v)

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

"""
    Vec2D

Alias to `SVector{2, Float32}`, used for uv mapping on shapes.
"""
const Vec2D = SVector{2, Float32}

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

"""
    VEC_X

A unitary [`Vec`](@ref) along the x-axis.
"""
const VEC_X = Vec(1f0, 0f0, 0f0)

"""
    VEC_Y

A unitary [`Vec`](@ref) along the y-axis.
"""
const VEC_Y = Vec(0f0, 1f0, 0f0)

"""
    VEC_Z

A unitary [`Vec`](@ref) along the z-axis.
"""
const VEC_Z = Vec(0f0, 0f0, 1f0)

"""
    NORMAL_X

A unitary and normalized [`Normal{true}`](@ref) along the x-axis.
"""
const NORMAL_X = Normal{true}(1f0, 0f0, 0f0)

"""
    NORMAL_Y

A unitary and normalized [`Normal{true}`](@ref) along the y-axis.
"""
const NORMAL_Y = Normal{true}(0f0, 1f0, 0f0)

"""
    NORMAL_Z

A unitary and normalized [`Normal{true}`](@ref) along the z-axis.
"""
const NORMAL_Z = Normal{true}(0f0, 0f0, 1f0)

"""
    NORMAL_X_false

A unitary and non-normalized [`Normal{false}`](@ref) along the x-axis.
"""
const NORMAL_X_false = Normal{false}(1f0, 0f0, 0f0)

"""
    NORMAL_Y_false

A unitary and non-normalized [`Normal{false}`](@ref) along the y-axis.
"""
const NORMAL_Y_false = Normal{false}(0f0, 1f0, 0f0)

"""
    NORMAL_Z_false

A unitary and non-normalized [`Normal{false}`](@ref) along the z-axis.
"""
const NORMAL_Z_false = Normal{false}(0f0, 0f0, 1f0)

"""
    ORIGIN

A [`Point`](@ref) representing the origin of the frame of reference.
"""
const ORIGIN = Point(0f0, 0f0, 0f0)

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

"""
    create_onb_from_z(input_normal::Normal)

Create an orthonormal base from a input [`Normal`](@ref) representing the z-axis using the
[Duff et al. 2017](https://graphics.pixar.com/library/OrthonormalB/paper.pdf) algorithm.

# Examples

```jldoctest
julia> n = Normal(0,0,5);

julia> nn = normalize(Normal(0,0,5));

julia> create_onb_from_z(n)
(Vec(1.0, -0.0, -0.0), Vec(-0.0, 1.0, -0.0), Vec(0.0, 0.0, 1.0))

julia> create_onb_from_z(nn)
(Vec(1.0, -0.0, -0.0), Vec(-0.0, 1.0, -0.0), Vec(0.0, 0.0, 1.0))
```

Note that `create_onb_from_z(n)` and `create_onb_from_z(nn)` give the same result.
"""
function create_onb_from_z(input_normal::Normal)
    normal = normalize(input_normal)

    sign = copysign(1f0, normal.z)

    a = -1f0 / (sign + normal.z)
    b = normal.x * normal.y * a

    e1 = Vec(1f0 + sign * normal.x * normal.x * a, sign * b, -sign * normal.x)
    e2 = Vec(b, sign + normal.y * normal.y * a, -normal.y)
    e3 = convert(Vec, normal)

    (e1, e2, e3)
end

"""
    normalized_dot(v1::AbstractVector, v2::AbstractVector)

Normalize `v1` and `v2` and then compute the dot product.
"""
normalized_dot(v1::AbstractVector, v2::AbstractVector) = normalize(v1) ⋅ normalize(v2)


text:
"""
    accumulated_edges(adj)

Return a vector which acts as a mapping table. The index is the vertex index,
value is accumulated numbers of edge (current vertex not included).
"""
accumulated_edges(adj::AbstractVector{<:AbstractVector{<:Integer}}) = [0, cumsum(map(length, adj))...]

function generate_cluster(M::AbstractArray{T,N}, accu_edge) where {T,N}
    num_V = length(accu_edge) - 1
    num_E = accu_edge[end]
    cluster = similar(M, Int, num_E)
    @inbounds for i = 1:num_V
        j = accu_edge[i]
        k = accu_edge[i+1]
        cluster[j+1:k] .= i
    end
    cluster
end

@non_differentiable accumulated_edges(x...)
@non_differentiable generate_cluster(x...)


text:
mutable struct VerticalDiffusion

    amo :: AdvancedMatrixOperators

    eT_length    :: Int64

    eT_send_T    :: AbstractArray{Float64, 2}
    T_send_eT    :: AbstractArray{Float64, 2}
    eT_I_eT      :: AbstractArray{Float64, 2}

    K_iso :: Float64
    K_cva :: Float64

    K_func :: Function

    wksp

    function VerticalDiffusion(
        amo   :: AdvancedMatrixOperators;
        K_iso :: Float64,
        K_cva :: Float64,
    )

        cvtDiagOp = (a,) -> spdiagm(0 => view(a, :))
        Nz, Ny, Nx = amo.bmo.T_dim

        # Create coversion matrix and its inverse
        T_num        = reshape(amo.T_mask_T * collect(1:amo.bmo.T_pts), Nz, Ny, Nx)
        active_T_num = T_num[ T_num .!= 0.0 ]
        eT_send_T    = amo.bmo.T_I_T[active_T_num, :]
        T_send_eT    = eT_send_T' |> sparse
        eT_I_eT      = eT_send_T * T_send_eT

        eT_length = size(eT_I_eT)[1]

        #K = spdiagm(0 => ones(Float64, amo.bmo.W_pts))

        K_func = (dbdz,) -> (dbdz > 0.0) ? K_iso : K_cva

        # iEBD = inverse of Euler Backward Diffusion ( I - Δt * ∇K∇ )
        wksp = (
            above_W    = zeros(Float64, amo.bmo.W_pts),
            below_W    = zeros(Float64, amo.bmo.W_pts),
            Δ_W        = zeros(Float64, amo.bmo.W_pts),
            K_W        = zeros(Float64, amo.bmo.W_pts),
            W_K_W      = spdiagm(0 => ones(Float64, amo.bmo.W_pts)),
            W_K∇_T     = spzeros(Float64, amo.bmo.W_pts, amo.bmo.T_pts),
            T_∇K∇_T    = spzeros(Float64, amo.bmo.T_pts, amo.bmo.T_pts),
            eT_∇K∇_T   = spzeros(Float64, eT_length, amo.bmo.T_pts),
            eT_∇K∇_eT  = spzeros(Float64, eT_length, eT_length),
            eT_iEBD_eT = spzeros(Float64, eT_length, eT_length),
            rhs_eT     = zeros(Float64, eT_length),
            lhs_eT     = zeros(Float64, eT_length),
        )

        return new(
            amo,
            eT_length,
            eT_send_T,
            T_send_eT,
            eT_I_eT,
            K_iso,
            K_cva,
            K_func,
            wksp,
        )

    end
end

function calOp_vdiff(
    vd   :: VerticalDiffusion,
    b    :: AbstractArray, # of grid :T
    HMXL :: AbstractArray, # of grid :sT
)

    amo = vd.amo
    bmo = amo.bmo
    gd  = amo.gd

    dbdz = vd.amo.W_∂z_T * view(b, :)
    HMXL = reshape(HMXL, 1, gd.Nx, gd.Ny)

    K_W = reshape( vd.K_func.(dbdz), bmo.W_dim...)
    K_W[gd.z_W .> - HMXL] .= vd.K_cva

    # below this W-pts the temperature is below freezing point.
    # reasoning please see stepColumn discussing about Q_FRZHEAT / frzmlt
    #below_frz_below = (bmo.W_UP_T * view(TEMP, :)) .< T_sw_frz
    #K_W[below_frz_below] .= vd.K_cva

    #K_W .= 0
    op_vdiff = sparse(vd.amo.T_DIVz_W * vd.amo.W_mask_W * spdiagm( 0 => view(K_W,:)) * vd.amo.W_∂z_T)

    dropzeros!(op_vdiff)

    return op_vdiff
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using SpaceExParser
using SpaceExParser: _get_coeffs
using MathematicalSystems
using MAT

H = readsxmodel("heli_modif.xml", raw_dict=true, ST=ConstrainedLinearControlContinuousSystem)

state_vars = convert.(SpaceExParser.Basic, collect(keys(H["variables"]))[1:28])
input_vars = convert.(SpaceExParser.Basic, collect(keys(H["variables"]))[29:34]);

A, B = _get_coeffs(H["flows"][1], 28, 6, state_vars, input_vars);

matwrite("heli.mat", Dict("A" => A, "B" => B))

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary