text:
type MMatrix{S1, S2, T, L} <: StaticMatrix{T}
    data::NTuple{L, T}

    function MMatrix(d::NTuple{L,T})
        check_MMatrix_params(Val{S1}, Val{S2}, T, Val{L})
        new(d)
    end

    function MMatrix(d::NTuple{L})
        check_MMatrix_params(Val{S1}, Val{S2}, T, Val{L})
        new(convert_ntuple(T, d))
    end

    function MMatrix()
        check_MMatrix_params(Val{S1}, Val{S2}, T, Val{L})
        new()
    end
end

@generated function check_MMatrix_params{S1,S2,L}(::Type{Val{S1}}, ::Type{Val{S2}}, T, ::Type{Val{L}})
    if !(T <: DataType) # I think the way types are handled in generated fnctions might have changed in 0.5?
        return :(error("MMatrix: Parameter T must be a DataType. Got $T"))
    end

    if !isa(S1, Int) || !isa(S2, Int) || !isa(L, Int) || S1 < 0 || S2 < 0 || L < 0
        return :(error("MMatrix: Sizes must be positive integers. Got $S1 × $S2 ($L elements)"))
    end

    if S1*S2 == L
        return nothing
    else
        str = "Size mismatch in MMatrix. S1 = $S1, S2 = $S2, but recieved $L elements"
        return :(error($str))
    end
end

@generated function (::Type{MMatrix{S1}}){S1,L}(x::NTuple{L})
    S2 = div(L, S1)
    if S1*S2 != L
        error("Incorrect matrix sizes. $S1 does not divide $L elements")
    end
    T = promote_tuple_eltype(x)

    return quote
        $(Expr(:meta, :inline))
        MMatrix{S1, $S2, $T, L}(x)
    end
end

@generated function (::Type{MMatrix{S1,S2}}){S1,S2,L}(x::NTuple{L})
    T = promote_tuple_eltype(x)

    return quote
        $(Expr(:meta, :inline))
        MMatrix{S1, S2, $T, L}(x)
    end
end

@generated function (::Type{MMatrix{S1,S2,T}}){S1,S2,T,L}(x::NTuple{L})
    return quote
        $(Expr(:meta, :inline))
        MMatrix{S1, S2, T, L}(x)
    end
end

@generated function (::Type{MMatrix{S1,S2,T}}){S1,S2,T}()
    return quote
        $(Expr(:meta, :inline))
        MMatrix{S1, S2, T, $(S1*S2)}()
    end
end

@inline convert{S1,S2,T}(::Type{MMatrix{S1,S2}}, a::StaticArray{T}) = MMatrix{S1,S2,T}(Tuple(a))
@inline convert{S1,S2,T}(::Type{MMatrix{S1,S2}}, a::AbstractArray{T}) = MMatrix{S1,S2,T}((a...))
@inline MMatrix(a::StaticMatrix) = MMatrix{size(typeof(a),1),size(typeof(a),2)}(Tuple(a))

# Some more advanced constructor-like functions
@inline eye{N}(::Type{MMatrix{N}}) = eye(MMatrix{N,N})
@inline eye{N,M}(::Type{MMatrix{N,M}}) = eye(MMatrix{N,M,Float64})
@inline zeros{N,M}(::Type{MMatrix{N,M}}) = zeros(MMatrix{N,M,Float64})
@inline ones{N,M}(::Type{MMatrix{N,M}}) = ones(MMatrix{N,M,Float64})

#####################
## MMatrix methods ##
#####################

@pure size{S1,S2}(::Union{MMatrix{S1,S2},Type{MMatrix{S1,S2}}}) = (S1, S2)
@pure size{S1,S2,T}(::Type{MMatrix{S1,S2,T}}) = (S1, S2)
@pure size{S1,S2,T,L}(::Type{MMatrix{S1,S2,T,L}}) = (S1, S2)

@propagate_inbounds function getindex{S1,S2,T}(m::MMatrix{S1,S2,T}, i::Integer)
    #@boundscheck if i < 1 || i > length(m)
    #    throw(BoundsError(m,i))
    #end

    # This is nasty... but it turns out Julia will literally copy the whole tuple to the stack otherwise!
    if isbits(T)
        unsafe_load(Base.unsafe_convert(Ptr{T}, Base.data_pointer_from_objref(m)), i)
    else
        # Not sure about this... slow option for now...
        m.data[i]
        #unsafe_load(Base.unsafe_convert(Ptr{Ptr{Void}}, Base.data_pointer_from_objref(m.data)), i)
    end
end

@propagate_inbounds setindex!{S1,S2,T}(m::MMatrix{S1,S2,T}, val, i::Integer) = setindex!(m, convert(T, val), i)
@propagate_inbounds function setindex!{S1,S2,T}(m::MMatrix{S1,S2,T}, val::T, i::Integer)
    #@boundscheck if i < 1 || i > length(m)
    #    throw(BoundsError(m,i))
    #end

    if isbits(T)
        unsafe_store!(Base.unsafe_convert(Ptr{T}, Base.data_pointer_from_objref(m)), val, i)
    else # TODO check that this isn't crazy. Also, check it doesn't cause problems with GC...
        # This one is thought to be unsafe (#27)
        # unsafe_store!(Base.unsafe_convert(Ptr{Ptr{Void}}, Base.data_pointer_from_objref(m.data)), Base.data_pointer_from_objref(val), i)
        error("setindex!() with non-isbits eltype is not supported by StaticArrays")
    end

    return val
end

@inline Tuple(v::MMatrix) = v.data

@inline function Base.unsafe_convert{N,M,T}(::Type{Ptr{T}}, m::MMatrix{N,M,T})
    Base.unsafe_convert(Ptr{T}, Base.data_pointer_from_objref(m))
end

macro MMatrix(ex)
    if !isa(ex, Expr)
        error("Bad input for @MMatrix")
    end
    if ex.head == :vect && length(ex.args) == 1 # 1 x 1
        return esc(Expr(:call, MMatrix{1, 1}, Expr(:tuple, ex.args[1])))
    elseif ex.head == :ref && length(ex.args) == 2 # typed, 1 x 1
        return esc(Expr(:call, Expr(:curly, :MMatrix, 1, 1, ex.args[1]), Expr(:tuple, ex.args[2])))
    elseif ex.head == :hcat # 1 x n
        s1 = 1
        s2 = length(ex.args)
        return esc(Expr(:call, MMatrix{s1, s2}, Expr(:tuple, ex.args...)))
    elseif ex.head == :typed_hcat # typed, 1 x n
        s1 = 1
        s2 = length(ex.args) - 1
        return esc(Expr(:call, Expr(:curly, :MMatrix, s1, s2, ex.args[1]), Expr(:tuple, ex.args[2:end]...)))
    elseif ex.head == :vcat
        if isa(ex.args[1], Expr) && ex.args[1].head == :row # n x m
            # Validate
            s1 = length(ex.args)
            s2s = map(i -> ((isa(ex.args[i], Expr) && ex.args[i].head == :row) ? length(ex.args[i].args) : 1), 1:s1)
            s2 = minimum(s2s)
            if maximum(s2s) != s2
                error("Rows must be of matching lengths")
            end

            exprs = [ex.args[i].args[j] for i = 1:s1, j = 1:s2]
            return esc(Expr(:call, MMatrix{s1, s2}, Expr(:tuple, exprs...)))
        else # n x 1
            return esc(Expr(:call, MMatrix{length(ex.args), 1}, Expr(:tuple, ex.args...)))
        end
    elseif ex.head == :typed_vcat
        if isa(ex.args[2], Expr) && ex.args[2].head == :row # typed, n x m
            # Validate
            s1 = length(ex.args) - 1
            s2s = map(i -> ((isa(ex.args[i+1], Expr) && ex.args[i+1].head == :row) ? length(ex.args[i+1].args) : 1), 1:s1)
            s2 = minimum(s2s)
            if maximum(s2s) != s2
                error("Rows must be of matching lengths")
            end

            exprs = [ex.args[i+1].args[j] for i = 1:s1, j = 1:s2]
            return esc(Expr(:call, Expr(:curly, :MMatrix,s1, s2, ex.args[1]), Expr(:tuple, exprs...)))
        else # typed, n x 1
            return esc(Expr(:call, Expr(:curly, :MMatrix, length(ex.args)-1, 1, ex.args[1]), Expr(:tuple, ex.args[2:end]...)))
        end
    elseif isa(ex, Expr) && ex.head == :comprehension
        if length(ex.args) != 1 || !isa(ex.args[1], Expr) || ex.args[1].head != :generator
            error("Expected generator in comprehension, e.g. [f(i,j) for i = 1:3, j = 1:3]")
        end
        ex = ex.args[1]
        if length(ex.args) != 3
            error("Use a 2-dimensional comprehension for @MMatrx")
        end

        rng1 = eval(current_module(), ex.args[2].args[2])
        rng2 = eval(current_module(), ex.args[3].args[2])
        f = gensym()
        f_expr = :($f = (($(ex.args[2].args[1]), $(ex.args[3].args[1])) -> $(ex.args[1])))
        exprs = [:($f($j1, $j2)) for j1 in rng1, j2 in rng2]

        return quote
            $(Expr(:meta, :inline))
            $(esc(f_expr))
            $(esc(Expr(:call, Expr(:curly, :MMatrix, length(rng1), length(rng2)), Expr(:tuple, exprs...))))
        end
    elseif isa(ex, Expr) && ex.head == :typed_comprehension
        if length(ex.args) != 2 || !isa(ex.args[2], Expr) || ex.args[2].head != :generator
            error("Expected generator in typed comprehension, e.g. Float64[f(i,j) for i = 1:3, j = 1:3]")
        end
        T = ex.args[1]
        ex = ex.args[2]
        if length(ex.args) != 3
            error("Use a 2-dimensional comprehension for @MMatrx")
        end

        rng1 = eval(current_module(), ex.args[2].args[2])
        rng2 = eval(current_module(), ex.args[3].args[2])
        f = gensym()
        f_expr = :($f = (($(ex.args[2].args[1]), $(ex.args[3].args[1])) -> $(ex.args[1])))
        exprs = [:($f($j1, $j2)) for j1 in rng1, j2 in rng2]

        return quote
            $(Expr(:meta, :inline))
            $(esc(f_expr))
            $(esc(Expr(:call, Expr(:curly, :MMatrix, length(rng1), length(rng2), T), Expr(:tuple, exprs...))))
        end
    elseif isa(ex, Expr) && ex.head == :call
        if ex.args[1] == :zeros || ex.args[1] == :ones || ex.args[1] == :rand || ex.args[1] == :randn
            if length(ex.args) == 3
                return quote
                    $(Expr(:meta, :inline))
                    $(ex.args[1])(MMatrix{$(esc(ex.args[2])),$(esc(ex.args[3]))})
                end
            elseif length(ex.args) == 4
                return quote
                    $(Expr(:meta, :inline))
                    $(ex.args[1])(MMatrix{$(esc(ex.args[3])), $(esc(ex.args[4])), $(esc(ex.args[2]))})
                end
            else
                error("@MMatrix expected a 2-dimensional array expression")
            end
        elseif ex.args[1] == :fill
            if length(ex.args) == 4
                return quote
                    $(Expr(:meta, :inline))
                    $(esc(ex.args[1]))($(esc(ex.args[2])), MMatrix{$(esc(ex.args[3])), $(esc(ex.args[4]))})
                end
            else
                error("@MMatrix expected a 2-dimensional array expression")
            end
        elseif ex.args[1] == :eye
            if length(ex.args) == 2
                return quote
                    $(Expr(:meta, :inline))
                    eye(MMatrix{$(esc(ex.args[2]))})
                end
            elseif length(ex.args) == 3
                # We need a branch, depending if the first argument is a type or a size.
                return quote
                    $(Expr(:meta, :inline))
                    if isa($(esc(ex.args[2])), DataType)
                        eye(MMatrix{$(esc(ex.args[3])), $(esc(ex.args[3])), $(esc(ex.args[2]))})
                    else
                        eye(MMatrix{$(esc(ex.args[2])), $(esc(ex.args[3]))})
                    end
                end
            elseif length(ex.args) == 4
                return quote
                    $(Expr(:meta, :inline))
                    eye(MMatrix{$(esc(ex.args[3])), $(esc(ex.args[4])), $(esc(ex.args[2]))})
                end
            else
                error("Bad eye() expression for @MMatrix")
            end
        else
            error("@MMatrix only supports the zeros(), ones(), fill(), rand(), randn() and eye() functions.")
        end
    else
        error("Bad input for @MMatrix")
    end
end


text:
struct NICER <: Mission end

const nicer = NICER()

export nicer

_mission_name(::NICER) = "nicer"

_mission_master_url(::NICER) = "https://heasarc.gsfc.nasa.gov/FTP/heasarc/dbase/tdat_files/heasarc_nicermastr.tdat.gz"

function _obs_path_server(::NICER, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame})
    obsid      = obs_row[:obsid]
    date_time  = string(obs_row[:time])
    date_year  = date_time[1:4]
    date_month = date_time[6:7]

    folder_path = string("/.nicer_archive/.nicer_$(date_year)$(date_month)a/",
                            "obs/$(date_year)_$(date_month)/$obsid")

    return folder_path
end

function _obs_files_cl(::NICER, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame})
    cl_dir = joinpath(_obs_path_local(nicer, obs_row; kind=:download), "xti/event_cl/")

    cl_files = Dict(:XTI=>joinpath(cl_dir, "ni$(obs_row[:obsid])_0mpu7_cl.evt.gz"))

    return cl_files
end

_mission_good_e_range(::NICER) = (0.2, 12.0)

_mission_instruments(::NICER) = [:XTI]

ext: jl


text:
using Scat
using BenchmarkTools
using Plots

function time_subsamplingSize(Qs::Vector{Int64}, N::Int64; nTrials=1)
    # define the scattering transform
    # first layer parameters and filter bank
    J=8
    σ0=0.1

    times_subsamp = zeros(N)
    times_nosubsamp = zeros(N)
    for i=1:nTrials
        for n=1:N
            println(n)
            x=ones(2^n)
            Transform = ScatteringTransform1d(2, Qs, [J,J])
            FilterBanks = get_FilterBanks(n, Qs, J, σ0)
            times_subsamp[n] += @elapsed Scat.scatter1d(x, Transform, FilterBanks, subsample=true)
            times_nosubsamp[n] += @elapsed Scat.scatter1d(x, Transform, FilterBanks, subsample=false)
        end
    end
    times_subsamp./=nTrials
    times_nosubsamp./=nTrials

    return times_subsamp, times_nosubsamp
end
times_ss, times_nss = time_subsamplingSize([12, 1], 18, nTrials=5)
scatter(collect(1:18).+log2(64), times_nss, yscale=:log10, label="No subsampling", legend=:topleft)
scatter!(collect(1:18).+log2(64), times_ss, yscale=:log10, label="With subsampling")
xlabel!("log2(Input size [bits])")
ylabel!("Time [s]")
savefig("time_subsampling_size.png")
mean(times_nss./times_ss)

function time_scatteringTransform_Q(Q::Array{Int64,2}, N::Int64)
    # define the scattering transform
    # first layer parameters and filter bank
    J=8
    σ0=0.1
    x=ones(2^N)

    times = zeros(size(Q,1))
    nFilters = zeros(size(Q,1))
    for i=1:size(Q,1)
        FilterBanks = get_FilterBanks(N, Q[i,:], J, σ0)
        Transform = ScatteringTransform1d(2, Q[i,:], [J,J])
        println(length(FilterBanks[1].Λ) * length(FilterBanks[2].Λ))
        nFilters[i] = length(FilterBanks[1].Λ)*length(FilterBanks[2].Λ)
        times[i] = @elapsed Scat.scatter1d(x, Transform, FilterBanks)
    end

    return nFilters, times
end

# vary the number of filters in the first layer
nQ=100
Q = ones(nQ, 2)
#Q[:, 1] .= Q[:, 1].*collect(1:100)
Q[:,1].= rand(1:12,nQ)
Q[:,2].=rand(1:4,nQ)
Q = convert.(Int, Q)
N=12
nFilters, times = time_scatteringTransform_Q(Q, N)
scatter(Q[:,1].*Q[:,2], times, legend=:none)
xlabel!("Q1 ⋅ Q2")
ylabel!("Time [s]")

scatter(nFilters, times, legend=:none)
xlabel!("Total Number of Filters")
ylabel!("Time [s]")
savefig("nFilters_times.png")

function time_subsamplingQ(Q::Array{Int64,2}, N::Int64)
    # define the scattering transform
    # first layer parameters and filter bank
    J=8
    σ0=0.1
    x=ones(2^N)

    times_ss = zeros(size(Q,1))
    times_nss = zeros(size(Q,1))
    nFilters = zeros(size(Q,1))
    for i=1:size(Q,1)
        FilterBanks = get_FilterBanks(N, Q[i,:], J, σ0)
        Transform = ScatteringTransform1d(2, Q[i,:], [J,J])
        println(length(FilterBanks[1].Λ) * length(FilterBanks[2].Λ))
        nFilters[i] = length(FilterBanks[1].Λ)*length(FilterBanks[2].Λ)
        times_ss[i] = @elapsed Scat.scatter1d(x, Transform, FilterBanks)
        times_nss[i] = @elapsed Scat.scatter1d(x, Transform, FilterBanks, subsample=false)
    end

    return nFilters, times_ss, times_nss
end
# vary the number of filters in the first layer
nQ=100
Q = ones(nQ, 2)
#Q[:, 1] .= Q[:, 1].*collect(1:100)
Q[:,1].= rand(1:12,nQ)
Q[:,2].=rand(1:4,nQ)
Q = convert.(Int, Q)
N=12
nFilters, times_ss, times_nss = time_subsamplingQ(Q, N)
scatter(nFilters, times_nss, label="No subsampling")
scatter!(nFilters, times_ss, legend=:topleft, label="With subsampling")
xlabel!("Q1 ⋅ Q2")
ylabel!("Time [s]")
savefig("time_subsampling_nFilters.png")

function time_scatteringTransform_J(J::Vector{Int64}, N::Int64)
    # define the scattering transform
    # first layer parameters and filter bank
    σ0=0.1
    Q=[12, 2]

    times = zeros(length(J))
    for i=1:length(J)
        println(J[i])
        x=ones(2^N)
        FilterBanks = get_FilterBanks(N, Q, J[i], σ0)
        Transform = ScatteringTransform1d(2, Q, [J[i], J[i]])
        times[i] = @elapsed Scat.scatter1d(x, Transform, FilterBanks)
    end

    return times
end
N=12
J = collect(4:20)
times = time_scatteringTransform_J(J, N)
scatter(J, sqrt.(times), legend=:none)
xlabel!("log2(maximum scale)")
ylabel!("√Time")
savefig("maxScale_time.png")

function time_subsamplingJ(J::Vector{Int64}, N::Int64; nTrials=1)
    # define the scattering transform
    # first layer parameters and filter bank
    σ0=0.1
    Q=[12, 2]

    times = zeros(length(J))
    for n=1:nTrials
        for i=1:length(J)
            println(J[i])
            x=ones(2^N)
            FilterBanks = get_FilterBanks(N, Q, J[i], σ0)
            Transform = ScatteringTransform1d(2, Q, [J[i], J[i]])
            times[i] += @elapsed Scat.scatter1d(x, Transform, FilterBanks)
        end
    end

    return times./nTrials
end
N=12
J = collect(4:16)
times = time_subsamplingJ(J, N, nTrials=15)
scatter(J, sqrt.(times), legend=:none)
xlabel!("log2(maximum scale)")
ylabel!("√Time")
savefig("maxScale_time.png")

function time_scatteringTransform_BreadthDepth(N)
    # define the scattering transform
    J=[8,8]
    σ0=0.1
    Q=[12,6]

    times_depth = zeros(N)
    times_breadth = zeros(N)
    for n=1:N
        println(n)
        x=ones(2^n)
        FilterBanks = get_FilterBanks(n, Q, J, σ0)
        FilterBankBlocks = get_FilterBanks(n, Q, J, σ0, typ=FilterBank1dBlock)
        Transform = ScatteringTransform1d(2, Q, J)
        times_depth[n] = @elapsed Scat.scatter1d(x, Transform, FilterBanks)
        times_breadth[n] = @elapsed Scat.scatter1d(x, Transform, FilterBankBlocks)
    end

    return times_depth, times_breadth
end
td, tb = time_scatteringTransform_BreadthDepth(14)
scatter(collect(1:14).+log2(64), td, yscale=:log10, label="Depth-first", legend=:topleft)
scatter!(collect(1:14).+log2(64), tb, yscale=:log10, label="Breadth-first")
xlabel!("log2(Input size) [bits]")
ylabel!("Time [s]")
savefig("depth_vs_breadth_time.png")


text:
### Extra utility methods

Host.hostmap(obj::API.Object, map::RT.RTMap{Name,MT,K,V,ME,F}) where {Name,MT,K,V,ME,F} =
    Host.hostmap(API.findmap(obj, String(Name)); K=K, V=V)
Host.hostmap(probe::AbstractProbe, map::RT.RTMap{Name,MT,K,V,ME,F}) where {Name,MT,K,V,ME,F} =
    Host.hostmap(API.findmap(only(objects(probe)), String(Name)); K=K, V=V)


text:
# RwaJur1

export 	RwaJur1Exposure,
		getExposureArray,
		RwaJur1Parameters,
		getCovarianceMatrix,
		getVaR,
		getPJur1,
		getRwaJur1,
		VERTICES_RWAJUR1

"""
Standardized vertices for mapping BRL fixed rate exposures (RWAJUR1)
"""
const VERTICES_RWAJUR1 = [21, 42, 63, 126, 252, 504, 756, 1008, 1260, 2520]

"""
    RwaJur1Exposure

Holds exposures in BRL fixed rate.

# Fields
* `date::Date`: date of evaluation
* `VXX::Float64`: exposure to vertex XX,
given by the sum of present value in BRL
of all cashflows that matures in XX business days (BD).

# Constructors
    RwaJur1Exposure(date::Date)

Create an instance with all exposures equal to zero.

	RwaJur1Exposure(date::Date, exposures::Array{Float64, 1})

Create an instance with exposures given by a vector of length 10.
Exposures are assigned in sequence from vertex 21 BD to vertex 2520 BD.

"""
type RwaJur1Exposure
	# Date of evaluation
	date::Date
	# Exposures in standard vertices: VXX, where XX is a maturity calculated in business days
	#V1::Float64
	V21::Float64
	V42::Float64
	V63::Float64
	V126::Float64
	V252::Float64
	V504::Float64
	V756::Float64
	V1008::Float64
	V1260::Float64
	V2520::Float64

	# Constructors
	RwaJur1Exposure(date::Date) = new(date, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

	function RwaJur1Exposure(date::Date, exposures::Array{Float64, 1})
		if size(exposures, 1) != 10
			error("exposures must have 10 values")
		end #if

		new(date,
		    exposures[1], exposures[2], exposures[3], exposures[4],
			exposures[5], exposures[6], exposures[7], exposures[8],
			exposures[9], exposures[10])
	end #function

end #type

"""
    getExposureArray(rExp::RwaJur1Exposure)

Return a vector with all exposures in sequence,
from vertex 21 BD to vertex 2520 BD.

"""
getExposureArray(rExp::RwaJur1Exposure) = [rExp.V21, rExp.V42, rExp.V63, rExp.V126, rExp.V252, rExp.V504, rExp.V756, rExp.V1008, rExp.V1260, rExp.V2520]

"""
	RWAJUR1Parameters

Holds BCB parameters to calculate fixed BRL stantardized capital (RWA JUR1)

"""
type RwaJur1Parameters
	date::Date
	M1::Float64
	rho::Float64
	rhoStress::Float64
	K::Float64
	KStress::Float64
	sigmaShortTerm::Float64
	sigmaMediumTerm::Float64
	sigmaLongTerm::Float64
	sigmaShortTermStress::Float64
	sigmaMediumTermStress::Float64
	sigmaLongTermStress::Float64
end #type

"""
	getCovarianceMatrix(par::RwaJur1Parameters, stress=false)

Return a variance-covariance matrix from RWA JUR1 parameters
supplied by BCB.

"""
function getCovarianceMatrix(par::RwaJur1Parameters, stress=false)

	if stress
		vecSigma = [ repmat([par.sigmaShortTermStress], 3);
					repmat([par.sigmaMediumTermStress], 3);
					repmat([par.sigmaLongTermStress], 4)]
	else
		vecSigma = [repmat([par.sigmaShortTerm], 3);
					repmat([par.sigmaMediumTerm], 3);
					repmat([par.sigmaLongTerm], 4)]
	end #if

	Rho = repmat([0.0], 10, 10)

	rho = stress? par.rhoStress : par.rho
	K = stress? par.KStress : par.K
	for i = 1:10
		for j = 1:10
			Rho[i, j] = rho + (1 - rho)^((max(VERTICES_RWAJUR1[i], VERTICES_RWAJUR1[j])/min(VERTICES_RWAJUR1[i], VERTICES_RWAJUR1[j]))^K)
		end #for j
	end #for i

	diagm(vecSigma) * Rho * diagm(vecSigma)

end #function

"""
	getVaR(rExp::RwaJur1Exposure, par::RwaJur1Parameters, stress=false)

Return portfolio Value-at-Risk (VaR) from its exposures and
from BCB RWA JUR1 parameters.
If `stress==true`, return stressed VaR (sVaR).

"""
function getVaR(rExp::RwaJur1Exposure, par::RwaJur1Parameters, stress=false)

	Sigma = getCovarianceMatrix(par, stress)
	X = getExposureArray(rExp) .* (VERTICES_RWAJUR1/252)

	2.33 * sqrt(X' * Sigma * X)[1] * sqrt(10)

end # function

"""
	getPJur1(rExp::RwaJur1Exposure, par::RwaJur1Parameters, VaRSMA60=0, sVaRSMA60=0)

Calculate PJUR1.

# Parameters
* `VaRSMA60=0`: VaR simple moving average, using last 60 business days
* `sVaRSMA60=0`: stressed VaR (sVaR) simple moving average, using last 60 business days

"""
function getPJur1(rExp::RwaJur1Exposure, par::RwaJur1Parameters, VaRSMA60=0, sVaRSMA60=0)

	VaR = getVaR(rExp, par, false)
	sVaR = getVaR(rExp, par, true)

	(max(par.M1 * VaRSMA60, VaR) + max(sVaRSMA60, sVaR))

end #function

"""
	getRwaJur1(rExp::RwaJur1Exposure, par::RwaJur1Parameters, VaRSMA60=0, sVaRSMA60=0)

Calculate risk-weighted assets for fixed BRL exposures (RWA JUR1).

# Parameters
* `VaRSMA60=0`:
VaR simple moving average, using last 60 business days.
See `getVaR()` for more info.
* `sVaRSMA60=0`:
stressed VaR (sVaR) simple moving average, using last 60 business days.
See `getVaR()` for more info.

"""
function getRwaJur1(rExp::RwaJur1Exposure, par::RwaJur1Parameters, VaRSMA60=0, sVaRSMA60=0)

	(1/getF(par.date)) * getPJur1(rExp, par, VaRSMA60, sVaRSMA60)

end #function


text:
function linearTrussStiffnessMatrix( E, A, nodalCoords  )

  diff = nodalCoords[ 2, : ] - nodalCoords[ 1, : ] ;

  length   = sqrt( diff' * diff ) ;
  c        = diff[1] / length ;
  s        = diff[3] / length ;

  Qloc2glo = [  c -s 0  0 ;
                s  c 0  0 ;
                0 0  c -s ;
                0 0  s  c ] ;

  Kloc     = E*A/length * [ 1  0 -1 0 ;
                            0  0  0 0 ;
                            -1 0  1 0 ;
                            0  0  0 0 ] ;
  Kglo     = Qloc2glo * Kloc * transpose(Qloc2glo) ;

  return Kglo
end
#


text:
######################################################################
#
#    Q/nZ interface
#
######################################################################

mutable struct QmodnZ <: GrpAb
  trivialmodulus::Bool
  n::fmpz
  d::fmpz

  function QmodnZ()
    z = new()
    z.n = fmpz(1)
    z.d = fmpz(1)
    z.trivialmodulus = true
    return z
  end

  function QmodnZ(n::fmpz)
    z = new()
    if isone(abs(n))
      return QmodnZ()
    else
      z = new()
      z.trivialmodulus = false
      z.n = abs(n)
      z.d = one(FlintZZ)
      return z
    end
  end

  function QmodnZ(n::fmpq)
    if isintegral(n)
      return QmodnZ(numerator(n))
    else
      z = new()
      z.trivialmodulus = false
      z.n = numerator(n)
      z.d = denominator(n)
      return z
    end
  end

end

function show(io::IO, G::QmodnZ)
  if G.trivialmodulus
    print(io, "Q/Z")
  else
    if isone(G.d)
      print(io, "Q/", G.n, "Z")
    else
      print(io, "Q/(", G.n, "/", G.d, ")Z")
    end
  end
end

modulus(T::QmodnZ) = T.n

struct QmodnZElem <: GrpAbElem
  elt::fmpq
  parent::QmodnZ

  function QmodnZElem(R::QmodnZ, a::fmpq)
    if R.trivialmodulus
      q = fmpq(mod(numerator(a), denominator(a)), denominator(a))
    else
      q = fmpq(mod(R.d * numerator(a), R.n * denominator(a)), R.d * denominator(a))
    end
    @assert isintegral(R.d * (q - a)) && iszero(mod(numerator(R.d * (q - a)), R.n))
    return new(q, R)
  end
end

function show(io::IO, a::QmodnZElem)
  G = parent(a)
  if G.trivialmodulus
    print(io, "$(a.elt) + Z")
  else
    if isone(G.d)
      print(io, "$(a.elt) + ", G.n, "Z")
    else
      print(io, "$(a.elt) + (", G.n, "/", G.d, ")Z")
    end
  end
end

function +(a::QmodnZElem, b::QmodnZElem)
  return QmodnZElem(a.parent, a.elt + b.elt)
end

function *(a::fmpz, b::QmodnZElem)
  return QmodnZElem(b.parent, a*b.elt)
end

function *(a::Integer, b::QmodnZElem)
  return QmodnZElem(b.parent, a*b.elt)
end

function divexact(a::QmodnZElem, b::fmpz)
  iszero(b) && throw(DivideError())
  return QmodnZElem(a.parent, a.elt // b)
end

function divexact(a::QmodnZElem, b::Integer)
  iszero(b) && throw(DivideError())
  return QmodnZElem(a.parent, a.elt // b)
end

function root_of_one(::Type{QmodnZ}, n::fmpz)
  return QmodnZElem(QmodnZ(), fmpq(1, n))
end

function inv(a::QmodnZElem)
  return QmodnZElem(a.parent, -(a.elt))
end

function parent(x::QmodnZElem)
  return x.parent
end

elem_type(::Type{QmodnZ}) = QmodnZElem

function order(a::QmodnZElem)
  if parent(a).trivialmodulus
    return denominator(a.elt)
  else
    throw(NotImplemented())
  end
end

(R::QmodnZ)(a::fmpz) = QmodnZElem(R, fmpq(a))
(R::QmodnZ)(a::Integer) = QmodnZElem(R, fmpq(a))
(R::QmodnZ)(a::fmpq) = QmodnZElem(R, a)
(R::QmodnZ)(a::Rational) = QmodnZElem(R, fmpq(a))

function Base.:(==)(a::QmodnZElem, b::QmodnZElem)
  if parent(a).trivialmodulus
    return isintegral(a.elt - b.elt)
  else
    z = a.elt - b.elt
    d = denominator(z)
    return isone(d) && iszero(mod(numerator(z), modulus(parent(a))))
  end
end

lift(a::QmodnZElem) = a.elt

################################################################################
#
#   Dual
#
################################################################################

mutable struct GrpAbFinGenToQmodnZ <: Map{GrpAbFinGen, QmodnZ,
                                            HeckeMap, GrpAbFinGenToQmodnZ}
  header::MapHeader{GrpAbFinGen, QmodnZ}

  function GrpAbFinGenToQmodnZ(G::GrpAbFinGen, QZ::QmodnZ, image)
    z = new()
    z.header = MapHeader(G, QZ, image)
    return z
  end
end

function kernel(mp::GrpAbFinGenToQmodnZ)
  G = domain(mp)
  QmodnZ = codomain(mp)
  n = exponent(G)
  C = abelian_group(fmpz[n])
  gensG = gens(G)
  imgs = Vector{GrpAbFinGenElem}()
  if QmodnZ.trivialmodulus
    for x in gensG
      imgx = lift(mp(x))
      o = divexact(n, denominator(imgx))
      push!(imgs, numerator(imgx)*o*C[1])
    end
  end
  mpk = hom(gensG, imgs)
  return kernel(mpk)
end

#TODO: technically, dual Z could be Q/Z ...
@doc Markdown.doc"""
    dual(G::GrpAbFinGen) -> GrpAbFinGen, Map

Computes the dual group, i.e. $hom(G, Q/Z)$ as an
abstract group. The map can be used to obtain actual homomorphisms.
It assumes that the group is finite.
"""
function dual(G::GrpAbFinGen)
  T, mT = torsion_subgroup(G)
  u = root_of_one(QmodnZ, exponent(T))
  return dual(G, u)
end

function dual(G::GrpAbFinGen, u::QmodnZElem)
  o = order(u)
  H = abelian_group(fmpz[o])
  QZ = parent(u)
  R, phi = hom(G, H)
  R::GrpAbFinGen
  ex = MapFromFunc(x -> x[1]*u, y -> H(fmpz[numerator(y.elt) * div(o, denominator(y.elt))]), H, parent(u))
  local mu
  let phi = phi, G = G, QZ = QZ, u = u
    function mu(r::GrpAbFinGenElem)
      f = phi(r)
      return GrpAbFinGenToQmodnZ(G, QZ, x -> f(x)[1]*u)
    end
  end

  local nu
  let ex = ex, phi = phi
    function nu(f::Map)
      g = GrpAbFinGenMap(f*inv(ex))
      return preimage(phi, g)
    end
  end
  return R, MapFromFunc(mu, nu, R, MapParent(R, parent(u), "homomorphisms"))
end


text:
# include("../../../src/plotSol.jl")
# include("../../../src/simulate.jl")
# include("../../../src/linearize.jl")
# include("../../../src/OrientationConversion.jl")
# include("../../../src/trim_kronLazy.jl");
using UdKaDyn, LinearAlgebra
# using Revise, UdKaDyn
# using BenchmarkTools

include("simplePendulumMinReal.jl")
##
# clearconsole()
j1 = Nothing
##
m = 1.0; l = 1.0; # Mass and length of bar
# Assuming bar revolves about Y axis
I1 = [1 0 0; 0 m * l^2 / 12 0; 0 0 1];

# Testing Dynamics with Revolute Joint
R1 = RigidBody(m, I1, 2)
RbI = InertialFrameAsRB()

# Suspended at an angle theta from vertical
theta = deg2rad(30)
x_temp = [l/2*sin(theta); 0.0; -l/2*cos(theta)];
x0R1 = ([x_temp;[1;zeros(3)];zeros(3);zeros(4)]);
initialiseRigidBody!(R1,x0R1)

# Axis about which bar is revolving
axisY = [0.0 1.0 0.0][:];

rj1 = [0.0 0.0 0.0][:]; # Joint Location in body frame of first body
rj2 = -R1.x[1:3];

j1 = Joint(RbI, R1, rj1, rj2, type="Revolute", axis=axisY)

# External Forces Definition
g = [0.0,0.0,-9.806]; # Gravity Vector.

j = [j1];

## Linearization
# x0, u0 = getXU_0(j);
# A,B = linearize(x0, u0, j, g);

## Simulation and checking
tEnd = (10.0);
tInt = 1e-2; # Only for fixed step solver, currently using adaptive step.
tSim, solFinal = simulate(tEnd, tInt, j, g=g)
solPend1  = solFinal[1];
println("Simulation complete., plotting results:");

##
# Check if joint location has moved
ωSol, jointLoc1 = checkRevJointIn(solPend1, rj2)
## Plots
close("all");
# errors
plotErrNorm(tSim,solPend1.β) # Drift in quaternion unit norm constraint
plotPos(tSim,jointLoc1) # Location of joint connecting inertial frame to pendulum. Should be zero.

# trajectory
plotQuat(tSim, solPend1.β)
plotPos(tSim, solPend1.r)
plotAngVel(tSim, ωSol)

## Compare with MinReal
rMin, vMin, ωMin = PendulumMinReal(m, l, theta, tInt, tEnd);
# # Plotting errors
plotPos(tSim,rMin-solPend1.r)
plotVel(tSim,vMin-solPend1.v)
plotAngVel(tSim,ωMin-ωSol)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
###################################
"""
    track = loadtrack("simlation/path/_output"::AbstractString)
    track = loadtrack("simlation/path/_output/fort.track"::AbstractString)

Function: fort.track reader
"""
function loadtrack(outdir::AbstractString)
    ## definition of filename
    fname = occursin("fort.track", basename(outdir)) ? outdir : joinpath(outdir, "fort.track")
    ## check
    isfile(fname) || error("File $fname is not found.")

    ## load
    datorg = readdlm(fname)
    datorg[datorg.>1e30] .= NaN

    ## return
    return VisClaw.Track(datorg[:,1], datorg[:,2], datorg[:,3], datorg[:,4])
end


text:
# TODO `solve = true`.
config = MOI.Test.TestConfig(solve = false)
_optimize!(mock) = MOIU.mock_optimize!(mock, zeros(MOI.get(mock, MOI.NumberOfVariables())))
for mock in mocks(_optimize!)
    Tests.dsos_options_pricing_test(mock, config)
    Tests.sdsos_options_pricing_test(mock, config)
    Tests.sos_options_pricing_test(mock, config)
end


text:
#=
State matrix `A` from the transmission line model [1]; we refer to that paper
for the interpretation of the physical parameters in this matrix. Here `A` is a
sparse, structured, square matrix of order 2η, where η ≥ 2 is a parameter
(optional, default: 2). For any choice of (positive) parameters, this matrix
is stable in the sense of Hurwitz (i.e. its eigenvalues have strictly negative
real-part). A parameter `scale` (optional, default: 1e-9) to rescale the physical
quantities is applied to `A`.

[1] M. Althoff, B. H. Krogh, and O. Stursberg. Analyzing Reachability of
Linear Dynamic Systems with Parametric Uncertainties. Modeling, Design,
and Simulation of Systems with Uncertainties. 2011.
=#
function transmission_line(; η=2, R₀=1.00, Rd₀=10.0, L₀=1e-10, C₀=1e-13 * 4.00, scale=1e-9)
    A₁₁ = zeros(η, η)
    A₁₂ = Bidiagonal(fill(-1/C₀, η), fill(1/C₀, η-1), :U)
    A₂₁ = Bidiagonal(fill(1/L₀, η), fill(-1/L₀, η-1), :L)
    A₂₂ = Diagonal(vcat(-Rd₀/L₀, fill(-R₀/L₀, η-1)))
    A  = [A₁₁ A₁₂; A₂₁ A₂₂]
    return A .* scale
end


text:
# This example demonstrates how to use custom h5geo data containers within JUDI.
# The idea of the example is to generate SEGY files using JUDI API
# and after that use generated SEGY files to prepare HDF5 containers with h5geo library.
# After the h5geo containers are created they are planned to be used to calculate FWI, RSTM, TWRI.
# The example creates 'ColadaSeismic/examles/tmp' folder and stores all the calculated data there.
# For now I tested that with JUDI 'time-tests' branch wich is not finished yet and
# it throws exception when JUDI.Options::save_data_to_disk is 'true'.
# Thus I distribute previously modeled SEGY files with 'master' branch

using ColadaSeismic, JUDI, SegyIO, LinearAlgebra, PyCall, PyPlot

# 'h5map_segy' maps SEGY to HDF5 h5geo::H5Seis object so that the same API is provided
# but no any data is read from SEGY
function h5map_segy(segy_files::Array{String})
  if length(segy_files) < 1
    @error "No SEGY files were provided"
    return
  end

  h5geo = pyimport("h5geopy._h5geo")
  if isnothing(h5geo)
    @error "Unable to import h5geo"
    return
  end

  container_name = splitext(segy_files[1])[1]
  container_name *= "_mapped.h5"
  container = h5geo.createSeisContainerByName(container_name, h5geo.CreationType.CREATE_UNDER_NEW_NAME)
  if isnothing(container)
    @error "Unable to create HDF5 container: $container_name"
    return
  end

  # parameters to create seis object
  p = h5geo.SeisParam()
  # Cant understand why in julia inheritance from H5BaseObject doesn't work
  # p.lengthUnits = "mm"
  # p.temporalUnits = "microsecond"
  p.domain = h5geo.Domain.TWT
  p.dataType = h5geo.SeisDataType.PRESTACK
  p.surveyType = h5geo.SurveyType.TWO_D
  p.trcChunk = 5
  p.segyFiles = segy_files
  p.mapSEGY = true

  @info "segy_files[1]: $(segy_files[1])"
  @info "segy_files[2]: $(segy_files[2])"

  seis_name = "mapped_segy"
  seis = container.createSeis(seis_name, p, h5geo.CreationType.CREATE_OR_OVERWRITE)
  if isnothing(seis)
    @error "Unable to create seis"
    return
  end

  seis.setLengthUnits("mm")
  seis.setTemporalUnits("microsecond")
  return seis
end

# 'h5read_segy' read SEGY files to HDF5 h5geo::H5Seis object
function h5read_segy(segy_files::Array{String})
  if length(segy_files) < 1
    @error "No SEGY files were provided"
    return
  end

  h5geo = pyimport("h5geopy._h5geo")
  if isnothing(h5geo)
    @error "Unable to import h5geo"
    return
  end

  container_name = splitext(segy_files[1])[1]
  container_name *= "_read.h5"
  container = h5geo.createSeisContainerByName(container_name, h5geo.CreationType.CREATE_UNDER_NEW_NAME)
  if isnothing(container)
    @error "Unable to create HDF5 container: $container_name"
    return
  end

  # parameters to create seis object
  p = h5geo.SeisParam()
  # Cant understand why in julia inheritance from H5BaseObject doesn't work
  # p.lengthUnits = "mm"
  # p.temporalUnits = "microsecond"
  p.domain = h5geo.Domain.TWT
  p.dataType = h5geo.SeisDataType.PRESTACK
  p.surveyType = h5geo.SurveyType.TWO_D
  # explicitely set 'nSamp' so that chunking along 1-st dim is equal to 'nSamp'
  p.nSamp = h5geo.getSEGYNSamp(segy_files[1], h5geo.getSEGYEndian(segy_files[1]))
  # 'nTrc' must be bigger that 'trcChunk'
  p.nTrc = 10
  p.trcChunk = 5

  @info "segy_files[1]: $(segy_files[1])"
  @info "segy_files[2]: $(segy_files[2])"
  @info "p.nSamp: $(p.nSamp)"

  seis_name = "read_segy"
  seis = container.createSeis(seis_name, p, h5geo.CreationType.CREATE_OR_OVERWRITE)
  if isnothing(seis)
    @error "Unable to create seis"
    return
  end

  seis.readSEGYTextHeader(segy_files[1])
  seis.readSEGYBinHeader(segy_files[1])
  seis.readSEGYTraces(segy_files)
  seis.setLengthUnits("mm")
  seis.setTemporalUnits("microsecond")
  return seis
end

# Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
v0 = ones(Float32,n) .+ 0.5f0
v[:,Int(round(end/2)):end] .= 3.5f0
rho = (v0 .+ .5f0) ./ 2

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 2	# number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
timeR = 1000f0   # receiver recording time [ms]
dtR = 2f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(range(400f0, stop=800f0, length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(200f0, stop=200f0, length=nsrc))

# source sampling and number of time steps
timeS = 1000f0  # ms
dtS = 2f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

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

segy_path = @__DIR__
segy_path *= "/tmp/"
segy_name = "shot"

# rm(segy_path, force=true, recursive=true)
# mkpath(segy_path)

# Write shots as segy files to disk (gives error in 'time-tests' branch)
# opt = Options(optimal_checkpointing=false, isic=false, subsampling_factor=2, dt_comp=1.0,
#               save_data_to_disk=true, file_path=segy_path, file_name=segy_name)
opt = Options(optimal_checkpointing=false, isic=false, subsampling_factor=2, dt_comp=1.0)

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*adjoint(Ps), q)

# Nonlinear modeling
dobs = Pr*F*adjoint(Ps)*q

# modeled SEGY files
files = readdir(segy_path, join=true)
segy_files = Array{String}(undef, 0)
for file in files
  if splitext(file)[2] == ".segy"
    push!(segy_files, file)
  end
end

# Ether map SEGY (only ieee32 supported) or read (any 4-byte format) it

# seis = h5map_segy(segy_files)
seis = h5read_segy(segy_files)
if isnothing(seis)
  @error "Unable to create seis"
  return
end

# To retrieve sorted data we need to prepare sorting from primary keys (PKey)
# 'h5geo' uses trace header names that differs from the ones that uses SegyIO
if !seis.hasPKeySort("SRCX")
  seis.addPKeySort("SRCX")
end

# create H5SeisCon and judiVector from it
con = H5SeisCon(seis=seis, pkey="SRCX")
h5dobs = judiVector(con, xkey="GRPX", ykey="GRPY", zkey="RGE")

# as long as SEGY stores coordinates in Int32 and h5geo stores it in Float64
# we need to set previously calculated 'recGeometry' to prevent round off mismatch error (needed only for this example)
h5dobs.geometry = recGeometry

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

# Continue colaculations  (FWI, RSTM etc) with 'judiVector{H5SeisCon}'

# Adjoint with h5dobs
pdeFull = Ps*adjoint(F)*adjoint(Pr)
qad = pdeFull*h5dobs

# Linearized modeling
dD = J*dm
# Adjoint jacobian
rtm = adjoint(J)*dD

# evaluate FWI objective function with h5dobs
f, g = fwi_objective(model0, q, h5dobs; options=opt)

imshow(g'); title("FWI (g')")
savefig(segy_path * "/g.png")

# TWRI with h5dobs
f, gm, gy = twri_objective(model0, q, h5dobs, nothing; options=opt, optionswri=TWRIOptions(params=:all))
f, gm = twri_objective(model0, q, h5dobs, nothing; options=Options(frequencies=[[.009, .011], [.008, .012]]),
                                                 optionswri=TWRIOptions(params=:m))

imshow(gm'); title("TWRI (gm')")
savefig(segy_path * "/gm.png")

# evaluate LSRTM objective function with h5dobs
fj, gj = lsrtm_objective(model0, q, dD, dm; options=opt)
fjn, gjn = lsrtm_objective(model0, q, h5dobs, dm; nlind=true, options=opt)

imshow(gj'); title("LSRTM (gj')")
savefig(segy_path * "/gj.png")
imshow(gjn'); title("LSRTM (gjn')")
savefig(segy_path * "/gjn.png")

# By extension, lsrtm_objective is the same as fwi_objecive when `dm` is zero
# And with computing of the residual. Small noise can be seen in the difference
# due to floating point roundoff errors with openMP, but running with
# OMP_NUM_THREAS=1 (no parllelism) produces the exact (difference == 0) same result
# gjn2 == g
fjn2, gjn2 = lsrtm_objective(model0, q, h5dobs, 0f0.*dm; nlind=true, options=opt)

imshow(gjn2'); title("LSRTM (gjn2')")
savefig(segy_path * "/gjn2.png")


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

mutable struct OnErrorDeployment <: SwaggerModel
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    deploymentName::Any # spec type: Union{ Nothing, String } # spec name: deploymentName

    function OnErrorDeployment(;type=nothing, deploymentName=nothing)
        o = new()
        validate_property(OnErrorDeployment, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(OnErrorDeployment, Symbol("deploymentName"), deploymentName)
        setfield!(o, Symbol("deploymentName"), deploymentName)
        o
    end
end # type OnErrorDeployment

const _property_map_OnErrorDeployment = Dict{Symbol,Symbol}(Symbol("type")=>Symbol("type"), Symbol("deploymentName")=>Symbol("deploymentName"))
const _property_types_OnErrorDeployment = Dict{Symbol,String}(Symbol("type")=>"String", Symbol("deploymentName")=>"String")
Base.propertynames(::Type{ OnErrorDeployment }) = collect(keys(_property_map_OnErrorDeployment))
Swagger.property_type(::Type{ OnErrorDeployment }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OnErrorDeployment[name]))}
Swagger.field_name(::Type{ OnErrorDeployment }, property_name::Symbol) =  _property_map_OnErrorDeployment[property_name]

const _allowed_OnErrorDeployment_type = ["LastSuccessful", "SpecificDeployment"]

function check_required(o::OnErrorDeployment)
    true
end

function validate_property(::Type{ OnErrorDeployment }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "OnErrorDeployment", :enum, val, _allowed_OnErrorDeployment_type)
    end
end


text:
import IterTools: groupby, imap
using Zygote: @nograd

struct BSplineBasis <: Basis1D
    knots::Vector{Float64}
    order::Int
    deriv::Int

    function BSplineBasis(knots, order, deriv=0; extend=true)
        if deriv >= order throw(ArgumentError("Differentiation order not supported")) end

        for (kn, kp) in zip(knots[2:end], knots[1:end-1])
            if kn < kp throw(ArgumentError("Knot vector must be nondecreasing")) end
        end

        if extend
            knots = [fill(knots[1], order-1); knots; fill(knots[end], order-1)]
        else
            for d in 1:order-1
                if knots[1+d] != knots[1] || knots[end-d] != knots[end]
                    throw(ArgumentError("Expected $order repeated knots on either end"))
                end
            end
        end

        new(knots, order, deriv)
    end

    BSplineBasis(lft::Real, rgt::Real, elements::Int, order::Int) =
        BSplineBasis(range(lft, stop=rgt, length=elements+1), order)
end

const BSpline = BasisFunction1D{BSplineBasis}

Base.length(b::BSplineBasis) = length(b.knots) - b.order
Base.size(b::BSplineBasis) = (length(b),)
Base.getindex(b::BSplineBasis, i) = BSpline(b, i, b.deriv)
nderivs(b::BSplineBasis) = b.order - 1

domain(b::BSplineBasis) = Interval(b.knots[1], b.knots[end])
domain(b::BSpline) = Interval(b.basis.knots[b.index], b.basis.knots[b.index+b.basis.order])

degree(b::BSplineBasis) = b.order - b.deriv - 1
degree(b::BSpline) = b.basis.order - b.deriv - 1

deriv(b::BSplineBasis) = BSplineBasis(b.knots, b.order, b.deriv+1; extend=false)
deriv(b::BSplineBasis, order) = BSplineBasis(b.knots, b.order, b.deriv+order; extend=false)

function supported(b::BSplineBasis, pt::T) where {T<:Real}
    kidx = b.order - 1 + searchsorted(b.knots[b.order:end], pt).stop
    stop = b.knots[kidx] == b.knots[end] ? kidx - b.order : kidx
    stop - b.order + 1 : stop
end

function supported(b::BSplineBasis, pts::Vector{T}) where {T<:Real}
    (min, max) = extrema(pts)
    @assert(min in domain(b) && max in domain(b), "Domain error")

    idxs = zeros(Int, length(pts))

    if !issorted(pts)
        for (i, pt) in enumerate(pts)
            idxs[i] = supported(b, pt).stop
        end
    else
        kidx = b.order
        for (i, pt) in enumerate(pts)
            kidx += searchsorted(b.knots[kidx:end], pt).stop - 1
            idxs[i] = b.knots[kidx] == b.knots[end] ? kidx - b.order : kidx
        end
    end

    imap(groupby(i -> i[2], enumerate(idxs))) do i
        (pts[i[1][1]:i[end][1]], i[1][2] - b.order + 1 : i[1][2])
    end
end

@nograd function evaluate_raw(b::BSplineBasis, pts::Vector{T}, deriv::Int, rng::UnitRange{Int}) where {T<:Real}
    # Basis values of order 1 (piecewise constants)
    bvals = zeros(Float64, (b.order, length(pts)))
    bvals[end,:] .= 1.0

    p = b.order
    bi = rng.start + p

    # Order increment
    for k in 0:p-deriv-2
        bvals[p-k:end, :] ./= b.knots[bi:bi+k] - b.knots[bi-k-1:bi-1]

        for (i, kp, kn) in zip(p-k-1:p-1, b.knots[bi-k-2:bi-2], b.knots[bi:bi+k])
            bvals[i,:] .*= (pts .- kp)
            bvals[i,:] .+= bvals[i+1,:] .* (kn .- pts)
        end
        bvals[end,:] .*= (pts .- b.knots[bi-1])
    end

    # Differentiation
    for k = p-deriv-1:p-2
        bvals[p-k:end, :] ./= b.knots[bi:bi+k] - b.knots[bi-k-1:bi-1]

        bvals[1:end-1,:] = -diff(bvals, dims=1)
        bvals .*= k + 1
    end

    bvals
end


text:
julia> z1 = 1.5 + 3im
julia> z2 = 1.5 + 1.5im
julia> z1 + z2
3.0 + 4.5im
julia> z1 - z2
0.0 + 1.5im
julia> z1 * z2
-2.25 + 6.75im
julia> z1 / z2
1.5 + 0.5im
julia> - z1
-1.5 - 3.0im
julia> conj(z1), z1'   # two ways to conjugate
(1.5 - 3.0im,1.5 - 3.0im)
julia> abs(z1)
3.3541019662496847
julia> z1^z2
-1.102482955327779 - 0.38306415117199305im
julia> real(z1)
1.5
julia> imag(z1)
3.0

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Printf
import MathOptInterface

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

# supported GAMS command line options
const MOI_SUPPORTED_CLOPTIONS = (
   "reslim",
   "iterlim",
   "holdfixed",
   "nodlim",
   "optca",
   "optcr",
   "solver",
   "threads",
   "trace",
   "traceopt",
   "license",
   "logoption",
   "lp",
   "mip",
   "rmip",
   "nlp",
   "dnlp",
   "cns",
   "minlp",
   "rminlp",
   "qcp",
   "miqcp",
   "rmiqcp"
)

# supported model types
const MOI_SUPPORTED_MODEL_TYPES = (
   GAMS.MODEL_TYPE_LP,
   GAMS.MODEL_TYPE_MIP,
   GAMS.MODEL_TYPE_RMIP,
   GAMS.MODEL_TYPE_NLP,
   GAMS.MODEL_TYPE_MINLP,
   GAMS.MODEL_TYPE_RMINLP,
   GAMS.MODEL_TYPE_QCP,
   GAMS.MODEL_TYPE_MIQCP,
   GAMS.MODEL_TYPE_RMIQCP,
   GAMS.MODEL_TYPE_MCP,
   GAMS.MODEL_TYPE_MPEC,
)

mutable struct VariableInfo
   name::String
   type::GAMSVarType
   lower_bound::Union{Float64, Nothing}
   lower_bound_dual_start::Union{Float64, Nothing}
   upper_bound::Union{Float64, Nothing}
   upper_bound_dual_start::Union{Float64, Nothing}
   start::Union{Float64, Nothing}
end

VariableInfo() = VariableInfo("", VARTYPE_FREE, nothing, nothing, nothing,
   nothing, nothing)

mutable struct ConstraintInfo{F, S}
   name::String
   func::F
   set::S
   dual_start::Union{Nothing, Float64}
end

ConstraintInfo(func, set) = ConstraintInfo("", func, set, nothing)

mutable struct Optimizer <: MOI.AbstractOptimizer
   gamswork::Union{GAMSWorkspace, Nothing}

   # problem attributes
   model_type::GAMSModelType
   n_binary::Int
   n_integer::Int
   n_semicont::Int
   n_semiint::Int
   m::Int
   m_lin::Int
   m_quad::Int
   m_nonlin::Int
   sense::MOI.OptimizationSense
   objective::Union{MOI.SingleVariable, MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}, Nothing}
   objvar::Bool
   variable_info::Vector{VariableInfo}
   linear_le_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}}
   linear_ge_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}}
   linear_eq_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}}
   quadratic_le_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}}
   quadratic_ge_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}}
   quadratic_eq_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}}
   sos1_constraints::Vector{ConstraintInfo{MOI.VectorOfVariables, MOI.SOS1{Float64}}}
   sos2_constraints::Vector{ConstraintInfo{MOI.VectorOfVariables, MOI.SOS2{Float64}}}
   nlp_data::Union{MOI.NLPBlockData, Nothing}
   complementarity_constraints::Vector{ConstraintInfo{MOI.VectorAffineFunction{Float64}, MOI.Complements}}

   # parameters
   sysdir::Union{String, Nothing}
   workdir::Union{String, Nothing}
   user_model_type::GAMSModelType
   gams_options::Dict{String, Any}
   solver_options::Dict{String, Any}

   # solution attributes
   solve_time::Float64
   solve_status::GAMSSolveStatus
   model_status::GAMSModelStatus
   sol::Union{GAMSSolution, Nothing}
   obj::Float64
   obj_est::Float64
end

function Optimizer(workspace::Union{Nothing, GAMSWorkspace} = nothing)
   gams_options = Dict{String, Any}()
   gams_options["threads"] = 1
   return Optimizer(workspace, MODEL_TYPE_UNDEFINED, 0, 0, 0, 0, 0, 0, 0, 0,
                    MOI.FEASIBILITY_SENSE, nothing, true, [], [], [], [], [], [],
                    [], [], [], nothing, [], nothing, nothing, MODEL_TYPE_UNDEFINED,
                    gams_options, Dict{String, Any}(), NaN, SOLVE_STATUS_UNDEFINED,
                    MODEL_STATUS_UNDEFINED, nothing, NaN, NaN)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "GAMS"

MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.RawParameter) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{MOI.VariableIndex}) = true
MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{MOI.VariableIndex}) = true
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{MOI.ConstraintIndex}) = true
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = true
MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
MOI.supports(::Optimizer, ::MOI.NLPBlockDual) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.ZeroOne}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.Integer}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.Semicontinuous{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.Semiinteger{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{MOI.SOS1{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{MOI.SOS2{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorAffineFunction{Float64}}, ::Type{MOI.Complements}) = true
MOI.supports(::Optimizer, ::MOI.ConstraintDualStart, ::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}, ::Union{Real, Nothing}) = true
MOI.supports(::Optimizer, ::MOI.ConstraintDualStart, ::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}, ::Union{Real, Nothing}) = true
MOI.supports(::Optimizer, ::MOI.ConstraintDualStart, ::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}, ::Union{Real, Nothing}) = true

MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = !copy_names

function MOI.copy_to(
   model::Optimizer,
   src::MOI.ModelLike;
   copy_names = false
)
   return MOIU.default_copy_to(model, src, copy_names)
end

function MOI.empty!(
   model::Optimizer
)
   model.model_type = MODEL_TYPE_UNDEFINED
   model.n_binary = 0
   model.n_integer = 0
   model.n_semicont = 0
   model.n_semiint = 0
   model.m = 0
   model.m_lin = 0
   model.m_quad = 0
   model.m_nonlin = 0
   model.sense = MOI.FEASIBILITY_SENSE
   model.objective = nothing
   empty!(model.variable_info)
   empty!(model.linear_le_constraints)
   empty!(model.linear_ge_constraints)
   empty!(model.linear_eq_constraints)
   empty!(model.quadratic_le_constraints)
   empty!(model.quadratic_ge_constraints)
   empty!(model.quadratic_eq_constraints)
   empty!(model.sos1_constraints)
   empty!(model.sos2_constraints)
   model.nlp_data = nothing
   empty!(model.complementarity_constraints)
   model.solve_time = NaN
   model.solve_status = SOLVE_STATUS_UNDEFINED
   model.model_status = MODEL_STATUS_UNDEFINED
   model.sol = nothing
   model.obj = NaN
   model.obj_est = NaN
end

function MOI.is_empty(
   model::Optimizer
)
   return model.model_type == MODEL_TYPE_UNDEFINED
      model.n_binary == 0 &&
      model.n_integer == 0 &&
      model.n_semicont == 0 &&
      model.n_semiint == 0 &&
      model.m == 0 &&
      model.m_lin == 0 &&
      model.m_quad == 0 &&
      model.m_nonlin == 0 &&
      model.sense == MOI.FEASIBILITY_SENSE &&
      isnothing(model.objective) &&
      isempty(model.variable_info) &&
      isempty(model.linear_le_constraints) &&
      isempty(model.linear_ge_constraints) &&
      isempty(model.linear_eq_constraints) &&
      isempty(model.quadratic_le_constraints) &&
      isempty(model.quadratic_ge_constraints) &&
      isempty(model.quadratic_eq_constraints) &&
      isempty(model.sos1_constraints) &&
      isempty(model.sos2_constraints) &&
      isnothing(model.nlp_data) &&
      isempty(model.complementarity_constraints)
end

has_start(var::VariableInfo) = ! isnothing(var.start)
has_upper_bound(var::VariableInfo) = ! isnothing(var.upper_bound)
has_lower_bound(var::VariableInfo) = ! isnothing(var.lower_bound)
is_fixed(var::VariableInfo) = has_lower_bound(var) && has_upper_bound(var) && var.lower_bound == var.upper_bound

has_upper_bound(model::Optimizer, vi::MOI.VariableIndex) = has_upper_bound(model.variable_info[vi.value])
has_lower_bound(model::Optimizer, vi::MOI.VariableIndex) = has_lower_bound(model.variable_info[vi.value])
is_fixed(model::Optimizer, vi::MOI.VariableIndex) = is_fixed(model.variable_info[vi.value])

offset_linear_le(model::Optimizer) = 0
offset_linear_ge(model::Optimizer) = length(model.linear_le_constraints)
offset_linear_eq(model::Optimizer) = offset_linear_ge(model) + length(model.linear_ge_constraints)
offset_quadratic_le(model::Optimizer) = offset_linear_eq(model) + length(model.linear_eq_constraints)
offset_quadratic_ge(model::Optimizer) = offset_quadratic_le(model) + length(model.quadratic_le_constraints)
offset_quadratic_eq(model::Optimizer) = offset_quadratic_ge(model) + length(model.quadratic_ge_constraints)
offset_nonlin(model::Optimizer) = offset_quadratic_eq(model) + length(model.quadratic_eq_constraints)
offset_complementarity(model::Optimizer) = 0

function MOI.set(
   model::Optimizer,
   ::MOI.NLPBlock,
   nlp_data::MOI.NLPBlockData
)
   model.nlp_data = nlp_data
   MOI.initialize(model.nlp_data.evaluator, [:ExprGraph])
   return
end

include("options.jl")
include("variables.jl")
include("objective.jl")
include("constraints.jl")
include("translate.jl")
include("solve.jl")


text:
"""
    y = leftof(p, a, b)

Return true if the point `p` is to the left of the line `a``b`.
For a n×2 list of points `p`, return a vector of results.
"""
function leftof(p, a, b)
    # Translate to the origin and do the check
    p = p - ones(size(p, 1)) * a
    b = b - a
    return b[1]*p[:,2] .> b[2]*p[:,1]
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
import PencilFFTs

import Oceananigans.Solvers: poisson_eigenvalues, solve_poisson_equation!

struct DistributedFFTBasedPoissonSolver{P, F, L, λ, S}
              plan :: P
         full_grid :: F
           my_grid :: L
       eigenvalues :: λ
           storage :: S
end

function DistributedFFTBasedPoissonSolver(arch, full_grid, local_grid)
    topo = (TX, TY, TZ) = topology(full_grid)

    λx = poisson_eigenvalues(full_grid.Nx, full_grid.Lx, 1, TX())
    λy = poisson_eigenvalues(full_grid.Ny, full_grid.Ly, 2, TY())
    λz = poisson_eigenvalues(full_grid.Nz, full_grid.Lz, 3, TZ())

    I, J, K = arch.local_index
    λx = λx[(J-1)*local_grid.Ny+1:J*local_grid.Ny, :, :]

    eigenvalues = (; λx, λy, λz)

    transform = PencilFFTs.Transforms.FFT!()
    proc_dims = (arch.ranks[2], arch.ranks[3])
    plan = PencilFFTs.PencilFFTPlan(size(full_grid), transform, proc_dims, MPI.COMM_WORLD)
    storage = PencilFFTs.allocate_input(plan)

    return DistributedFFTBasedPoissonSolver(plan, full_grid, local_grid, eigenvalues, storage)
end

function solve_poisson_equation!(solver::DistributedFFTBasedPoissonSolver)
    λx, λy, λz = solver.eigenvalues

    # Apply forward transforms.
    solver.plan * solver.storage

    # Solve the discrete Poisson equation.
    RHS = ϕ = solver.storage[2]
    @. ϕ = - RHS / (λx + λy + λz)

    # Setting DC component of the solution (the mean) to be zero. This is also
    # necessary because the source term to the Poisson equation has zero mean
    # and so the DC component comes out to be ∞.
    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        ϕ[1, 1, 1] = 0
    end

    # Apply backward transforms.
    solver.plan \ solver.storage

    return nothing
end


text:
using CrypticCrosswords
const CC = CrypticCrosswords

@testset "Insertions" begin
    out = String[]
    buffer = Vector{Char}()
    CC.insertions!(out, buffer, "a", "b")
    @test isempty(out)

    empty!(out)
    CC.insertions!(out, buffer, "a", "bc")
    @test out == ["bac"]

    empty!(out)
    CC.insertions!(out, buffer, "", "bc")
    @test isempty(out)
end


text:
# sorts generators and corresponding coefficients from basis
# by their leading monomial in increasing ordering
#
# Used only once to sort input generators
function sort_gens_by_lead_increasing!(basis::Basis, ht::MonomialHashtable)
    gens = basis.gens
    exps = ht.exponents

    inds = collect(1:basis.ntotal)

    if ht.ord == :degrevlex
        cmp  = (x, y) -> exponent_isless_drl(
                                @inbounds(exps[gens[x][1]]),
                                @inbounds(exps[gens[y][1]]))
    elseif ht.ord == :lex
        cmp  = (x, y) -> exponent_isless_lex(
                                @inbounds(exps[gens[x][1]]),
                                @inbounds(exps[gens[y][1]]))
    else
        cmp  = (x, y) -> exponent_isless_dl(
                                @inbounds(exps[gens[x][1]]),
                                @inbounds(exps[gens[y][1]]))
    end
    sort!(inds, lt=cmp)

    # use array assignment insted of elemewise assignment
    # Reason: slice array assignment uses fast setindex! based on unsafe copy
    basis.gens[1:basis.ntotal]   = basis.gens[inds]
    basis.coeffs[1:basis.ntotal] = basis.coeffs[inds]
end

# sorts generators and corresponding coefficients from basis
# by their leading monomial in increasing ordering
#
# Used only once to sort output generators
function sort_gens_by_lead_increasing_in_reduce!(basis::Basis, ht::MonomialHashtable)
    gens = basis.gens
    exps = ht.exponents
    nnr  = basis.nonred

    inds = collect(1:basis.nlead)

    if ht.ord == :degrevlex
        cmp  = (x, y) -> exponent_isless_drl(
                                @inbounds(exps[gens[nnr[x]][1]]),
                                @inbounds(exps[gens[nnr[y]][1]]))
    elseif ht.ord == :lex
        cmp  = (x, y) -> exponent_isless_lex(
                                @inbounds(exps[gens[nnr[x]][1]]),
                                @inbounds(exps[gens[nnr[y]][1]]))
    else
        cmp  = (x, y) -> exponent_isless_dl(
                                @inbounds(exps[gens[nnr[x]][1]]),
                                @inbounds(exps[gens[nnr[y]][1]]))
    end
    sort!(inds, lt=cmp)

    # use array assignment insted of elemewise assignment
    # Reason: slice array assignment uses fast setindex! based on unsafe copy
    basis.gens[nnr[1:basis.nlead]]   = basis.gens[nnr[inds]]
    basis.coeffs[nnr[1:basis.nlead]] = basis.coeffs[nnr[inds]]
end

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

# Sorts pairs from pairset in range [from, from+size]
# by lcm total degrees in increasing order
#
# Used in update once per one f4 iteration to sort pairs in pairset
# Also used in normal selection strategy also once per iteration
function sort_pairset_by_degree!(ps::Pairset, from::Int, sz::Int)
    part = ps.pairs[from:from + sz]

    sort!(part, by=p -> p.deg)

    ps.pairs[from:from + sz] = part
end

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

# sorts first npairs pairs from pairset by increasing of
# lcm exponent wrt the given monomial ordering
#
# Used only in normal selection strategy once per f4 iteration
function sort_pairset_by_lcm!(
        pairset::Pairset, npairs::Int, ht::MonomialHashtable)

    part = pairset.pairs[1:npairs]
    exps = ht.exponents

    if ht.ord == :degrevlex
        sort!(part, lt=(x, y) -> exponent_isless_drl(
                                        @inbounds(exps[x.lcm]),
                                        @inbounds(exps[y.lcm])))
    elseif ht.ord == :lex
        sort!(part, lt=(x, y) -> exponent_isless_lex(
                                        @inbounds(exps[x.lcm]),
                                        @inbounds(exps[y.lcm])))
    else
        sort!(part, lt=(x, y) -> exponent_isless_dl(
                                        @inbounds(exps[x.lcm]),
                                        @inbounds(exps[y.lcm])))
    end

    pairset.pairs[1:npairs] = part
end

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

# sorts generators selected in normal strategy function
# by their ordering in the current basis (identity sort)
function sort_generators_by_position!(gens::Vector{Int}, load::Int)
    # @info "Sorting first $load gens"
    part = gens[1:load]

    sort!(part)

    gens[1:load] = part
end

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

#=
    :degrevlex exponent vector comparison
=#
function exponent_isless_drl(ea::Vector{UInt16}, eb::Vector{UInt16})

    if @inbounds ea[end] < eb[end]
        return true
    elseif @inbounds ea[end] != eb[end]
        return false
    end

    i = length(ea) - 1
    @inbounds while i > 1 && ea[i] == eb[i]
        i -= 1
    end

    return ea[i] < eb[i] ? false : true
end

#=
    :deglex exponent vector comparison
=#
function exponent_isless_dl(ea::Vector{UInt16}, eb::Vector{UInt16})

    if @inbounds ea[end] < eb[end]
        return true
    elseif @inbounds ea[end] != eb[end]
        return false
    end

    i = 1
    @inbounds while i < length(ea) && ea[i] == eb[i]
        i += 1
    end
    return ea[i] < eb[i] ? true : false
end

#=
    :lex exponent vector comparison
=#
function exponent_isless_lex(ea::Vector{UInt16}, eb::Vector{UInt16})
    i = 1
    @inbounds while i < length(ea) - 1 && ea[i] == eb[i]
        i += 1
    end
    return ea[i] < eb[i] ? true : false
end

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

function matrix_row_decreasing_cmp(a, b)
    #= a, b - rows as arrays of exponent hashes =#

    @inbounds va = a[1]
    @inbounds vb = b[1]

    if va > vb
        return false
    end
    if va < vb
        return true
    end

    # same column index => compare density of rows

    va = length(a)
    vb = length(b)
    if va > vb
        return false
    end
    if va < vb
        return true
    end

    # hmm, equal rows?
    return false
end

function matrix_row_increasing_cmp(a, b)
    #= a, b - rows as arrays of exponent hashes =#

    @inbounds va = a[1]
    @inbounds vb = b[1]

    if va > vb
        return false
    end
    if va < vb
        return true
    end

    # same column index => compare density of rows

    va = length(a)
    vb = length(b)
    if va > vb
        return false
    end
    if va < vb
        return true
    end

    # hmm, equal rows?
    return false
end

# by column index and density
function sort_matrix_rows_decreasing!(matrix)
    #= smaller means  pivot being more left  =#
    #= and density being smaller             =#

    # @info "before up sort"
    #println(matrix.uprows)
    #println(matrix.up2coef)

    inds = collect(1:matrix.nup)
    cmp  = (x, y) ->  matrix_row_decreasing_cmp(
                                    @inbounds(matrix.uprows[x]),
                                    @inbounds(matrix.uprows[y]))
    sort!(inds, lt=cmp)

    matrix.uprows[1:matrix.nup]  = matrix.uprows[inds]
    matrix.up2coef[1:matrix.nup] = matrix.up2coef[inds]

    # @info "after up sort"
    #println(matrix.uprows)
    #println(matrix.up2coef)

    matrix
end

# by column index and density
function sort_matrix_rows_increasing!(matrix)
    #= smaller means  pivot being more right =#
    #= and density being larger =#

    # @info "before low sort"
    #println(matrix.lowrows)
    #println(matrix.low2coef)

    inds = collect(1:matrix.nlow)
    cmp  = (x, y) ->  matrix_row_increasing_cmp(
                                        @inbounds(matrix.lowrows[x]),
                                        @inbounds(matrix.lowrows[y]))
    sort!(inds, lt=cmp)

    matrix.lowrows[1:matrix.nlow]  = matrix.lowrows[inds]
    matrix.low2coef[1:matrix.nlow] = matrix.low2coef[inds]

    # @info "after low sort"
    #println(matrix.lowrows)
    #println(matrix.low2coef)

    matrix
end

function sort_columns_by_hash!(col2hash, symbol_ht)
    hd = symbol_ht.hashdata
    es = symbol_ht.exponents
    len = symbol_ht.explen

    function cmp(a, b, exponent_isless)
        @inbounds ha = hd[a]
        @inbounds hb = hd[b]
        if ha.idx != hb.idx
            return ha.idx > hb.idx
        end
        @inbounds ea = es[a]
        @inbounds eb = es[b]
        !exponent_isless(ea, eb)
    end

    if symbol_ht.ord == :degrevlex
        ordcmp = (x, y) -> cmp(x, y, exponent_isless_drl)
    elseif symbol_ht.ord == :lex
        ordcmp = (x, y) -> cmp(x, y, exponent_isless_lex)
    else
        ordcmp = (x, y) -> cmp(x, y, exponent_isless_dl)
    end

    sort!(col2hash, lt = ordcmp)
end


text:
"""
    mutable struct Graph{T<:Integer} <: AGraph
        ne::Int
        fadjlist::Vector{Vector{T}}
    end

A simple graph type based on an adjacency list.

    Graph{T}(n=0)
    Graph(n=0) = Graph{Int}(n)

Construct a `Graph` with `n` vertices and no edges.

    Graph{T}(adjmx::AbstractMatrix; upper=false, selfedges=true)

Construct a `Graph{T}` from the adjacency matrix `adjmx`, placing an edge in
correspondence to each nonzero element of `adjmx`.
If `selfedges=false` the diagonal elements of `adjmx` are ignored.
If `upper=true` only the upper triangular part of `adjmx` is considered.
"""
mutable struct Graph{T<:Integer} <: AGraph
    ne::Int
    fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst)
                                # fadjlist is sorted
end

function Graph{T}(n::Integer = 0) where T<:Integer
    fadjlist = [Vector{T}() for _=1:n]
    return Graph{T}(0, fadjlist)
end

Graph{T}(n::Integer, m::Integer; seed = -1) where {T<:Integer} =
    erdos_renyi(n, m, Graph{T}; seed=seed)
Graph(n::T, m::T; kws...) where {T<:Integer} = Graph{T}(n, m; kws...)
Graph(n::T) where {T<:Integer} = Graph{T}(n)
Graph() = Graph{Int}()

"""
    mutable struct DiGraph{T<:Integer} <: ADiGraph
        ne::Int
        fadjlist::Vector{Vector{T}}
        badjlist::Vector{Vector{T}}
    end

A simple digraph type based on two adjacency lists (forward and backward).

    DiGraph{T}(n=0)
    DiGraph(n=0) = DiGraph{Int}(n)

Construct a `DiGraph` with `n` vertices and no edges.

    DiGraph{T}(adjmx::AbstractMatrix; selfedges=true)

Construct a `DiGraph{T}` from the adjacency matrix `adjmx`, placing an edge in
correspondence to each nonzero element of `adjmx`.
If `selfedges=false` the diagonal elements of `adjmx` are ignored.
"""
mutable struct DiGraph{T<:Integer} <: ADiGraph
    ne
    fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst)
    badjlist::Vector{Vector{T}} # [dst]: (src, src, src)
                                # fadjlist and badjlist are sorted
end

function DiGraph{T}(n::Integer = 0) where T<:Integer
    fadjlist = [Vector{T}() for _=1:n]
    badjlist = [Vector{T}() for _=1:n]
    return DiGraph{T}(0, fadjlist, badjlist)
end

DiGraph{T}(n::Integer, m::Integer; seed = -1) where {T<:Integer} =
    erdos_renyi(n, m, DiGraph{T}; seed=seed)
DiGraph(n::T, m::T; kws...) where {T<:Integer} = DiGraph{T}(n, m; kws...)
DiGraph(n::T) where {T<:Integer} = DiGraph{T}(n)
DiGraph() = DiGraph{Int}()

const GraphOrDiGraph{T} = Union{Graph{T}, DiGraph{T}}

edgetype(::Type{DiGraph{T}}) where {T} = Edge{T}
edgetype(::Type{Graph{T}}) where {T} = Edge{T}
graphtype(::Type{DiGraph{T}}) where {T} = Graph{T}
graphtype(::Type{DiGraph}) = Graph{Int}
graphtype(::Type{Graph}) = Graph{Int}
digraphtype(::Type{DiGraph}) = DiGraph{Int}
digraphtype(::Type{Graph}) = DiGraph{Int}
digraphtype(::Type{Graph{T}}) where {T} = DiGraph{T}
vertextype(::Type{Graph{T}}) where {T} = T
vertextype(::Type{DiGraph{T}}) where {T} = T

nv(g::GraphOrDiGraph{T}) where {T} = T(length(g.fadjlist))
ne(g::GraphOrDiGraph) = g.ne

pop_vertex!(g::Graph) = (clean_vertex!(g, nv(g)); pop!(g.fadjlist); nv(g)+1)
pop_vertex!(g::DiGraph) = (clean_vertex!(g, nv(g)); pop!(g.fadjlist); pop!(g.badjlist); nv(g)+1)

function add_edge!(g::Graph, s::Integer, d::Integer)
    E = edgetype(g)
    e = s <= d ? E(s, d) : E(d, s)
    (s in vertices(g) && d in vertices(g)) || return (false, e)
    inserted = _insert_and_dedup!(g.fadjlist[s], d)
    if inserted
        g.ne += 1
    end
    if s != d
        inserted = _insert_and_dedup!(g.fadjlist[d], s)
    end
    return (inserted, e)
end

function rem_edge!(g::Graph, u, v)
    i = searchsorted(g.fadjlist[u], v)
    length(i) > 0 || return false   # edge not in graph
    deleteat!(g.fadjlist[u], i[1])
    if u != v     # not a self loop
        i2 = searchsorted(g.fadjlist[v], u)[1]
        deleteat!(g.fadjlist[v], i2)
    end
    g.ne -= 1
    return true # edge successfully removed
end

function add_vertex!(g::Graph{T}) where T
    push!(g.fadjlist, Vector{T}())
    return nv(g)
end

function copy(g::DiGraph{T}) where T
    return DiGraph{T}(g.ne, deepcopy(g.fadjlist), deepcopy(g.badjlist))
end

function add_edge!(g::DiGraph, s::Integer, d::Integer)
    E = edgetype(g)
    e = E(s, d)
    (s in vertices(g) && d in vertices(g)) || return (false, e)
    inserted = _insert_and_dedup!(g.fadjlist[s], d)
    g.ne = ifelse(inserted, g.ne+1, g.ne)
    inserted && !_insert_and_dedup!(g.badjlist[d], s)
    return (inserted, e)
end

function rem_edge!(g::DiGraph, u, v)
    has_edge(g,u,v) || return false
    i = searchsorted(g.fadjlist[u],v)[1]
    deleteat!(g.fadjlist[u], i)
    i = searchsorted(g.badjlist[v], u)[1]
    deleteat!(g.badjlist[v], i)
    g.ne -= 1
    return true
end

function add_vertex!(g::DiGraph{T}) where T
    push!(g.badjlist, Vector{T}())
    push!(g.fadjlist, Vector{T}())
    return nv(g)
end

function reverse(g::DiGraph{T}) where T
    gnv = nv(g)
    gne = ne(g)
    h = DiGraph{T}(gnv)
    h.fadjlist = deepcopy(g.badjlist)
    h.badjlist = deepcopy(g.fadjlist)
    h.ne = gne

    return h
end

function reverse!(g::DiGraph)
    g.fadjlist, g.badjlist = g.badjlist, g.fadjlist
    return g
end

out_neighbors(g::GraphOrDiGraph,v) = g.fadjlist[v]
in_neighbors(g::DiGraph,v) = g.badjlist[v]

edge(g::DiGraph{T}, u, v) where {T} = Edge{T}(u, v)
edge(g::Graph{T}, u, v) where {T} = Edge{T}(u, v)
# edge(g::Graph, u, v) = u <= v ? Edge(u, v) : Edge(v, u)

#### fallbaks override #######

digraph(g::Graph{T}) where {T} = DiGraph{T}(2ne(g), deepcopy(g.fadjlist), deepcopy(g.fadjlist))

out_adjlist(g::GraphOrDiGraph) = g.fadjlist
in_adjlist(g::DiGraph) = g.badjlist

==(g::G, h::G) where {G<:GraphOrDiGraph} = nv(g) == nv(h) &&
                ne(g) == ne(h) && g.fadjlist == h.fadjlist

function has_edge(g::Graph, u, v)
    (u > nv(g) || v > nv(g)) && return false
    if degree(g, u) > degree(g, v)
        u, v = v, u
    end
    return length(searchsorted(neighbors(g,u), v)) > 0
end

function has_edge(g::DiGraph, u, v)
    (u > nv(g) || v > nv(g)) && return false
    if out_degree(g,u) < in_degree(g,v)
        return length(searchsorted(out_neighbors(g,u), v)) > 0
    else
        return length(searchsorted(in_neighbors(g,v), u)) > 0
    end
end

# UNSAFE METHODS

function unsafe_add_edge!(g::Graph, s, d)
    push!(g.fadjlist[s], d)
    s != d && push!(g.fadjlist[d], s)
    g.ne += 1
end

#TODO should check for duplicate edges.
## sort!(union!(x)) ??
# shold verify i~j  => j~i
function rebuild!(g::Graph)
    for neigs in g.fadjlist
        sort!(neigs)
    end
end

function unsafe_add_edge!(g::DiGraph, s, d)
    push!(g.fadjlist[s], d)
    push!(g.badjlist[d], s)
    g.ne += 1
end

function rebuild!(g::DiGraph)
    for neigs in g.fadjlist
        sort!(neigs)
    end
    for neigs in g.badjlist
        sort!(neigs)
    end
end

function swap_vertices!(g::Graph, u::Integer, v::Integer)
    if u != v
        #TODO copying to avoid problems with self edges
        # maybe can copy only one of the two
        neigu = deepcopy(g.fadjlist[u])
        neigv = deepcopy(g.fadjlist[v])

        for j in neigu
            adj = g.fadjlist[j]
            kj = searchsortedfirst(adj, u)
            adj[kj] = v
            sort!(adj)
        end
        for j in neigv
            adj = g.fadjlist[j]
            kj = searchsortedfirst(adj, v)
            adj[kj] = u
            sort!(adj)
        end

        g.fadjlist[u], g.fadjlist[v] = g.fadjlist[v], g.fadjlist[u]
    end
end

function swap_vertices!(g::DiGraph, u::Integer, v::Integer)
    if u != v
        #TODO copying to avoid problems with self edges
        # maybe can copy only one of the two
        neigu = deepcopy(g.fadjlist[u])
        neigv = deepcopy(g.fadjlist[v])
        neiguin = deepcopy(g.badjlist[u])
        neigvin = deepcopy(g.badjlist[v])

        for j in neigu
            adj = g.badjlist[j]
            kj = searchsortedfirst(adj, u)
            adj[kj] = v
            sort!(adj)
        end

        for j in neigv
            adj = g.badjlist[j]
            kj = searchsortedfirst(adj, v)
            adj[kj] = u
            sort!(adj)
        end

        for j in neiguin
            adj = g.fadjlist[j]
            kj = searchsortedfirst(adj, u)
            adj[kj] = v
            sort!(adj)
        end

        for j in neigvin
            adj = g.fadjlist[j]
            kj = searchsortedfirst(adj, v)
            adj[kj] = u
            sort!(adj)
        end

        g.fadjlist[u], g.fadjlist[v] = g.fadjlist[v], g.fadjlist[u]
        g.badjlist[u], g.badjlist[v] = g.badjlist[v], g.badjlist[u]
    end
end


text:
module TestBalls

println("\nCollidingBalls: Demonstrating the use of allInstances to set up contact force between any number of balls")

using Modia
using ModiaMath

const k=10000
const d=0
const radius=0.05

springDamper(r, v) = k*r + d*v

function getForce(r, v, positions,
  velocities, contactLaw)
  # Calculate the forces on this body
  force = zeros(2)
  for i in 1:length(positions)
    pos = positions[i]
    vel = velocities[i]
    if r != pos
      delta = r - pos
      deltaV = v - vel
      f = if norm(delta) < 2*radius;
        -contactLaw((norm(delta)-
          2*radius)*delta/norm(delta), deltaV) else
        zeros(2) end
      force += f
    end
  end
  return force
end

@model Ball begin
  r = Var()
  v = Var()
  f = Var()
  m = 1.0
  event = Boolean(size=(), start=false)
@equations begin
  der(r) = v
  m*der(v) = f
  f = getForce(r, v, allInstances(r),
#    allInstances(v), (r,v) -> (k*r + d*v))
    allInstances(v), springDamper)
    event = Modia.Synchronous.positive(f[1])
  end
end

@model Balls1 begin
  b1 = Ball(r = Var(start=[0.0,2]), v = Var(start=[1,0]))
  b2 = Ball(r = Var(start=[0.5,2]), v = Var(start=[-1,0]))
  b3 = Ball(r = Var(start=[1.0,2]), v = Var(start=[0,0]))
end

#=
result = checkSimulation(Balls1, 1.5, "b3.r", 2.0, expandArrayIncidence=true, storeEliminated=false)
ModiaMath.plot(result, ("b1.r", "b2.r", "b3.r"), figure=1)
=#

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

@model Balls2 begin
  b1 = Ball(r = Var(start=[0.0,2]), v = Var(start=[1,0]))
  b2 = Ball(r = Var(start=[0.5,2]), v = Var(start=[0,0]))
  b3 = Ball(r = Var(start=[1.0,2]), v = Var(start=[-1,0]))
end

#=
result = checkSimulation(Balls2, 1.5, "b3.r", 2.0, expandArrayIncidence=true, storeEliminated=false)
ModiaMath.plot(result, ("b1.r", "b2.r", "b3.r"), figure=2)
=#

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

@model Balls3 begin
  b1 = Ball(r = Var(start=[1.0,0]), v = Var(start=[1,0]))
  b2 = Ball(r = Var(start=[2.0,0]), v = Var(start=[0,0]))
  b3 = Ball(r = Var(start=[3.0,0]), v = Var(start=[0,0]))
  b4 = Ball(r = Var(start=[4.0,0]), v = Var(start=[0,0]))
  b5 = Ball(r = Var(start=[5.0,0]), v = Var(start=[0,0]))
  b6 = Ball(r = Var(start=[10.0,0]), v = Var(start=[-1,0]))
end

result = simulate(Balls3, 8, expandArrayIncidence=true, storeEliminated=false)
ModiaMath.plot(result, Tuple(["b$i.r" for i in 1:5]), figure=1, heading="Colliding balls", grid=false, figure=19)

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

#=
nBalls=4
@model Balls begin
#  balls = [Ball(r = Var(start=[i,2]), v = Var(start=if i==1; [1,0] else [0,0] end)) for i in 1:nBalls]
  balls = [Ball() for i in 1:nBalls]
end

result = simulate(Balls, 1.5, expandArrayIncidence=true, storeEliminated=false, logTranslation=true, logOnFile=false)
=#

end


text:
pos2box(x::Float64; w=2) = floor(Int,x/w)
pos2box(p::Vector{Float64}, dimv::Val; w=2) = ntuple(i->pos2box(p[i];w=w), dimv)

function find_neighbors(cellbox, root::Int; s=2)
	bools = isneighbor.(cellbox; rootbox=cellbox[root], s=s)
	bools[root] = false
	return bools
end
isneighbor(neighbox::NTuple{dim,Int}; rootbox::NTuple{dim,Int}, s=2) where dim = all(-s < neighbox[i] - rootbox[i] < s for i=1:dim)

function pushing!(tumor, root, cellbox, dimv::Val)
    queue = [root]
    while !isempty(queue)
        root = popfirst!(queue)
        r1 = tumor[root].position

        for n in shuffle!(findall(find_neighbors(cellbox, root; s=2)))
            r2 = tumor[n].position
            d = r2.-r1
            if norm(d) < 2.
                r2 .+= d*(2.2/norm(d) - 1)#/2
                #r1 .-= d*(2.2/norm(d) - 1)/2
                cellbox[n] = pos2box(r2, dimv)
                #cellbox[root] = pos2box(r1;dimv=Val(dim))
                push!(queue, n)
				#push!(queue, n)
            end
        end
    end
end

function pushing_recursive!(tumor, root, cellbox, dimv::Val)
    r1 = tumor[root].position

    N=length(tumor)
    for n in shuffle!(findall(find_neighbors(cellbox, root)))
        r2 = tumor[n].position
        d = r2.-r1
        if norm(d) < 2.
            r2 .+= d*(2.2/norm(d) - 1)
            cellbox[n] = pos2box(r2, dimv)

            pushing_recursive!(tumor, n, cellbox)
        end
    end
end

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
export InputLayer

@doc raw"""
The input layer specifies the input dimension to a flow model.

```julia
layer = InputLayer(3)
```
"""
struct InputLayer <: AbstractLayerPlaceholder
    dim :: Int
    function InputLayer(dim::Int)
        @assert dim > 1 "The specified input dimension should be larger than 1."
        return new(dim)
    end
end

getdim(layer::InputLayer) = layer.dim


text:
@testset "Refine with variable splitting factor (k)" begin
    @testset "E = $E" for E in 3:7
        @testset "rep = $rep" for rep in 1:10
        println("E=",E, "\trep=",rep)
        ###########################################
        # A TRIANGULATION FROM A RANDOM POINT CLOUD
        ###########################################
        n_points = 10
        tau = 1
        embedding = embedding_example(n_points+E*2, E, tau)
        points, simplex_inds = triangulate(embedding[1:size(embedding, 1), :])

        # Take the images of the simplices in the triangulation as the forward projections
        # in time of the original simplices. Append first point of triangulation to
        # 'fake' an invariant set where the last point falls within the convex hull
        # of the triangulation
        impoints = vcat(points[2:end, :], embedding[1, :].')
        ########################
        # ORIGINAL TRIANGULATION
        ########################
        centroids, radii = centroids_radii2(points, simplex_inds)
        centroids_im, radii_im = centroids_radii2(impoints, simplex_inds)
        volumes = simplex_volumes(points, simplex_inds)
        volumes_im = simplex_volumes(impoints, simplex_inds)

        t = Triangulation(points, impoints, simplex_inds, centroids, radii, volumes, centroids_im, radii_im, volumes_im)

        volume_before = sum(t.volumes)
        imvolume_before = sum(t.volumes_im)

        ###########################################################################
        # VERIFY THAT DESIRED TARGET RADIUS GIVES DESIRED REDUCTION IN SIMPLEX SIZE
        ###########################################################################
        target_radius = min(minimum(t.radii), minimum(t.radii_im)) * 0.5
        # query = query_refinement(t, target_radius)
        #
        # @testset "Predicted maximum radius is smaller than target_radius" begin
        #    @test maximum(query.resulting_radii) <= target_radius
        # end

        # ####################################
        # # PERFORM TRIANGULATION REFINEMENT
        # ####################################
        refine_variable_k!(t, target_radius)

        volume_after = sum(t.volumes)
        imvolume_after = sum(t.volumes_im)

        @testset "Splitting preserved volumes" begin
            @test volume_before ≈ volume_after atol = 0.00001
            @test imvolume_before ≈ imvolume_after atol = 0.00001
        end

    end
    end
end

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
using COSMO, Random, Test

include("./UnitTests/COSMOTestUtils.jl")

@testset "All Unit Tests" begin

  # Split tests into native interface unit tests and MOI unit tests
  include("./run_cosmo_tests.jl")
  include("./UnitTests/moi_wrapper.jl")

end
nothing


text:
# Utils for modifying plot style variables

function PushStyleColor(idx, col::UInt32)
    LibCImPlot.PushStyleColorU32(idx, col)
end

function PushStyleColor(idx, col::ImVec4)
    LibCImPlot.PushStyleColorVec4(idx, col)
end

function PushStyleVar(idx, val::AbstractFloat)
    LibCImPlot.PushStyleVarFloat(idx, Float32(val))
end

function PushStyleVar(idx, val::Integer)
    LibCImPlot.PushStyleVarInt(idx, Cint(val))
end


text:
## Connection and Linearization

## A connection between polynomial systems `P` and `Q` is a function `C(m,n)` satisfying
## p_n = ∑C(m,n) q_m or  pᵢ = ∑ Cʲᵢqⱼ (if unicode is helpful)
##
## From this a polynomial in P = ∑ aᵢ⋅pᵢ can be represented in a Q through
##
##  ∑ᵢ aᵢ⋅pᵢ = ∑ᵢ aᵢ  ∑ⱼ Cʲᵢqⱼ  = ∑ᵢ (∑ⱼ (aᵢ⋅Cʲᵢ)) qⱼ or ∑ⱼ (∑ᵢ aᵢCʲᵢ) qⱼ
##
## The choice of which depends on how the connection coefficients are presented.
##
## A connection between P and Q defines a `convert` method `Base.convert(Q,  p::P)`
## Conversion can also be achieved through polynomial evaluation. For example
## if `x=variable(Q)` then `Base.convert(Q,p) = p(x)` will also convert `p` into
## the basis Q. This is helpful when Q is a standard basis polynomial.
##
## However, if Q is not a standard basis polynomial, it isn't helpful unless
## one can multiply monomomials within Q. For standard basis polynomials this
## is known, as  xⁱ⋅xʲ = xⁱ⁺ʲ.
##
## For general Q though, one needs a linearization: a function α with
##
##  qᵢqⱼ = ∑_k α(k,i,j) q_k
##
## with a given α, multiplication can be defined. The linearization can be defined
## through composition: If P is the standard basis, connect Q to P, multiply ini P, then connect
## to Q.
##
## In short, if there is a linearization then multiplication can be defined and conversion defined through
## polynomial evaluation
## If there is a connection defined then conversion is defined and multiplication can be defined through composition
##

## For `AbstractCCOP` polynomials, Koepf and Schmersau present recursive formulas for the connection coefficients
## between Q and the standard basis, the standard and Q, and between P and Q when σ and σ̂ are the same. (This is the case when a polynomial system is parameterized, like Laguerre, Gegenbauer, Bessel, and  Jacobi.)

ConvertibleTypes = Union{AbstractCOP,Polynomials.StandardBasisPolynomial,FallingFactorial}
kn(::Type{P}, n) where {P<:Union{Polynomials.StandardBasisPolynomial,FallingFactorial}} =
    one(eltype(one(P)))
k1k0(::Type{P}, n) where {P<:Union{Polynomials.StandardBasisPolynomial,FallingFactorial}} =
    one(eltype(one(P)))
knk_1(::Type{P}, n) where {P<:Union{Polynomials.StandardBasisPolynomial,FallingFactorial}} =
    one(eltype(one(P)))

## If  the connection  coefficients for  (P,Q) satisfy c₀⋅ C̃ⁱⱼ + c₁⋅C̃ⁱ⁺¹ⱼ  + c₂⋅*C̃ⁱ⁺²ⱼ = 0 for some computable
##  values c₀,c₁,c₂ then this  will implement  `convert(Q,p::P)`
function _convert_cop(::Type{Q}, p::P) where {P<:ConvertibleTypes,Q<:ConvertibleTypes}
    d = degree(p)
    T, S = eltype(Q), eltype(p)  #
    R = typeof(one(T) * p[0] / 1)
    as = zeros(R, 1 + d)

    λⱼⱼ = one(R) * k0(P) / k0(Q)

    for j in 0:d
        λ = λⱼⱼ

        λⱼⱼ *= k1k0(P, j) / k1k0(Q, j)

        pⱼ = p[j]
        iszero(pⱼ) && continue

        C̃ʲⱼ = one(R)

        as[1 + j] += pⱼ * (λ * C̃ʲⱼ)

        C̃ⁱ⁺²ⱼ, C̃ⁱ⁺¹ⱼ = zero(R), C̃ʲⱼ  # for recursion starting  with C̃ʲ⁺¹ⱼ, C̃ʲⱼ
        for i in (j - 1):-1:0
            λ *= k1k0(Q, i)

            # c₀⋅ C̃ⁱⱼ + c₁⋅C̃ⁱ⁺¹ⱼ  + c₂⋅*C̃ⁱ⁺²ⱼ) = 0
            c₀, c₁, c₂ = connection_m(P, Q, i, j)
            C̃ⁱⱼ = iszero(c₀) ? zero(R) : -one(R) * (c₁ * C̃ⁱ⁺¹ⱼ + c₂ * C̃ⁱ⁺²ⱼ) / c₀
            as[1 + i] += pⱼ * (λ * C̃ⁱⱼ)

            C̃ⁱ⁺²ⱼ, C̃ⁱ⁺¹ⱼ = C̃ⁱ⁺¹ⱼ, C̃ⁱⱼ
        end
    end
    X = Polynomials.indeterminate(Q, p)
    ⟒(Q){R,X}(as)
end

# ## Use FastTransform for conversion, as  possible
# ## FastTransforms.kind2string.(0:15)
# # Legendre<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     T <: AbstractFloat,
#     Q <: Union{Legendre, OrthonormalLegendre},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}}
#     ps = coeffs(p)
#     Q(cheb2leg(ps, normcheb=isorthonormal(P), normleg=isorthonormal(Q)))
# end

# # Chebyshev<--Legendre
# function _convert(::Type{Q}, p::P) where {
#     T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Legendre{T}, OrthonormalLegendre{T}}
# }
#     ps =  coeffs(p)
#     Q(leg2cheb(ps, normleg=isorthonormal(P), normcheb=isorthonormal(Q)) )
# end

# # ultraspherical<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     β, α,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{β}, OrthonormalGegenbauer{β}},
#     P <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}}
# }

#     ps =  coeffs(p)
#     Q( ultra2ultra(ps, α, β, norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Jacobi<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     γ, δ, α, β,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Jacobi{α, β, T}, OrthonormalJacobi{α, β, T}}
# }

#     ps =  coeffs(p)
#     Q( jac2jac(ps, α, β, γ, δ,  norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Laguerre<--Laguerre
# function _convert(::Type{Q}, p::P) where {
#     α, β,  T <: AbstractFloat,
#     Q <: Union{Laguerre{β}, OrthonormalLaguerre{β}},
#     P <: Union{Laguerre{α, T}, OrthonormalLaguerre{α, T}}
# }

#     ps =  coeffs(p)
#     Q( lag2lag(ps, α, β,  norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Jacobi<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     γ, δ, α,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}}
# }

#     ps =  coeffs(p)
#     Q( ultra2jac(ps, α, γ, δ,  normultra=isorthonormal(P), normjac=isorthonormal(Q)) )

# end

# # ultraspherical<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     α, γ, δ,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}},
#     P <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}}
# }

#     ps =  coeffs(p)
#     Q( jac2ultra(ps,γ, δ, α, normjac=isorthonormal(P), normultra=isorthonormal(Q)) )

# end

# # Jacobi<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     γ, δ,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}
# }

#     ps =  coeffs(p)
#     Q( cheb2jac(ps,γ, δ, normcheb=isorthonormal(P), normjac=isorthonormal(Q)) )

# end

# # Chebyshev<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     γ, δ,  T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Jacobi{γ, δ,T}, OrthonormalJacobi{γ, δ,T}}
# }

#     ps =  coeffs(p)
#     Q( jac2cheb(ps,γ, δ, normjac = isorthonormal(P), normcheb=isorthonormal(Q)) )

# end

# # ultraspherical<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     α,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{α}, OrthonormalGegenbauer{α}},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}
# }

#     ps =  coeffs(p)
#     Q( cheb2ultra(ps,α, normcheb=isorthonormal(P), normultra=isorthonormal(Q)) )

# end

# # Chebyshev<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     α,  T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Gegenbauer{α,T}, OrthonormalGegenbauer{α,T}}
# }
#     ps =  coeffs(p)
#     Q( ultra2cheb(ps,α, normultra=isorthonormal(P), normcheb=isorthonormal(Q)) )

# end

##  Koepf and Schmersa Thm 2, p11 connection for P to Q when σ = σ̂
function connection_m(::Type{P}, ::Type{Q}, m, n) where {P<:AbstractCCOP,Q<:AbstractCCOP}
    a, b, c, d, e = abcde(P)
    ā, b̄, c̄, d̄, ē = abcde(Q)

    (a == ā && b == b̄ && c == c̄) || throw(ArgumentError("This connection assume σ = σ̄"))
    a², b², c², d², d̄², e², ē², m², n² = a^2, b^2, c^2, d^2, d̄^2, e^2, ē^2, m^2, n^2

    c0 =
        -(m - n) ⋅ (a ⋅ m + d - a + a ⋅ n) ⋅ (d̄ + 2a ⋅ m) ⋅ (d̄ + a + 2a ⋅ m) ⋅
        (d̄ + 3a + 2a ⋅ m)
    c0 *= (d̄ + 2a ⋅ m + 2a)^2

    c1 = -d ⋅ b ⋅ n ⋅ d̄ + 2d ⋅ a ⋅ m² ⋅ b + d ⋅ b ⋅ d̄ + 2d ⋅ a ⋅ m ⋅ b + 2d ⋅ ē ⋅ n ⋅ a
    c1 +=
        d ⋅ d̄ ⋅ ē + 2d ⋅ d̄ ⋅ b ⋅ m - m ⋅ b ⋅ d̄² - e ⋅ d̄² - 4a² ⋅ m² ⋅ e -
        m² ⋅ a ⋅ b ⋅ d̄ + b ⋅ n ⋅ d̄ ⋅ a - 2e ⋅ d̄ ⋅ a
    c1 +=
        -4a² ⋅ m ⋅ e - 4e ⋅ d̄ ⋅ a ⋅ m + 2m² ⋅ a² ⋅ ē + 2ē ⋅ a² ⋅ n² - 2ē ⋅ a² ⋅ n -
        m ⋅ a ⋅ b ⋅ d̄ + 2m ⋅ d̄ ⋅ ē ⋅ a
    c1 += 2m ⋅ ē ⋅ a² - b ⋅ n² ⋅ d̄ ⋅ a
    c1 *= (d̄ + 2a ⋅ m + 2a) ⋅ (m + 1) ⋅ (d̄ + a + 2a ⋅ m) ⋅ (d̄ + 3a + 2a ⋅ m)

    c2 =
        -(d̄ + 2a ⋅ m) ⋅ (m + 1) ⋅ (-a ⋅ m - 2a + a ⋅ n - d̄ + d) ⋅ (a ⋅ m + a ⋅ n + a + d̄)
    c2a =
        a ⋅ b² ⋅ m² - 4a² ⋅ m² ⋅ c - 8a² ⋅ m ⋅ c + 2a ⋅ m ⋅ b² - 4a ⋅ d̄ ⋅ m ⋅ c +
        m ⋅ b² ⋅ d̄ - 4a ⋅ d̄ ⋅ c - a ⋅ ē² + a ⋅ b² - c ⋅ d̄²
    c2a += b ⋅ ē ⋅ d̄ - 4 ⋅ a² ⋅ c + b² ⋅ d̄
    c2 *= c2a ⋅ (m + 2)

    return (c0, c1, c2)
end

##  Koepf and Schmersa Thm 4, p14 connection for P, Q=x^n
function connection_m(
    ::Type{P},
    ::Type{Q},
    m,
    n,
) where {P<:AbstractCCOP,Q<:Polynomials.StandardBasisPolynomial}
    a, b, c, d, e = abcde(P)

    c0 = (m - n) ⋅ (a ⋅ n + d - a + a ⋅ m)

    c1 = (m + 1) ⋅ (b ⋅ m + e)

    c2 = c ⋅ (m + 1) ⋅ (m + 2)

    return (c0, c1, c2)
end

##  Koepf and Schmersa Thm 5, p16 connection for P=x^n, Q
function connection_m(
    ::Type{P},
    ::Type{Q},
    m,
    n,
) where {P<:Polynomials.StandardBasisPolynomial,Q<:AbstractCCOP}
    ā, b̄, c̄, d̄, ē = abcde(Q)
    ā², b̄², c̄², d̄², ē², m², n² = ā^2, b̄^2, c̄^2, d̄^2, ē^2, m^2, n^2

    c0 =
        (n - m) ⋅ (d̄ + 2ā ⋅ m) ⋅ (d̄ + 3ā + 2ā ⋅ m) ⋅ (d̄ + ā + 2ā ⋅ m) ⋅
        (d̄ + 2ā ⋅ m + 2ā)^2

    c1 =
        (
            d̄ ⋅ ē + b̄ ⋅ d̄ + 2d̄ ⋅ b̄ ⋅ m + 2ā ⋅ m² ⋅ b̄ + 2ā ⋅ m ⋅ b̄ + 2ē ⋅ ā ⋅ n -
            d̄ ⋅ b̄ ⋅ n
        ) ⋅ (d̄ + 2ā ⋅ m + 2ā)
    c1 *= (m + 1) ⋅ (d̄ + 3ā + 2ā ⋅ m) ⋅ (d̄ + ā + 2ā ⋅ m)

    c2 = -(m + 2)
    c2a = -4ā² ⋅ c̄ ⋅ m²
    c2a +=
        ā ⋅ b̄² ⋅ m² + 2ā ⋅ b̄² ⋅ m - 4ā ⋅ c̄ ⋅ m ⋅ d̄ - 8ā² ⋅ c̄ ⋅ m + m ⋅ b̄² ⋅ d̄ -
        ā ⋅ ē^2 - d̄² ⋅ c̄ + b̄ ⋅ ē ⋅ d̄ - 4ā² ⋅ c̄
    c2a += -4ā ⋅ c̄ ⋅ d̄ + ā ⋅ b̄² + b̄² ⋅ d̄
    c2 *= c2a ⋅ (ā ⋅ m + ā ⋅ n + ā + d̄) ⋅ (m + 1) ⋅ (d̄ + 2ā ⋅ m)

    return (c0, c1, c2)
end

## With an explicit connectioni
function connection(
    ::Type{P},
    q::Q,
) where {P<:Polynomials.AbstractPolynomial,Q<:Polynomials.AbstractPolynomial}
    n = degree(q)
    T = eltype(q)
    T = eltype(one(eltype(one(P))) * one(q) / 1)
    cs = zeros(T, n + 1)

    for k in 0:n
        for (i, val) in Connection{P,Q}(n, k)
            cs[1 + k] = muladd(q[i], val, cs[1 + k])
        end
    end
    X = Polynomials.indeterminate(P, q)
    ⟒(P){T,X}(cs)
end


text:
function test_svd()
    println("-- Testing SVD-based recommender")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    recommender = SVD(data, 2)
    build!(recommender)

    # dimensionality reduction should preserve user-user/item-item similarities
    # i.e., recommendation list should be same as ItemKNN
    rec = recommend(recommender, 1, 4, [i for i in 1:8])
    @test first(rec[1]) == 8
    @test first(rec[2]) == 2
    @test first(rec[3]) == 5
    @test first(rec[4]) == 6
end

test_svd()


text:
type move
    move_number::Int64 # -1 if irrelevant
    move_type::String # "move", "drop" or "resign"
    sourcex::Int8 # -1 if no source
    sourcey::Int8 # -1 if no source
    targetx::Int8 # -1 if no target
    targety::Int8 # -1 if no target
    option::String # "" (empty string) if no option
    i_am_cheating::Bool # false by default, true if cheating
    targetx2::Int8 # -1 if no target2
    targety2::Int8 # -1 if no target2
    targetx3::Int8 # -1 if no target3
    targety3::Int8 # -1 if no target3
end

noMove = move(-1, "move", -1, -1, -1, -1, "", false, -1, -1, -1, -1)

type piece
    original::String # original piece name ("" for empty space) - constant throughout life of program
    name::String  # piece name ("" for empty space) - changes upon promotion
    x::Int8    # x-coordinate (0 for empty space)
    y::Int8    # y-coordinate (0 for empty space)
    side::Int8 # 0 for white; 1 for black; -1 for empty space
    promoted::Bool # false by default, true if promoted
end

empty = piece("", "", 0, 0, -1, false)

max_stars_repo_name: DarienFlamont/Shogi

max_issues_repo_name: DarienFlamont/Shogi

max_forks_repo_name: DarienFlamont/Shogi