text:
"""Wraps the data read from the text files with time series"""
struct RawTimeSeries
    initial_time::Dates.DateTime
    data::Dict
    length::Int
end

"""Describes how to construct time_series from raw time series data files."""
mutable struct TimeSeriesFileMetadata
    "User description of simulation"
    simulation::String
    "String version of abstract type for the component associated with the time series."
    category::String
    "Calling module should determine the actual type."
    "Name of time_series component"
    component_name::String
    "User-defined name"
    name::String
    "Controls normalization of time series.
     Use 1.0 for pre-normalized data.
     Use 'Max' to divide the time series by the max value in the column.
     Use any float for a custom scaling factor."
    normalization_factor::Union{String, Float64}
    "Path to the time series data file"
    data_file::String
    "Resolution of the data being parsed in seconds"
    resolution::Dates.Period
    percentiles::Vector{Float64}
    time_series_type::DataType
    "Calling module must set."
    component::Union{Nothing, InfrastructureSystemsComponent}
    "Applicable when data are scaling factors. Accessor function on component to apply to
    values."
    scaling_factor_multiplier::Union{Nothing, String}
    scaling_factor_multiplier_module::Union{Nothing, String}
end

function TimeSeriesFileMetadata(;
    simulation,
    category,
    component_name,
    name,
    normalization_factor,
    data_file,
    resolution,
    percentiles,
    time_series_type_module,
    time_series_type,
    scaling_factor_multiplier = nothing,
    scaling_factor_multiplier_module = nothing,
)
    return TimeSeriesFileMetadata(
        simulation,
        category,
        component_name,
        name,
        normalization_factor,
        data_file,
        resolution,
        percentiles,
        get_type_from_strings(time_series_type_module, time_series_type),
        nothing,
        scaling_factor_multiplier,
        scaling_factor_multiplier_module,
    )
end

"""Reads time_series metadata and fixes relative paths to the data files."""
function read_time_series_file_metadata(file_path::AbstractString)
    if endswith(file_path, ".json")
        metadata = open(file_path) do io
            metadata = Vector{TimeSeriesFileMetadata}()
            data = JSON3.read(io, Array)
            for item in data
                parsed_resolution = Dates.Millisecond(Dates.Second(item["resolution"]))
                normalization_factor = item["normalization_factor"]
                if !isa(normalization_factor, AbstractString)
                    normalization_factor = Float64(normalization_factor)
                end
                scaling_factor_multiplier =
                    get(item, "scaling_factor_multiplier", nothing)
                scaling_factor_multiplier_module =
                    get(item, "scaling_factor_multiplier_module", nothing)
                push!(
                    metadata,
                    TimeSeriesFileMetadata(;
                        simulation = item["simulation"],
                        category = item["category"],
                        component_name = item["component_name"],
                        name = item["name"],
                        normalization_factor = normalization_factor,
                        data_file = item["data_file"],
                        resolution = parsed_resolution,
                        # Use default values until CDM data is updated.
                        percentiles = get(item, "percentiles", []),
                        time_series_type_module = get(
                            item,
                            "module",
                            "InfrastructureSystems",
                        ),
                        time_series_type = get(item, "type", "SingleTimeSeries"),
                        scaling_factor_multiplier = scaling_factor_multiplier,
                        scaling_factor_multiplier_module = scaling_factor_multiplier_module,
                    ),
                )
            end
            return metadata
        end
    elseif endswith(file_path, ".csv")
        csv = DataFrames.DataFrame(CSV.File(file_path))
        metadata = Vector{TimeSeriesFileMetadata}()
        for row in eachrow(csv)
            category = row.category
            scaling_factor_multiplier = get(row, :scaling_factor_multiplier, nothing)
            scaling_factor_multiplier_module =
                get(row, :scaling_factor_multiplier_module, nothing)
            push!(
                metadata,
                TimeSeriesFileMetadata(;
                    simulation = row.simulation,
                    category = row.category,
                    component_name = row.component_name,
                    name = row.name,
                    resolution = Dates.Millisecond(Dates.Second(row.resolution)),
                    normalization_factor = row.normalization_factor,
                    data_file = row.data_file,
                    percentiles = [],
                    time_series_type_module = get(row, :module, "InfrastructureSystems"),
                    time_series_type = get(row, :type, "SingleTimeSeries"),
                    scaling_factor_multiplier = scaling_factor_multiplier,
                    scaling_factor_multiplier_module = scaling_factor_multiplier_module,
                ),
            )
        end

    else
        error("file not supported")
    end

    directory = dirname(file_path)
    for ts_metadata in metadata
        ts_metadata.data_file = abspath(joinpath(directory, ts_metadata.data_file))
    end

    return metadata
end

@scoped_enum NormalizationTypes MAX = 1

const NormalizationFactor = Union{Float64, NormalizationTypes}

function handle_normalization_factor(
    data::AbstractDict,
    normalization_factor::NormalizationFactor,
)
    for (k, v) in data
        data[k] = handle_normalization_factor(v, normalization_factor)
    end
    return data
end

get_max_value(ta::TimeSeries.TimeArray) = maximum(TimeSeries.values(ta))
get_max_value(ta::Vector) = maximum(ta)

function handle_normalization_factor(
    ta::Union{TimeSeries.AbstractTimeSeries, AbstractArray},
    normalization_factor::NormalizationFactor,
)
    if normalization_factor isa NormalizationTypes
        if normalization_factor == NormalizationTypes.MAX
            max_value = get_max_value(ta)
            ta = ta ./ max_value
        else
            error("support for normalization_factor=$normalization_factor not implemented")
        end
    else
        if normalization_factor != 1.0
            ta = ta ./ normalization_factor
        end
    end

    return ta
end

struct TimeSeriesParsedInfo
    simulation::String
    component::InfrastructureSystemsComponent
    name::String  # Component field on which time series data is based.
    normalization_factor::NormalizationFactor
    data::RawTimeSeries
    percentiles::Vector{Float64}
    file_path::String
    resolution::Dates.Period
    scaling_factor_multiplier::Union{Nothing, Function}

    function TimeSeriesParsedInfo(
        simulation,
        component,
        name,
        normalization_factor,
        data,
        percentiles,
        file_path,
        resolution,
        scaling_factor_multiplier = nothing,
    )
        new(
            simulation,
            component,
            name,
            normalization_factor,
            data,
            percentiles,
            abspath(file_path),
            resolution,
            scaling_factor_multiplier,
        )
    end
end

function TimeSeriesParsedInfo(metadata::TimeSeriesFileMetadata, raw_data::RawTimeSeries)
    if (
        metadata.scaling_factor_multiplier === nothing &&
        metadata.scaling_factor_multiplier_module !== nothing
    ) || (
        metadata.scaling_factor_multiplier !== nothing &&
        metadata.scaling_factor_multiplier_module === nothing
    )
        throw(
            DataFormatError(
                "scaling_factor_multiplier and scaling_factor_multiplier_module must both be set or not set",
            ),
        )
    end

    if metadata.scaling_factor_multiplier === nothing
        multiplier_func = nothing
    else
        multiplier_func = get_type_from_strings(
            metadata.scaling_factor_multiplier_module,
            metadata.scaling_factor_multiplier,
        )
    end

    if metadata.normalization_factor isa String
        if lowercase(metadata.normalization_factor) == "max"
            normalization_factor = NormalizationTypes.MAX
        else
            factor = metadata.normalization_factor
            throw(DataFormatError("unsupported normalization_factor {factor}"))
        end
    elseif metadata.normalization_factor == 0.0
        throw(DataFormatError("unsupported normalization_factor value of 0.0"))
    else
        normalization_factor = metadata.normalization_factor
    end

    return TimeSeriesParsedInfo(
        metadata.simulation,
        metadata.component,
        metadata.name,
        normalization_factor,
        raw_data,
        metadata.percentiles,
        metadata.data_file,
        metadata.resolution,
        multiplier_func,
    )
end

function make_time_array(info::TimeSeriesParsedInfo)
    return make_time_array(info.data, get_name(info.component), info.resolution)
end

function make_time_array(raw::RawTimeSeries, component_name, resolution)
    series_length = raw.length
    ini_time = raw.initial_time
    timestamps = range(ini_time; length = series_length, step = resolution)
    return TimeSeries.TimeArray(timestamps, raw.data[component_name])
end

struct TimeSeriesParsingCache
    time_series_infos::Vector{TimeSeriesParsedInfo}
    data_files::Dict{String, RawTimeSeries}
end

function TimeSeriesParsingCache()
    return TimeSeriesParsingCache(Vector{TimeSeriesParsedInfo}(), Dict{String, Any}())
end

function _add_time_series_info!(
    cache::TimeSeriesParsingCache,
    metadata::TimeSeriesFileMetadata,
)
    if !haskey(cache.data_files, metadata.data_file)
        cache.data_files[metadata.data_file] = read_time_series(metadata)
        @debug "Added time series file" _group = LOG_GROUP_TIME_SERIES metadata.data_file
    end
    return cache.data_files[metadata.data_file]
end


text:
@testset "show" begin
    root = "assets"

    sec = DemoSection(joinpath(root, "section", "default"))
    @test_reference joinpath("references", "section.txt") sec by=ignore_CR

    page = DemoPage(joinpath(root, "page", "default"))
    @test_reference joinpath("references", "page.txt") page by=ignore_CR
end

ext: jl


text:
cd(@__DIR__); include("setups/grid23x22.jl")
gr(dpi = 200)

##
frame = sgwt_frame(W; nf = 6)
x = 242
for j = 1:6
    plt = heatmap(reshape(frame[:, x, j], (Nx, Ny))', c = :viridis, ratio = 1,
            frame = :none, xlim = [1, Nx], size = (500, 400))
    savefig(plt, "../figs/Grid$(Nx)x$(Ny)_SGWT_frame_j$(j-1)_x$(x).png")
end


text:
using NCTiles,NCDatasets,NetCDF

# Point to interpolated 2D Data
examplesdir = joinpath("data","ex1")
selectfields = ["Chl050"]
indir = joinpath(examplesdir,"diags_interp",selectfields[1])
savedir = joinpath(examplesdir,"interp_ncfiles")
if ~ispath(savedir); mkpath(savedir); end
fnames = joinpath.(Ref(indir),filter(x -> occursin(".data",x), readdir(indir)))

# Hard coding these for now
prec = Float32
lon=-179.75:0.5:179.75; lat=-89.75:0.5:89.75;
n1,n2 = (length(lon),length(lat))
nsteps = 2
timeinterval = 3
time_steps = timeinterval/2:timeinterval:timeinterval*nsteps
longname = "Average chlorophyll concentration (top 50m)"
units = "mg Chl"

README = readlines(joinpath(examplesdir,"README"))

# Using NCDatasets

# Set up data for writing to NCfiles
# Define dimensions
dims = [NCvar("lon_c","degrees_east",size(lon),lon,Dict("long_name" => "longitude"),NCDatasets),
        NCvar("lat_c","degrees_north",size(lat),lat,Dict("long_name" => "latitude"),NCDatasets),
        NCvar("tim","days since 1992-01-01",Inf,time_steps,Dict(("long_name" => "tim","standard_name" => "time")),NCDatasets)
        ]

# Define field- BinData contains the filenames where the data sits so it's only loaded when needed
fielddata = BinData(fnames,prec,(n1,n2))
field = NCvar(selectfields[1],units,dims,fielddata,Dict("long_name" => longname),NCDatasets)

# Create the NetCDF file and populate with dimension and field info
ds,fldvar,dimlist = createfile(joinpath(savedir,"ex1_NCDatasets.nc"),field,README)

# Add field and dimension data
addData(fldvar,field)
addDimData.(Ref(ds),field.dims)

# Close the file
close(ds)

# Using NetCDF

# Set up data for writing to NCfiles
# Define dimensions
dims = [NCvar("lon_c","degrees east",size(lon),lon,Dict("long_name" => "longitude"),NetCDF),
        NCvar("lat_c","degrees north",size(lat),lat,Dict("long_name" => "latitude"),NetCDF),
        NCvar("tim","days from 1992-01-01",Inf,collect(time_steps),Dict(("long_name" => "time","standard_name" => "time")),NetCDF)
        ]

# Define field- BinData contains the filenames where the data sits so it's only loaded when needed
fielddata = BinData(fnames,prec,(n1,n2))
field = NCvar(selectfields[1],units,dims,fielddata,Dict("long_name" => longname),NetCDF)

using NetCDF

# Create the NetCDF file and populate with dimension and field info, as well as dimension data
ncfile,fldvar,dimlist = createfile(joinpath(savedir,"ex1_NetCDF.nc"),field,README)

# Add field data
addData(fldvar,field)

# Close the file
NetCDF.close(ncfile)


text:
using Xsum, Test, Random

Random.seed!(314159)
setprecision(65536)

@testset "xsum" begin
    for n in (100, 10^4)
        a = randn(n)
        exact = Float64(sum(big.(a)))
        @test exact == xsum(a)

        A = [a a]
        @test exact == xsum(@view A[:,1])
        @test exact == xsum(@view Array(A')[1,:])

        @test 2exact == xsum(2x for x in a) == xsum(x -> 2x, a)

        @test Complex(2exact, 4exact) == xsum(Complex(2x,4x) for x in a)
    end
end


text:
export DecSize, DecRank, Dec, VectorDec, MatrixDec
export checkndims, checklength, checkndims, checksize, checkrank, checkranks
export declength, decndims, decsize, decranks, decrank
export dec, dec!, vector, decrankselect!, decrankselect, factor!, factor, block!, block, decvcat, dechcat, decdcat
export decscale!, decreverse!, decmodetranspose!, decfill!, decrand!
export deczeros, decones, decrand
export decappend!, decprepend!, decpush!, decpushfirst!, decpop!, decpopfirst!, decinsert!, decdeleteat!
export decinsertidentity!
export decskp!, decskp, decmp, deckp, decadd, dechp
export decqr!, decsvd!

const DecSize = Matrix{Int}
const DecRank = Vector{Int}
const Dec{T,N} = Vector{Factor{T,N}} where {T<:Number,N}
const VectorDec{T} = Vector{VectorFactor{T}} where {T<:Number}
const MatrixDec{T} = Vector{MatrixFactor{T}} where {T<:Number}

function checkndims(d::Int)
	if d < 0
		throw(ArgumentError("d should be nonnegative"))
	end
end

function checklength(L::Int)
	if L < 0
		throw(ArgumentError("the number of factors should be nonnegative"))
	end
end

function checkndims(d::Vector{Int})
	L = length(d)
	if L == 0
		throw(ArgumentError("d is empty"))
	end
	if any(d .≤ 0)
		throw(BoundsError("the elements of d should be positive"))
	end
	if length(unique(d)) ≠ 1
		throw(ArgumentError("the values in the dimension vector are not identical"))
	end
end

function checksize(n::DecSize; len::Int=0, dim::Int=0)
	if length(n) == 0
		throw(ArgumentError("the size matrix is empty"))
	end
	if any(n .≤ 0)
		throw(ArgumentError("the elements of the size matrix should be positive"))
	end
	if len < 0
		throw(ArgumentError("the number of factors should be positive"))
	end
	if dim < 0
		throw(ArgumentError("the number of mode dimensions should be positive"))
	end
	if dim > 0 && size(n, 1) ≠ dim
		throw(ArgumentError("the number of rows in n is incorrect"))
	end
	if len > 0 && size(n, 2) ≠ len
		throw(ArgumentError("the number of columns in n is incorrect"))
	end
end

function checkrank(r::DecRank; len::Int=0)
	if length(r) < 2
		throw(ArgumentError("the rank vector should contain at least two elements"))
	end
	if any(r .≤ 0)
		throw(ArgumentError("the elements of the rank vector should be positive"))
	end
	if len > 0 && length(r) ≠ len+1
		throw(ArgumentError("the number of elements in the rank vector is incorrect"))
	end
end

function checkranks(p::Vector{Int}, q::Vector{Int}; len::Int=0)
	if length(p) ≠ length(q)
		throw(ArgumentError("the rank vectors should have the same length"))
	end
	if any([p q] .≤ 0)
		throw(ArgumentError("the elements of p and q should be positive"))
	end
	if p[2:end] ≠ q[1:end-1]
		throw(DimensionMismatch("the ranks are inconsistent"))
	end
	if len > 0 && length(p) ≠ len+1
		throw(ArgumentError("the number of elements in the rank vectors is incorrect"))
	end
end

function declength(U::Dec{T,N}) where {T<:Number,N}
	L = length(U)
	checklength(L)
	return L
end

function decndims(U::Dec{T,N}) where {T<:Number,N}
	# L = declength(U)
	# d = [ ndims(U[ℓ])-2 for ℓ ∈ 1:L ]
	d = [ factorndims(V) for V ∈ U ]
	checkndims(d)
	return d[1]
end

function decsize(U::Dec{T,N}) where {T<:Number,N}
	L = length(U)
	d = decndims(U)
	n = [ size(U[ℓ], 1+k) for k ∈ 1:d, ℓ ∈ 1:L ]
	n = n[:,:]
	checksize(n)
	return n
end

function decranks(U::Dec{T,N}) where {T<:Number,N}
	L = length(U)
	d = decndims(U)
	p = [ size(U[ℓ], 1) for ℓ ∈ 1:L ]
	q = [ size(U[ℓ], d+2) for ℓ ∈ 1:L ]
	return p,q
end

function decrank(U::Dec{T,N}) where {T<:Number,N}
	p,q = decranks(U)
	try checkranks(p,q) catch e
		isa(e, DimensionMismatch) && throw(DimensionMismatch("the factors of the decomposition have inconsistent ranks"))
	end
	return [p..., q[end]]
end

function dec(::Type{T}, d::Int, L::Int) where {T<:Number}
	checkndims(d)
	checklength(L)
	U = Vector{Array{T,d+2}}(undef, L)
	return U
end

dec(d::Int, L::Int) = dec(Float64, d, L)

function dec(::Type{T}, d::Int) where {T<:Number}
	checkndims(d)
	return Vector{Array{T,d+2}}(undef, 0)
end

dec(d::Int) = dec(Float64, d, L)

function dec(U::Dec{T,N}) where {T<:Number,N}
	return U
end

function dec(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number}
	checksize(n[:,:])
	nmat = isa(n, DecSize)
	nlen = size(n, 2)
	rvec = isa(r, DecRank)
	rlen = length(r)
	if any(r .≤ 0)
		throw(ArgumentError("the rank should be a positive integer or a vector of such"))
	end
	if rvec && (rlen < 2)
		throw(ArgumentError("when the rank parameter is specified as a vector, it should contain at least two entries"))
	end
	if rvec && nmat && rlen ≠ nlen+1
		throw(ArgumentError("when the size parameter is specified as a matrix with L columns and the rank parameter is specified as a vector, the latter should contain L+1 entries"))
	end
	if first < 0
		throw(ArgumentError("the first rank, when specified separately, should be positive"))
	end
	if last < 0
		throw(ArgumentError("the last rank, when specified separately, should be positive"))
	end
	if rvec && (first > 0 || last > 0)
		throw(ArgumentError("the first and last ranks are allowed to be specified separately only if the rank parameter is specified as an integer"))
	end
	if len < 0
		throw(ArgumentError("the number of factors, when specified, should be positive"))
	end
	if len > 0
		if nmat && len ≠ nlen
			throw(ArgumentError("when the number of rows in the size matrix is larger than one and the number of factors is specified, the two should be equal"))
		end
		if rvec && len+1 ≠ rlen
			throw(ArgumentError("when the rank parameter is specified as a vector with L+1 entries and the number of factors len is specified, it is required that len=L"))
		end
	end
	L = max(nlen, rlen-1, len)
	if !nmat
		n = repeat(n[:,:], outer=(1,L))
	end
	if !rvec
		r = r*ones(Int, L+1)
		if first > 0
			r[1] = first
		end
		if last > 0
			r[L+1] = last
		end
	end
	checkrank(r; len=L)
	U = [ Array{T}(undef, r[ℓ], n[:,ℓ]..., r[ℓ+1]) for ℓ ∈ 1:L ]
	return U
end

dec(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = dec(Float64, n, r; first=first, last=last, len=len)

function dec!(U::Factor{T,N}) where {T<:Number,N}
	return [U]
end

function dec(U::Factor{T,N}; len::Int=1) where {T<:Number,N}
	if len < 0
		throw(ArgumentError("the number of factors, when specified, should be positive"))
	end
	if len == 0
		return dec(T, N-2)
	end
	if len > 1
		p,q = factorranks(U)
		if p ≠ q
			throw(ArgumentError("the two ranks of U should be equal when the number of factors is specified as larger than one"))
		end
	end
	return [ copy(U) for ℓ ∈ 1:len ]
end

function vector(U::Dec{T,N}) where {T<:Number,N}
	return U
end

function decrankselect!(U::Dec{T,N}, α::Indices, β::Indices) where {T<:Number,N}
	# if isa(α, Int) || isa(β, Int)
	# 	throw(ArgumentError("for consistency with Base.selectdim, scalar α and β are not accepted; use α:α or β:β instead of α or β to select a subtensor of the factor whose first or second rank is one"))
	# end
	L = declength(U)
	r = decrank(U); p,q = r[1],r[L+1]
	α = indvec(α; min=1, max=p)
	β = indvec(β; min=1, max=q)
	if length(α) == 0
		throw(ArgumentError("the range for the first rank is empty"))
	end
	if α ⊈ 1:p
		throw(ArgumentError("the range for the first rank is incorrect"))
	end
	if length(β) == 0
		throw(ArgumentError("the range for the second rank is empty"))
	end
	if β ⊈ 1:q
		throw(ArgumentError("the range for the second rank is incorrect"))
	end
	U[L] = factorrankselect(U[L], :, β)
	U[1] = factorrankselect(U[1], α, :)
	return U
end

function decrankselect(U::Dec{T,N}, α::Indices, β::Indices) where {T<:Number,N}
	V = deepcopy(U)
	decrankselect!(V, α, β)
	return V
end

function factor!(U::Dec{T,N}; major::String="last") where {T<:Number,N}
	decskp!(U; major=major); U = U[1]
	return U
end

function factor(U::Dec{T,N}; major::String="last") where {T<:Number,N}
	V = deepcopy(U)
	V = factor!(V; major=major)
	return V
end

function block!(U::Dec{T,N}, α::Int, β::Int; major::String="last") where {T<:Number,N}
	L = declength(U)
	r = decrank(U); p = r[1]; q = r[L+1]
	if α ∉ 1:p
		throw(ArgumentError("the first rank index is out of range"))
	end
	if β ∉ 1:q
		throw(ArgumentError("the second rank index is out of range"))
	end
	decrankselect!(U, α:α, β:β)
	decskp!(U; major=major); U = U[1]
	U = block(U, 1, 1)
	return U
end

function block(U::Dec{T,N}, α::Int, β::Int; major::String="last") where {T<:Number,N}
	V = deepcopy(U)
	V = block!(V, α, β, major=major)
	return V
end

function decvcat(U::Dec{T,N}, V::Dec{T,N}, W::Vararg{Dec{T,N},M}) where {T<:Number,N,M}
	L = declength(U); m = decsize(U); p = decrank(U)
	W = (V,W...)
	for V ∈ W
		if declength(V) ≠ L
			throw(ArgumentError("the decompositions are incompatible in the number of factors"))
		end
		if decsize(V) ≠ m
			throw(ArgumentError("the decompositions are incompatible in mode size"))
		end
		q = decrank(V)
		if q[L+1] ≠ p[L+1]
			throw(ArgumentError("the decompositions are incompatible in the last rank"))
		end
	end
	return [ [ factordcat(U[ℓ], [ V[ℓ] for V ∈ W ]...) for ℓ ∈ 1:L-1 ]..., factorvcat(U[L], [ V[L] for V ∈ W ]...) ]
end

function dechcat(U::Dec{T,N}, V::Dec{T,N}, W::Vararg{Dec{T,N},M}) where {T<:Number,N,M}
	L = declength(U); m = decsize(U); p = decrank(U)
	W = (V,W...)
	for V ∈ W
		if declength(V) ≠ L
			throw(ArgumentError("the decompositions are incompatible in the number of factors"))
		end
		if decsize(V) ≠ m
			throw(ArgumentError("the decompositions are incompatible in mode size"))
		end
		q = decrank(V)
		if q[1] ≠ p[1]
			throw(ArgumentError("the decompositions are incompatible in the first rank"))
		end
	end
	return [ factorhcat(U[1], [ V[1] for V ∈ W ]...), [ factordcat(U[ℓ], [ V[ℓ] for V ∈ W ]...) for ℓ ∈ 2:L ]...]
end

function decdcat(U::Dec{T,N}, V::Dec{T,N}, W::Vararg{Dec{T,N},M}) where {T<:Number,N,M}
	L = declength(U); m = decsize(U)
	W = (V,W...)
	for V ∈ W
		if declength(V) ≠ L
			throw(ArgumentError("the decompositions are incompatible in the number of factors"))
		end
		if decsize(V) ≠ m
			throw(ArgumentError("the decompositions are incompatible in mode size"))
		end
	end
	return [ factordcat(U[ℓ], [ V[ℓ] for V ∈ W ]...) for ℓ ∈ 1:L ]
end

function decscale!(U::Dec{T,N}, α::T) where {T<:Number,N}
	L = declength(U)
	U[L] *= α
	return U
end

function decreverse!(W::Dec{T,N}) where {T<:Number,N}
	L = declength(W)
	reverse!(W)
	for ℓ ∈ 1:L
		W[ℓ] = factorranktranspose(W[ℓ])
	end
	return W
end

function decmodetranspose!(U::Dec{T,N}, τ::Union{NTuple{K,Int},Vector{Int}}) where {T<:Number,N,K}
	d = N-2
	if d == 0
		throw(ArgumentError("the decomposition should have at least one mode dimension"))
	end
	if length(τ) ≠ d || !isperm(τ)
		throw(ArgumentError("τ is not a valid permutation of the mode dimensions of U"))
	end
	isa(τ, Vector{Int}) || (τ = collect(τ))
	L = declength(U)
	for ℓ ∈ 1:L
		U[ℓ] = factormodetranspose(U[ℓ], τ)
	end
	return U
end

decmodetranspose!(U::Dec{T,N}) where {T<:Number,N} = decmodetranspose!(U, collect(decndims(U):-1:1))

function decfill!(U::Dec{T,N}, v::T) where {T<:Number,N}
	L = declength(U)
	for ℓ ∈ 1:L
		fill!(U[ℓ], v)
	end
	return U
end

function decrand!(rng::AbstractRNG, U::Dec{T,N}) where {T<:Number,N}
	L = declength(U)
	for ℓ ∈ 1:L
		rand!(rng, U[ℓ])
	end
	return U
end
decrand!(U::Dec{T,N}) where {T<:Number,N} = decrand!(Random.GLOBAL_RNG, U)

deczeros(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = decfill!(dec(T, n, r; first=first, last=last, len=len), convert(T, 0))
deczeros(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = deczeros(Float64, n, r; first=first, last=last, len=len)

decones(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = decfill!(dec(T, n, r; first=first, last=last, len=len), convert(T, 1))
decones(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = deczeros(Float64, n, r; first=first, last=last, len=len)

decrand(rng::AbstractRNG, ::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = decrand!(rng, dec(T, n, r; first=first, last=last, len=len))
decrand(rng::AbstractRNG, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = decrand(rng, Float64, n, r; first=first, last=last, len=len)
decrand(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = decrand(Random.GLOBAL_RNG, T, n, r; first=first, last=last, len=len)
decrand(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = decrand(Random.GLOBAL_RNG, Float64, n, r; first=first, last=last, len=len)

function decappend!(U::Dec{T,N}, V::Dec{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if decndims(U) ≠ decndims(V)
		throw(DimensionMismatch("U and V are inconsistent in the number of dimensions"))
	end
	p,q = decranks(U)
	r,s = decranks(V)
	if rankprecheck
		try checkranks(p,q) catch
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
		try checkranks(r,s) catch
			throw(ArgumentError("the factors of V have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		append!(p, r); append!(q, s);
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return append!(U, V)
end

function decprepend!(U::Dec{T,N}, V::Dec{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if decndims(U) ≠ decndims(V)
		throw(DimensionMismatch("U and V have different numbers of dimensions"))
	end
	p,q = decranks(U); r,s = decranks(V)
	if rankprecheck
		try checkranks(p,q) catch
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
		try checkranks(r,s) catch
			throw(ArgumentError("the factors of V have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		prepend!(p, r); prepend!(q, s);
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return prepend!(U, V)
end

function decpush!(U::Dec{T,N}, V::Factor{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if decndims(U) ≠ factorndims(V)
		throw(DimensionMismatch("U and V have different numbers of dimensions"))
	end
	p,q = decranks(U)
	r,s = factorranks(V)
	if rankprecheck
		try checkranks(p,q) catch
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		r,s = factorranks(V); push!(p, r); push!(q, s)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return push!(U, V)
end

function decpushfirst!(U::Dec{T,N}, V::Factor{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if decndims(U) ≠ factorndims(V)
		throw(DimensionMismatch("U and V have different numbers of dimensions"))
	end
	p,q = decranks(U)
	r,s = factorranks(V)
	if rankprecheck
		try
			checkranks(p,q)
		catch e
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		r,s = factorranks(V); pushfirst!(p, r); pushfirst!(q, s)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return pushfirst!(U, V)
end

function decpop!(U::Dec{T,N}) where {T<:Number,N}
	return pop!(U)
end

function decpopfirst!(U::Dec{T,N}) where {T<:Number,N}
	return popfirst!(U)
end

function decinsert!(U::Dec{T,N}, ℓ::Int, V::Factor{T,N}; path::String="", rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	L = declength(U)
	if path ∉ ("forward","backward")
		throw(ArgumentError("path should be either \"forward\" or \"backward\""))
	end
	if ℓ ∉ 1:L
		throw(ArgumentError("ℓ is required to be from 1:L, where L is the number of factors in U"))
	end
	(path == "forward") && (ℓ = ℓ+1)
	(ℓ == 1) && return decpushfirst!(U, V; rankprecheck=rankprecheck, rankpostcheck=rankpostcheck)
	(ℓ == L+1) && return decpush!(U, V; rankprecheck=rankprecheck, rankpostcheck=rankpostcheck)
	p,q = decranks(U)
	if rankprecheck
		try checkranks(p,q) catch
			throw(DimensionMismatch("the factors of U have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		r,s = factorranks(V); insert!(p, ℓ, r); insert!(q, ℓ, s)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return insert!(U, ℓ, V)
end

function decdeleteat!(U::Dec{T,N}, Λ::Union{Int,Vector{Int},NTuple{M,Int} where M}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	L = declength(U)
	if isa(Λ, Vector{Int}) && unique(Λ) ≠ Λ
		throw(ArgumentError("the entries of Λ should be unique"))
	end
	isa(Λ, Vector{Int}) || (Λ = collect(Λ))
	if Λ ⊈ 1:L
		throw(ArgumentError("Λ should be an element or a subset of 1:L with unique entries, where L is the number of factors in U"))
	end
	p,q = decranks(U)
	if rankprecheck
		try checkranks(p,q) catch
			throw(DimensionMismatch("the factors of U have inconsistent ranks"))
		end
	end
	if rankpostcheck
		deleteat!(p, Λ); deleteat!(q, Λ)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U are inconsistent for this operation"))
		end
	end
	return deleteat!(U, Λ)
end

function decinsertidentity!(U::Dec{T,N}, ℓ::Int; path::String="", rankprecheck::Bool=true) where {T<:Number,N}
	L = declength(U)
	if path ∉ ("forward","backward")
		throw(ArgumentError("path should be either \"forward\" or \"backward\""))
	end
	if ℓ ∉ 1:L
		throw(ArgumentError("ℓ is required to be from 1:L, where L is the number of factors in U"))
	end
	p,q = decranks(U)
	d = decndims(U)
	r = (path == "forward") ? q[ℓ] : p[ℓ]
	V = Matrix{T}(I, r, r); V = factor(V, ones(Int, d), [])
	(path == "forward") && (ℓ = ℓ+1)
	(ℓ == 1) && return decpushfirst!(U, V; rankprecheck=rankprecheck, rankpostcheck=false)
	(ℓ == L+1) && return decpush!(U, V; rankprecheck=rankprecheck, rankpostcheck=false)
	p,q = decranks(U)
	if rankprecheck
		try checkranks(p,q) catch
			throw(DimensionMismatch("the factors of U have inconsistent ranks"))
		end
	end
	return insert!(U, ℓ, V)
end

function decskp!(W::Dec{T,N}, Λ::Indices; path::String="", major::String="last") where {T<:Number,N}
	if path ∉ ("","forward","backward")
		throw(ArgumentError("the value of the keyword argument path should be \"\" (default, accepted only for empty Λ and for Λ=:), \"forward\" or \"backward\""))
	end
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	L = declength(W); decrank(W)
	if L == 0
		throw(ArgumentError("the decomposition is empty"))
	end
	if isa(Λ, Colon) && path ≠ ""
		throw(ArgumentError("when Λ is a colon, path should be omitted or specfied as \"\" (default)"))
	end
	isa(Λ, Colon) && (Λ = collect(2:L); path = "backward")
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	(length(Λ) == 0) && return W
	if path ∉ ("forward","backward")
		throw(ArgumentError("when Λ is neither empty nor a colon, path should be either \"forward\" or \"backward\""))
	end
	(path == "backward") || (Λ .+= 1; path = "backward")
	if unique(Λ) ≠ Λ
		throw(ArgumentError("Λ has duplicate entries"))
	end
	if Λ ⊈ 2:L
		throw(ArgumentError("Λ, when nonempty, should be a colon or an Int/Vector/NTuple/UnitRange of Int, from 1:L-1 for path=\"forward\" and from 2:L for path=\"backward\", where L is the number of factors in W"))
	end
	sort!(Λ; rev=true)
	for ℓ ∈ Λ
		W[ℓ-1] = factorcontract(W[ℓ-1], W[ℓ]; major=major)
		decdeleteat!(W, ℓ; rankprecheck=false, rankpostcheck=false)
	end
	return W
end

decskp!(W::Dec{T,N}; path::String="", major::String="last") where {T<:Number,N} = decskp!(W, :; path=path, major=major)

function decskp(W::Dec{T,N}, Λ::Indices; path::String="", major::String="last") where {T<:Number,N}
	if path ∉ ("","forward","backward")
		throw(ArgumentError("the value of the keyword argument path should be \"\" (default, accepted only for empty Λ and for Λ=:), \"forward\" or \"backward\""))
	end
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	L = declength(W); decrank(W)
	if L == 0
		throw(ArgumentError("the decomposition is empty"))
	end
	if isa(Λ, Colon) && path ≠ ""
		throw(ArgumentError("when Λ is a colon, path should be omitted or specfied as \"\" (default)"))
	end
	isa(Λ, Colon) && (Λ = collect(2:L); path = "backward")
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	if length(Λ) == 0
		throw(ArgumentError("Λ is empty"))
	end
	if (minimum(Λ):maximum(Λ)) ⊈ Λ
		throw(ArgumentError("the entries of Λ should form a set of contiguous integers"))
	end
	if path ∉ ("forward","backward")
		throw(ArgumentError("when Λ is neither empty nor a colon, path should be either \"forward\" or \"backward\""))
	end
	(path == "backward") || (Λ .+= 1; path = "backward")
	if unique(Λ) ≠ Λ
		throw(ArgumentError("Λ has duplicate entries"))
	end
	if Λ ⊈ 2:L
		throw(ArgumentError("Λ, when nonempty, should be a colon or a Vector/NTuple/UnitRange of Int, from 1:L-1 for path=\"forward\" and from 2:L for path=\"backward\", where L is the number of factors in W"))
	end
	M = length(Λ)
	sort!(Λ; rev=true)
	U = W[Λ[1]]
	for λ ∈ 2:M
		U = factorcontract(W[Λ[λ]], U; major=major)
	end
	return W
end

decskp(W::Dec{T,N}; path::String="", major::String="last") where {T<:Number,N} = decskp(W, :; path=path, major=major)

function decmp(U₁::Dec{T,N₁}, σ₁::Indices, U₂::Dec{T,N₂}, σ₂::Indices) where {T<:Number,N₁,N₂}
	n₁ = decsize(U₁); d₁ = decndims(U₁); L₁ = declength(U₁)
	n₂ = decsize(U₂); d₂ = decndims(U₂); L₂ = declength(U₂)
	if isa(σ₁, Vector{Any}) && length(σ₁) > 0
		throw(ArgumentError("if σ₁ is passed as a vector, it should be a vector of the type Vector{Int} or an empty vector of the type Vector{Any}"))
	end
	if isa(σ₂, Vector{Any}) && length(σ₂) > 0
		throw(ArgumentError("if σ₂ is passed as a vector, it should be a vector of the type Vector{Int} or an empty vector of the type Vector{Any}"))
	end
	isa(σ₁, Vector{Int}) || (σ₁ = indvec(σ₁; max=d₁))
	isa(σ₂, Vector{Int}) || (σ₂ = indvec(σ₂; max=d₂))
	if L₁ ≠ L₂
		throw(ArgumentError("U₁ and U₂ differ in the number of factors"))
	end
	L = L₁
	if length(σ₁) ≠ length(σ₂)
		throw(ArgumentError("the specified sets of modes of σ₁ and σ₂ are inconsistent"))
	end
	if σ₁ ⊈ 1:d₁ || unique(σ₁) ≠ σ₁
		throw(ArgumentError("the set of modes of U₁ is specified incorrectly"))
	end
	if σ₂ ⊈ 1:d₂ || unique(σ₂) ≠ σ₂
		throw(ArgumentError("the set of modes of U₂ is specified incorrectly"))
	end
	if n₁[σ₁,:] ≠ n₂[σ₂,:]
		throw(ArgumentError("U₁ and U₂ are inconsistent with respect to the specified modes"))
	end
	U = [ factormp(U₁[ℓ], σ₁, U₂[ℓ], σ₂) for ℓ ∈ 1:L ]
	return U
end

# function deckp(U::Dec{T,N}, V::Vararg{Dec{T,N},M}) where {T<:Number,N,M}
# 	n = decsize(U); L = declength(U)
# 	d = 1+length(V)
# 	for k ∈ 2:d
# 		if declength(V[k-1]) ≠ L
# 			throw(ArgumentError("the input decompositions differ in the number of factors"))
# 		end
# 	end
# 	return [ factorkp(U[ℓ], [ V[k-1][ℓ] for k ∈ 2:d ]...) for ℓ ∈ 1:L ]
# end

function deckp(U::Union{Dec{T,N},Tuple{Dec{T,N},Int}}, V::Vararg{Union{Dec{T,N},Tuple{Dec{T,N},Int}},M}) where {T<:Number,N,M}
	V = (U,V...)
	nf = length(V)
	U = Vector{Tuple{Dec{T,N},Int}}(undef, nf)
	for k ∈ 1:nf
		W = V[k]
		if isa(W, Tuple)
			if W[2] < 0
				throw(ArgumentError("all the specified exponents should be nonnegative"))
			end
			U[k] = W[1],W[2]
		else
			U[k] = W,1
		end
	end
	n = decsize(U[1][1]); L = declength(U[1][1])
	for k ∈ 2:nf
		if declength(U[k][1]) ≠ L
			throw(ArgumentError("the input decompositions differ in the number of factors"))
		end
	end
	return [ factorkp([ (U[k][1][ℓ],U[k][2]) for k ∈ 1:nf ]...) for ℓ ∈ 1:L ]
end

function decadd(U::Dec{T,N}, V::Dec{T,N}) where {T<:Number,N}
	m = decsize(U); L = declength(U)
	p = decrank(U); q = decrank(V)
	if declength(V) ≠ L
		throw(ArgumentError("U and V differ in the number of factors"))
	end
	if decsize(V) ≠ m
		throw(ArgumentError("U and V are inconsistent in mode size"))
	end
	if q[1] ≠ p[1]
		throw(ArgumentError("the decompositions are incompatible in the first rank"))
	end
	if q[L+1] ≠ p[L+1]
		throw(ArgumentError("the decompositions are incompatible in the last rank"))
	end
	(L == 1) && return U .+ V
	W = [ factorhcat(U[1], V[1]), [ factordcat(U[ℓ], V[ℓ]) for ℓ ∈ 2:L-1 ]..., factorvcat(U[L], V[L]) ]
	return W
end

function dechp(U::Dec{T,N}, V::Dec{T,N}) where {T<:Number,N}
	m = decsize(U); L = declength(U)
	if declength(V) ≠ L
		throw(ArgumentError("U and V differ in the number of factors"))
	end
	if decsize(V) ≠ m
		throw(ArgumentError("U and V are inconsistent in mode size"))
	end
	W = [ factorhp(U[ℓ], V[ℓ]) for ℓ ∈ 1:L ]
	return W
end

function decqr!(W::Dec{T,N}, Λ::Indices; pivot::Bool=false, path::String="") where {T<:FloatRC,N}
	L = declength(W); decrank(W)
	if L == 0
		throw(ArgumentError("the decomposition is empty"))
	end
	if path ∉ ("","forward","backward")
		throw(ArgumentError("path should be either \"\" (default, accepted only when path can be deduced from Λ), \"forward\" or \"backward\""))
	end
	if path == "" && (isa(Λ, Colon) || length(unique(Λ)) == 1)
		throw(ArgumentError("path cannot be deduced from Λ and should be specified as either \"forward\" or \"backward\""))
	end
	isa(Λ, Colon) && (Λ = (path == "forward") ? (1:L) : (L:-1:1))
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	(length(Λ) == 0) && return W
	if unique(Λ) ≠ Λ
		throw(ArgumentError("Λ has duplicate entries"))
	end
	if (minimum(Λ):maximum(Λ)) ⊈ Λ
		throw(ArgumentError("the entries of Λ should form a set of contiguous integers"))
	end
	if Λ ⊈ 1:L
		throw(ArgumentError("Λ is out of range"))
	end
	fw,bw = issorted(Λ, rev=false),issorted(Λ, rev=true)
	if !fw && !bw
		throw(ArgumentError("Λ is not sorted in ascending or descending order, so it is inconsistent with any forward or backward path"))
	end
	@assert !(fw && bw && path == "")
	fw && path == "" && (path = "forward")
	bw && path == "" && (path = "backward")
	if path == "forward" && !fw
		throw(ArgumentError("Λ is not sorted in ascending order, so it is inconsistent with any forward path"))
	end
	if path == "backward" && !bw
		throw(ArgumentError("Λ is not sorted in descending order, so it is inconsistent with any backward path"))
	end
	M = length(Λ)
	for λ ∈ 1:M
		ℓ = Λ[λ]
		W[ℓ],R = factorqr!(W[ℓ], Val(pivot); rev=(path == "backward"))
		if λ < M
			ν = Λ[λ+1]
			W[ν] = factorcontract(R, W[ν], rev=(path == "backward"))
		else
			decinsert!(W, ℓ, R; path=path, rankprecheck=false, rankpostcheck=true)
		end
	end
	return W
end

decqr!(W::Dec{T,N}; pivot::Bool=false, path::String="") where {T<:FloatRC,N} = decqr!(W, :; pivot=pivot, path=path)

function decsvd!(W::Dec{T,N}, Λ::Indices, n::Union{Colon,DecSize}; path::String="", aTol::Float2=0.0, aTolDistr::Float2=0.0, rTol::Float2=0.0, rTolDistr::Float2=0.0, rank::Int2=0, major::String="last") where {T<:FloatRC,N}
	# assumes that the decomposition is orthogonal
	L = declength(W); decrank(W)
	if L == 0
		throw(ArgumentError("the decomposition is empty"))
	end
	if path ∉ ("","forward","backward")
		throw(ArgumentError("path should be either \"\" (default, accepted only when path can be deduced from Λ), \"forward\" or \"backward\""))
	end
	if path == "" && (isa(Λ, Colon) || length(unique(Λ)) == 1)
		throw(ArgumentError("path cannot be deduced from Λ and should be specified as either \"forward\" or \"backward\""))
	end
	isa(Λ, Colon) && (Λ = (path == "forward") ? (1:L) : (L:-1:1))
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	(length(Λ) == 0) && return W
	if (minimum(Λ):maximum(Λ)) ⊈ Λ
		throw(ArgumentError("the entries of Λ should form a set of contiguous integers"))
	end
	if Λ ⊈ 1:L
		throw(ArgumentError("Λ is out of range"))
	end
	fw,bw = issorted(Λ, rev=false),issorted(Λ, rev=true)
	if !fw && !bw
		throw(ArgumentError("Λ is not sorted in ascending or descending order, so it is inconsistent with any forward or backward path"))
	end
	fw && path == "" && (path = "forward")
	bw && path == "" && (path = "backward")
	if path == "forward" && !fw
		throw(ArgumentError("Λ is not sorted in ascending order, so it is inconsistent with any forward path"))
	end
	if path == "backward" && !bw
		throw(ArgumentError("Λ is not sorted in descending order, so it is inconsistent with any backward path"))
	end
	m = decsize(W)
	isa(n, Colon) && (n = m[:,Λ])
	if size(n,1) ≠ size(m,1)
		throw(ArgumentError("the number of rows in n should be equal to the number of dimensions in each factor of W"))
	end
	if size(n,2) ≠ length(Λ)
		throw(ArgumentError("the number of columns in n should be equal to the number of elements in Λ"))
	end
	for ℓ ∈ unique(Λ)
		λ = (Λ .== ℓ)
		if m[:,ℓ] ≠ prod(n[:,λ]; dims=2)[:]
			throw(ArgumentError("n and Λ are incompatible with the size of the factors of W"))
		end
	end
	K = length(Λ)
	#
	if any(aTol .< 0)
		throw(ArgumentError("aTol should be a nonnegative Float64 or a vector of such"))
	end
	if isa(aTol, Float64)
		aTol = aTol*ones(K)/sqrt(K)
	elseif length(aTol) ≠ K
		throw(ArgumentError("aTol, passed as a vector, has incorrect length"))
	end

	if any(aTolDistr .< 0)
		throw(ArgumentError("aTolDistr should be a nonnegative Float64 or a vector of such"))
	end
	if isa(aTolDistr, Float64)
		aTolDistr = aTolDistr*ones(K)/sqrt(K)
	elseif length(aTolDistr) ≠ K
		throw(ArgumentError("aTolDistr, passed as a vector, has incorrect length"))
	end

	if any(rTol .< 0)
		throw(ArgumentError("rTol should be a nonnegative Float64 or a vector of such"))
	end
	if isa(rTol, Float64)
		rTol = rTol*ones(K)/sqrt(K)
	elseif length(rTol) ≠ K
		throw(ArgumentError("rTol, passed as a vector, has incorrect length"))
	end

	if any(rTolDistr .< 0)
		throw(ArgumentError("rTolDistr should be a nonnegative Float64 or a vector of such"))
	end
	if isa(rTolDistr, Float64)
		rTolDistr = rTolDistr*ones(K)/sqrt(K)
	elseif length(rTolDistr) ≠ K
		throw(ArgumentError("rTolDistr, passed as a vector, has incorrect length"))
	end

	if any(rank .< 0)
		throw(ArgumentError("rank should be a nonnegative Int or a vector of such"))
	end
	if isa(rank, Int)
		rank = rank*ones(Int, K)
	elseif length(rank) ≠ K
		throw(ArgumentError("rank, passed as a vector, has incorrect length"))
	end
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end

	ε = zeros(Float64, K); δ = zeros(Float64, K)
	σ = Vector{Vector{Float64}}(undef, K)
	ρ = zeros(Int, K)
	μ = 0.0
	aTolAcc = 0.0; rTolAcc = 0.0
	for λ ∈ 1:K
		ℓ = Λ[λ]
		if λ == 1
			ε₁ = [aTol[λ],aTolDistr[λ]]; ε₁ = ε₁[ε₁ .> 0]; ε₁ = isempty(ε₁) ? 0.0 : minimum(ε₁)
			δ₁ = [rTol[λ],rTolDistr[λ]]; δ₁ = δ₁[δ₁ .> 0]; δ₁ = isempty(δ₁) ? 0.0 : minimum(δ₁)
			U,V,ε[1],δ[1],μ,ρ[1],σ[1] = factorsvd!(W[ℓ], n[:,λ], :; aTol=ε₁, rTol=δ₁, rank=rank[λ], rev=(path == "backward"), major=major)
		else
			(aTolDistr[λ] > 0) && (aTolDistr[λ] = sqrt(aTolDistr[λ]^2+aTolAcc^2); aTolAcc = 0.0)
			(rTolDistr[λ] > 0) && (rTolDistr[λ] = sqrt(rTolDistr[λ]^2+rTolAcc^2); rTolAcc = 0.0)
			ε₁ = [aTol[λ],aTolDistr[λ],μ*rTol[λ],μ*rTolDistr[λ]]; ε₁ = ε₁[ε₁ .> 0]
			ε₁ = isempty(ε₁) ? 0.0 : minimum(ε₁)
			U,V,ε[λ],_,_,ρ[λ],σ[λ] = factorsvd!(W[ℓ], n[:,λ], :; aTol=ε₁, rank=rank[λ], rev=(path == "backward"), major=major)
			δ[λ] = (μ > 0) ? ε[λ]/μ : 0.0
		end
		W[ℓ] = U
		if λ < K && Λ[λ+1] ≠ ℓ
			ν = Λ[λ+1]
			W[ν] = factorcontract(V, W[ν], rev=(path == "backward"), major=major)
		else
			decinsert!(W, ℓ, V; path=path, rankprecheck=false, rankpostcheck=true)
			(path == "forward") && (Λ[λ:end] .+= 1)
		end
	end
	return W,ε,δ,μ,ρ,σ
end

decsvd!(W::Dec{T,N}, Λ::Indices; path::String="", aTol::Float2=0.0, aTolDistr::Float2=0.0, rTol::Float2=0.0, rTolDistr::Float2=0.0, rank::Int2=0, major::String="last") where {T<:FloatRC,N} = decsvd!(W, Λ, :; path=path, aTol=aTol, aTolDistr=aTolDistr, rTol=rTol, rTolDistr=rTolDistr, rank=rank, major=major)

decsvd!(W::Dec{T,N}; path::String="", aTol::Float2=0.0, aTolDistr::Float2=0.0, rTol::Float2=0.0, rTolDistr::Float2=0.0, rank::Int2=0, major::String="last") where {T<:FloatRC,N} = decsvd!(W, :, :; path=path, aTol=aTol, aTolDistr=aTolDistr, rTol=rTol, rTolDistr=rTolDistr, rank=rank, major=major)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Abstract optimizer types #
# ========================== #
"""
    AbstractStructuredOptimizer

Abstract supertype for structure-exploiting optimizers.
"""
abstract type AbstractStructuredOptimizer <: MOI.AbstractOptimizer end
"""
    AbstractStructuredOptimizerAttribute

Abstract supertype for attribute objects that can be used to set or get attributes (properties) of the structure-exploiting optimizer.
"""
abstract type AbstractStructuredOptimizerAttribute <: MOI.AbstractOptimizerAttribute end
"""
    AbstractSampledOptimizer

Abstract supertype for sample-based optimizers.
"""
abstract type AbstractSampledOptimizer <: MOI.AbstractOptimizer end
"""
    AbstractSampledOptimizerAttribute

Abstract supertype for attribute objects that can be used to set or get attributes (properties) of the sample-based optimizer.
"""
abstract type AbstractSampledOptimizerAttribute <: MOI.AbstractOptimizerAttribute end

"""
    AbstractCrash

Abstract supertype for crash methods.
"""
abstract type AbstractCrash end

StochasticProgramOptimizerType = Union{MOI.AbstractOptimizer, AbstractStructuredOptimizer}

# StochasticProgramOptimizer #
# ========================== #
"""
    StochasticProgramOptimizer

Wrapper type around both the optimizer_constructor provided to a stochastic program and the resulting optimizer object. Used to conviently distinguish between standard MOI optimizers and structure-exploiting optimizers when instantiating the stochastic program.
"""
mutable struct StochasticProgramOptimizer{Opt <: StochasticProgramOptimizerType}
    optimizer_constructor
    optimizer::Opt

    function StochasticProgramOptimizer(::Nothing)
        universal_fallback = MOIU.UniversalFallback(MOIU.Model{Float64}())
        caching_optimizer = MOIU.CachingOptimizer(universal_fallback, MOIU.AUTOMATIC)
        return new{StochasticProgramOptimizerType}(nothing, caching_optimizer)
    end

    function StochasticProgramOptimizer(optimizer_constructor, optimizer::MOI.AbstractOptimizer)
        universal_fallback = MOIU.UniversalFallback(MOIU.Model{Float64}())
        caching_optimizer = MOIU.CachingOptimizer(universal_fallback, MOIU.AUTOMATIC)
        MOIU.reset_optimizer(caching_optimizer, optimizer)
        Opt = MOI.AbstractOptimizer
        return new{Opt}(optimizer_constructor, caching_optimizer)
    end

    function StochasticProgramOptimizer(optimizer_constructor, optimizer::AbstractStructuredOptimizer)
        Opt = AbstractStructuredOptimizer
        return new{Opt}(optimizer_constructor, optimizer)
    end
end

function StochasticProgramOptimizer(optimizer_constructor)
    optimizer = MOI.instantiate(optimizer_constructor; with_bridge_type = bridge_type(optimizer_constructor))
    return StochasticProgramOptimizer(optimizer_constructor, optimizer)
end

bridge_type(optimizer::Type{<:AbstractStructuredOptimizer}) = nothing
bridge_type(optimizer::Type{<:AbstractSampledOptimizer}) = nothing
bridge_type(optimizer::Type{<:MOI.AbstractOptimizer}) = Float64
bridge_type(optimizer::MOI.OptimizerWithAttributes) = bridge_type(optimizer.optimizer_constructor)
bridge_type(optimizer::Function) = bridge_type(typeof(optimizer()))

function has_provided_optimizer(sp_optimizer::StochasticProgramOptimizer)
    return sp_optimizer.optimizer_constructor !== nothing
end

function check_provided_optimizer(sp_optimizer::StochasticProgramOptimizer)
    if !has_provided_optimizer(sp_optimizer)
        throw(NoOptimizer())
    end
end

function master_optimizer(sp_optimizer::StochasticProgramOptimizer)
    return master_optimizer(sp_optimizer, sp_optimizer.optimizer)
end

function subproblem_optimizer(sp_optimizer::StochasticProgramOptimizer)
    return subproblem_optimizer(sp_optimizer, sp_optimizer.optimizer)
end

function master_optimizer(sp_optimizer::StochasticProgramOptimizer, ::MOI.AbstractOptimizer)
    return sp_optimizer.optimizer_constructor
end

function subproblem_optimizer(sp_optimizer::StochasticProgramOptimizer, optimizer::MOI.AbstractOptimizer)
    return master_optimizer(sp_optimizer, optimizer)
end

function master_optimizer(sp_optimizer::StochasticProgramOptimizer, optimizer::AbstractStructuredOptimizer)
    try
        return MOI.get(optimizer, MasterOptimizer())
    catch err
        return nothing
    end
end

function subproblem_optimizer(sp_optimizer::StochasticProgramOptimizer, optimizer::AbstractStructuredOptimizer)
    try
        return MOI.get(optimizer, SubProblemOptimizer())
    catch err
        return nothing
    end
end

function reset_optimizer!(sp_optimizer::StochasticProgramOptimizer{Opt}, optimizer) where Opt <: StochasticProgramOptimizerType
    @warn "Only optimizers of type $Opt can be set. Consider reinstantiating stochastic program."
    return nothing
end

function reset_optimizer!(sp_optimizer::StochasticProgramOptimizer, optimizer::MOI.AbstractOptimizer)
    MOIU.reset_optimizer(sp_optimizer.optimizer, optimizer)
    return nothing
end

function reset_optimizer!(sp_optimizer::StochasticProgramOptimizer, optimizer::AbstractStructuredOptimizer)
    if sp_optimizer.optimizer != nothing && sp_optimizer.optimizer isa AbstractStructuredOptimizer
        restore_structure!(sp_optimizer.optimizer)
    end
    sp_optimizer.optimizer = optimizer
    return nothing
end

function reset_optimizer!(sp_optimizer::StochasticProgramOptimizer, optimizer::AbstractSampledOptimizer)
    sp_optimizer.optimizer = optimizer
    return nothing
end

function set_optimizer!(sp_optimizer::StochasticProgramOptimizer, optimizer_constructor)
    optimizer = MOI.instantiate(optimizer_constructor; with_bridge_type = bridge_type(optimizer_constructor))
    reset_optimizer!(sp_optimizer, optimizer)
    sp_optimizer.optimizer_constructor = optimizer_constructor
    return nothing
end

function restore_structure!(sp_optimizer::StochasticProgramOptimizer)
    restore_structure!(sp_optimizer, sp_optimizer.optimizer)
    return nothing
end

function restore_structure!(sp_optimizer::StochasticProgramOptimizer, ::MOI.AbstractOptimizer)
    # Nothing to do
    return nothing
end

function restore_structure!(sp_optimizer::StochasticProgramOptimizer, optimizer::AbstractStructuredOptimizer)
    restore_structure!(optimizer)
    return nothing
end


text:
# QualType
QualType(ptr::CXType_, quals::Unsigned) = QualType(clang_QualType_constructFromTypePtr(ptr), quals)

getCanonicalTypeInternal(ptr::CXType_) = QualType(clang_Type_getCanonicalTypeInternal(ptr))
getTypePtr(x::AbstractQualType)::CXType_ = clang_QualType_getTypePtr(x)
getTypePtrOrNull(x::AbstractQualType)::CXType_ = clang_QualType_getTypePtrOrNull(x)

isCanonical(x::AbstractQualType) = clang_QualType_isCanonical(x)

isNull(x::AbstractQualType) = clang_QualType_isNull(x)

isConstQualified(x::AbstractQualType) = clang_QualType_isConstQualified(x)
isRestrictQualified(x::AbstractQualType) = clang_QualType_isRestrictQualified(x)
isVolatileQualified(x::AbstractQualType) = clang_QualType_isVolatileQualified(x)

hasQualifiers(x::AbstractQualType) = clang_QualType_hasQualifiers(x)

isLocalConstQualified(x::AbstractQualType) = clang_QualType_isLocalConstQualified(x)
isLocalRestrictQualified(x::AbstractQualType) = clang_QualType_isLocalRestrictQualified(x)
isLocalVolatileQualified(x::AbstractQualType) = clang_QualType_isLocalVolatileQualified(x)

hasLocalQualifiers(x::AbstractQualType) = clang_QualType_hasLocalQualifiers(x)

withConst(x::AbstractQualType) = QualType(clang_QualType_withConst(x))
withRestrict(x::AbstractQualType) = QualType(clang_QualType_withRestrict(x))
withVolatile(x::AbstractQualType) = QualType(clang_QualType_withVolatile(x))

addConst(x::AbstractQualType) = QualType(clang_QualType_addConst(x))
addRestrict(x::AbstractQualType) = QualType(clang_QualType_addRestrict(x))
addVolatile(x::AbstractQualType) = QualType(clang_QualType_addVolatile(x))

getCanonicalType(x::AbstractQualType) = QualType(clang_QualType_getCanonicalType(x))
getLocalUnqualifiedType(x::AbstractQualType) = QualType(clang_QualType_getLocalUnqualifiedType(x))
getUnqualifiedType(x::AbstractQualType) = QualType(clang_QualType_getUnqualifiedType(x))

dump(x::AbstractQualType) = clang_QualType_dump(x)

function getAsString(x::AbstractQualType)
    str = clang_QualType_getAsString(x)
    s = unsafe_string(str)
    clang_QualType_disposeString(str)
    return s
end

getPointeeType(ty_ptr::CXType_) = clang_Type_getPointeeType(ty_ptr)
getPointeeType(ty::AbstractQualType) = QualType(getPointeeType(getTypePtr(ty)))

isVoidType(ty_ptr::CXType_) = clang_Type_isVoidType(ty_ptr)
isVoidType(ty::AbstractQualType) = isVoidType(getTypePtr(ty))

isBooleanType(ty_ptr::CXType_) = clang_Type_isBooleanType(ty_ptr)
isBooleanType(ty::AbstractQualType) = isBooleanType(getTypePtr(ty))

isPointerType(ty_ptr::CXType_) = clang_Type_isPointerType(ty_ptr)
isPointerType(ty::AbstractQualType) = isPointerType(getTypePtr(ty))

isFunctionPointerType(ty_ptr::CXType_) = clang_Type_isFunctionPointerType(ty_ptr)
isFunctionPointerType(ty::AbstractQualType) = isFunctionPointerType(getTypePtr(ty))

isFunctionType(ty_ptr::CXType_) = clang_Type_isFunctionType(ty_ptr)
isFunctionType(ty::AbstractQualType) = isFunctionType(getTypePtr(ty))

isMemberFunctionPointerType(ty_ptr::CXType_) = clang_Type_isMemberFunctionPointerType(ty_ptr)
isMemberFunctionPointerType(ty::AbstractQualType) = isMemberFunctionPointerType(getTypePtr(ty))

isReferenceType(ty_ptr::CXType_) = clang_Type_isReferenceType(ty_ptr)
isReferenceType(ty::AbstractQualType) = isReferenceType(getTypePtr(ty))

isCharType(ty_ptr::CXType_) = clang_Type_isCharType(ty_ptr)
isCharType(ty::AbstractQualType) = isCharType(getTypePtr(ty))

isEnumeralType(ty_ptr::CXType_) = clang_Type_isEnumeralType(ty_ptr)
isEnumeralType(ty::AbstractQualType) = isEnumeralType(getTypePtr(ty))

isArrayType(ty_ptr::CXType_) = clang_Type_isArrayType(ty_ptr)
isArrayType(ty::AbstractQualType) = isArrayType(getTypePtr(ty))

# BuiltinTypes
isa_BuiltinType(ty_ptr::CXType_) = clang_isa_BuiltinType(ty_ptr)
isa_BuiltinType(ty::AbstractQualType) = isa_BuiltinType(getTypePtr(ty))

isa_BuiltinType_Void(ty_ptr::CXType_) = clang_isa_BuiltinType_Void(ty_ptr)
isa_BuiltinType_Void(ty::AbstractQualType) = isa_BuiltinType_Void(getTypePtr(ty))

isa_BuiltinType_Bool(ty_ptr::CXType_) = clang_isa_BuiltinType_Bool(ty_ptr)
isa_BuiltinType_Bool(ty::AbstractQualType) = isa_BuiltinType_Bool(getTypePtr(ty))

isa_BuiltinType_Char_U(ty_ptr::CXType_) = clang_isa_BuiltinType_Char_U(ty_ptr)
isa_BuiltinType_Char_U(ty::AbstractQualType) = isa_BuiltinType_Char_U(getTypePtr(ty))

isa_BuiltinType_Char_S(ty_ptr::CXType_) = clang_isa_BuiltinType_Char_S(ty_ptr)
isa_BuiltinType_Char_S(ty::AbstractQualType) = isa_BuiltinType_Char_S(getTypePtr(ty))

isa_BuiltinType_WChar_U(ty_ptr::CXType_) = clang_isa_BuiltinType_WChar_U(ty_ptr)
isa_BuiltinType_WChar_U(ty::AbstractQualType) = isa_BuiltinType_WChar_U(getTypePtr(ty))

isa_BuiltinType_WChar_S(ty_ptr::CXType_) = clang_isa_BuiltinType_WChar_S(ty_ptr)
isa_BuiltinType_WChar_S(ty::AbstractQualType) = isa_BuiltinType_WChar_S(getTypePtr(ty))

isa_BuiltinType_Char8(ty_ptr::CXType_) = clang_isa_BuiltinType_Char8(ty_ptr)
isa_BuiltinType_Char8(ty::AbstractQualType) = isa_BuiltinType_Char8(getTypePtr(ty))

isa_BuiltinType_Char16(ty_ptr::CXType_) = clang_isa_BuiltinType_Char16(ty_ptr)
isa_BuiltinType_Char16(ty::AbstractQualType) = isa_BuiltinType_Char16(getTypePtr(ty))

isa_BuiltinType_Char32(ty_ptr::CXType_) = clang_isa_BuiltinType_Char32(ty_ptr)
isa_BuiltinType_Char32(ty::AbstractQualType) = isa_BuiltinType_Char32(getTypePtr(ty))

isa_BuiltinType_SChar(ty_ptr::CXType_) = clang_isa_BuiltinType_SChar(ty_ptr)
isa_BuiltinType_SChar(ty::AbstractQualType) = isa_BuiltinType_SChar(getTypePtr(ty))

isa_BuiltinType_Short(ty_ptr::CXType_) = clang_isa_BuiltinType_Short(ty_ptr)
isa_BuiltinType_Short(ty::AbstractQualType) = isa_BuiltinType_Short(getTypePtr(ty))

isa_BuiltinType_Int(ty_ptr::CXType_) = clang_isa_BuiltinType_Int(ty_ptr)
isa_BuiltinType_Int(ty::AbstractQualType) = isa_BuiltinType_Int(getTypePtr(ty))

isa_BuiltinType_Long(ty_ptr::CXType_) = clang_isa_BuiltinType_Long(ty_ptr)
isa_BuiltinType_Long(ty::AbstractQualType) = isa_BuiltinType_Long(getTypePtr(ty))

isa_BuiltinType_LongLong(ty_ptr::CXType_) = clang_isa_BuiltinType_LongLong(ty_ptr)
isa_BuiltinType_LongLong(ty::AbstractQualType) = isa_BuiltinType_LongLong(getTypePtr(ty))

isa_BuiltinType_Int128(ty_ptr::CXType_) = clang_isa_BuiltinType_Int128(ty_ptr)
isa_BuiltinType_Int128(ty::AbstractQualType) = isa_BuiltinType_Int128(getTypePtr(ty))

isa_BuiltinType_UChar(ty_ptr::CXType_) = clang_isa_BuiltinType_UChar(ty_ptr)
isa_BuiltinType_UChar(ty::AbstractQualType) = isa_BuiltinType_UChar(getTypePtr(ty))

isa_BuiltinType_UShort(ty_ptr::CXType_) = clang_isa_BuiltinType_UShort(ty_ptr)
isa_BuiltinType_UShort(ty::AbstractQualType) = isa_BuiltinType_UShort(getTypePtr(ty))

isa_BuiltinType_UInt(ty_ptr::CXType_) = clang_isa_BuiltinType_UInt(ty_ptr)
isa_BuiltinType_UInt(ty::AbstractQualType) = isa_BuiltinType_UInt(getTypePtr(ty))

isa_BuiltinType_ULong(ty_ptr::CXType_) = clang_isa_BuiltinType_ULong(ty_ptr)
isa_BuiltinType_ULong(ty::AbstractQualType) = isa_BuiltinType_ULong(getTypePtr(ty))

isa_BuiltinType_ULongLong(ty_ptr::CXType_) = clang_isa_BuiltinType_ULongLong(ty_ptr)
isa_BuiltinType_ULongLong(ty::AbstractQualType) = isa_BuiltinType_ULongLong(getTypePtr(ty))

isa_BuiltinType_UInt128(ty_ptr::CXType_) = clang_isa_BuiltinType_UInt128(ty_ptr)
isa_BuiltinType_UInt128(ty::AbstractQualType) = isa_BuiltinType_UInt128(getTypePtr(ty))

isa_BuiltinType_Float(ty_ptr::CXType_) = clang_isa_BuiltinType_Float(ty_ptr)
isa_BuiltinType_Float(ty::AbstractQualType) = isa_BuiltinType_Float(getTypePtr(ty))

isa_BuiltinType_Double(ty_ptr::CXType_) = clang_isa_BuiltinType_Double(ty_ptr)
isa_BuiltinType_Double(ty::AbstractQualType) = isa_BuiltinType_Double(getTypePtr(ty))

isa_BuiltinType_LongDouble(ty_ptr::CXType_) = clang_isa_BuiltinType_LongDouble(ty_ptr)
isa_BuiltinType_LongDouble(ty::AbstractQualType) = isa_BuiltinType_LongDouble(getTypePtr(ty))

isa_BuiltinType_Float128(ty_ptr::CXType_) = clang_isa_BuiltinType_Float128(ty_ptr)
isa_BuiltinType_Float128(ty::AbstractQualType) = isa_BuiltinType_Float128(getTypePtr(ty))

isa_BuiltinType_Half(ty_ptr::CXType_) = clang_isa_BuiltinType_Half(ty_ptr)
isa_BuiltinType_Half(ty::AbstractQualType) = isa_BuiltinType_Half(getTypePtr(ty))

isa_BuiltinType_BFloat16(ty_ptr::CXType_) = clang_isa_BuiltinType_BFloat16(ty_ptr)
isa_BuiltinType_BFloat16(ty::AbstractQualType) = isa_BuiltinType_BFloat16(getTypePtr(ty))

isa_BuiltinType_Float16(ty_ptr::CXType_) = clang_isa_BuiltinType_Float16(ty_ptr)
isa_BuiltinType_Float16(ty::AbstractQualType) = isa_BuiltinType_Float16(getTypePtr(ty))

isa_BuiltinType_NullPtr(ty_ptr::CXType_) = clang_isa_BuiltinType_NullPtr(ty_ptr)
isa_BuiltinType_NullPtr(ty::AbstractQualType) = isa_BuiltinType_NullPtr(getTypePtr(ty))

# ComplexType
isa_ComplexType(ty_ptr::CXType_) = clang_isa_ComplexType(ty_ptr)
isa_ComplexType(ty::AbstractQualType) = isa_ComplexType(getTypePtr(ty))

# PointerType
getPointeeType(x::PointerType) = QualType(clang_PointerType_getPointeeType(getTypePtr(x)))

# ReferenceType
getPointeeType(x::ReferenceType) = QualType(clang_ReferenceType_getPointeeType(getTypePtr(x)))

# LValueReferenceType
isa_LValueReferenceType(ty_ptr::CXType_) = clang_isa_LValueReferenceType(ty_ptr)
isa_LValueReferenceType(ty::AbstractQualType) = isa_LValueReferenceType(getTypePtr(ty))

# RValueReferenceType
isa_RValueReferenceType(ty_ptr::CXType_) = clang_isa_RValueReferenceType(ty_ptr)
isa_RValueReferenceType(ty::AbstractQualType) = isa_RValueReferenceType(getTypePtr(ty))

# MemberPointerType
isa_MemberPointerType(ty_ptr::CXType_) = clang_isa_MemberPointerType(ty_ptr)
isa_MemberPointerType(ty::AbstractQualType) = isa_MemberPointerType(getTypePtr(ty))

getPointeeType(x::MemberPointerType) = QualType(clang_MemberPointerType_getPointeeType(getTypePtr(x)))

getClass(x::MemberPointerType)::CXType_ = clang_MemberPointerType_getClass(getTypePtr(x))

# ConstantArrayType
isa_ConstantArrayType(ty_ptr::CXType_) = clang_isa_ConstantArrayType(ty_ptr)
isa_ConstantArrayType(ty::AbstractQualType) = isa_ConstantArrayType(getTypePtr(ty))

# IncompleteArrayType
isa_IncompleteArrayType(ty_ptr::CXType_) = clang_isa_IncompleteArrayType(ty_ptr)
isa_IncompleteArrayType(ty::AbstractQualType) = isa_IncompleteArrayType(getTypePtr(ty))

# VariableArrayType
isa_VariableArrayType(ty_ptr::CXType_) = clang_isa_VariableArrayType(ty_ptr)
isa_VariableArrayType(ty::AbstractQualType) = isa_VariableArrayType(getTypePtr(ty))

# DependentSizedArrayType
isa_DependentSizedArrayType(ty_ptr::CXType_) = clang_isa_DependentSizedArrayType(ty_ptr)
isa_DependentSizedArrayType(ty::AbstractQualType) = isa_DependentSizedArrayType(getTypePtr(ty))

# FunctionType
getReturnType(x::AbstractFunctionType) = QualType(clang_FunctionType_getReturnType(getTypePtr(x)))

# FunctionNoProtoType
isa_FunctionNoProtoType(ty_ptr::CXType_) = clang_isa_FunctionNoProtoType(ty_ptr)
isa_FunctionNoProtoType(ty::AbstractQualType) = isa_FunctionNoProtoType(getTypePtr(ty))

# FunctionProtoType
isa_FunctionProtoType(ty_ptr::CXType_) = clang_isa_FunctionProtoType(ty_ptr)
isa_FunctionProtoType(ty::AbstractQualType) = isa_FunctionProtoType(getTypePtr(ty))

getNumParams(x::FunctionProtoType) = clang_FunctionProtoType_getNumParams(getTypePtr(x))
getParamType(x::FunctionProtoType, i::Integer) = QualType(clang_FunctionProtoType_getParamType(getTypePtr(x), i))

# TypedefType
isa_TypedefType(ty_ptr::CXType_) = clang_isa_TypedefType(ty_ptr)
isa_TypedefType(ty::AbstractQualType) = isa_TypedefType(getTypePtr(ty))

desugar(x::TypedefType) = QualType(clang_TypedefType_desugar(getTypePtr(x)))

# TagType
isa_TagType(ty_ptr::CXType_) = clang_isa_TagType(ty_ptr)
isa_TagType(ty::AbstractQualType) = isa_TagType(getTypePtr(ty))

# RecordType
isa_RecordType(ty_ptr::CXType_) = clang_isa_RecordType(ty_ptr)
isa_RecordType(ty::AbstractQualType) = isa_RecordType(getTypePtr(ty))

# EnumType
isa_EnumType(ty_ptr::CXType_) = clang_isa_EnumType(ty_ptr)
isa_EnumType(ty::AbstractQualType) = isa_EnumType(getTypePtr(ty))

getDecl(x::EnumType) = EnumDecl(clang_EnumType_getDecl(x))
getIntegerType(x::EnumType) = getIntegerType(getDecl(x))
getName(x::EnumType) = getName(getDecl(x))

# TemplateTypeParmType
isa_TemplateTypeParmType(ty_ptr::CXType_) = clang_isa_TemplateTypeParmType(ty_ptr)
isa_TemplateTypeParmType(ty::AbstractQualType) = isa_TemplateTypeParmType(getTypePtr(ty))

# SubstTemplateTypeParmType
isa_SubstTemplateTypeParmType(ty_ptr::CXType_) = clang_isa_SubstTemplateTypeParmType(ty_ptr)
isa_SubstTemplateTypeParmType(ty::AbstractQualType) = isa_SubstTemplateTypeParmType(getTypePtr(ty))

# SubstTemplateTypeParmPackType
isa_SubstTemplateTypeParmPackType(ty_ptr::CXType_) = clang_isa_SubstTemplateTypeParmPackType(ty_ptr)
isa_SubstTemplateTypeParmPackType(ty::AbstractQualType) = isa_SubstTemplateTypeParmPackType(getTypePtr(ty))

# TemplateSpecializationType
isa_TemplateSpecializationType(ty_ptr::CXType_) = clang_isa_TemplateSpecializationType(ty_ptr)
isa_TemplateSpecializationType(ty::AbstractQualType) = isa_TemplateSpecializationType(getTypePtr(ty))

function isCurrentInstantiation(x::TemplateSpecializationType)
    return clang_TemplateSpecializationType_isCurrentInstantiation(getTypePtr(ty))
end

function isTypeAlias(x::TemplateSpecializationType)
    return clang_TemplateSpecializationType_isTypeAlias(getTypePtr(ty))
end

function getAliasedType(x::TemplateSpecializationType)
    return QualType(clang_TemplateSpecializationType_getAliasedType(getTypePtr(ty)))
end

function getNumArgs(x::TemplateSpecializationType)
    return clang_TemplateSpecializationType_getNumArgs(getTypePtr(ty))
end

function getArg(x::TemplateSpecializationType, i::Integer)
    @check_ptrs x
    return TemplateArgument(clang_TemplateSpecializationType_getArg(x, i))
end

function isSugared(x::TemplateSpecializationType)
    return clang_TemplateSpecializationType_isSugared(getTypePtr(ty))
end

function desugar(x::TemplateSpecializationType)
    return QualType(clang_TemplateSpecializationType_desugar(getTypePtr(ty)))
end

getName(x::TemplateSpecializationType) = getName(getAsTemplateDecl(getTemplateName(x)))

# ElaboratedType
isa_ElaboratedType(ty_ptr::CXType_) = clang_isa_ElaboratedType(ty_ptr)
isa_ElaboratedType(ty::AbstractQualType) = isa_ElaboratedType(getTypePtr(ty))

desugar(x::ElaboratedType) = clang_ElaboratedType_desugar(getTypePtr(x))

# DependentNameType
isa_DependentNameType(ty_ptr::CXType_) = clang_isa_DependentNameType(ty_ptr)
isa_DependentNameType(ty::AbstractQualType) = isa_DependentNameType(getTypePtr(ty))

# DependentTemplateSpecializationType
isa_DependentTemplateSpecializationType(ty_ptr::CXType_) = clang_isa_DependentTemplateSpecializationType(ty_ptr)
isa_DependentTemplateSpecializationType(ty::AbstractQualType) = isa_DependentTemplateSpecializationType(getTypePtr(ty))


text:
using Cardio
using Documenter

makedocs(;
    modules=[Cardio],
    authors="Marius Pille <marius.pille@hu-berlin.de> and contributors",
    repo="https://github.com/mapi1/Cardio.jl/blob/{commit}{path}#L{line}",
    sitename="Cardio.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mapi1.github.io/Cardio.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Detection" => "Detection.md",
        "BRS" => "BRS.md",
        "HRV" => "HRV.md",
        "Utility" => "Utility.md",
    ],
)

deploydocs(;
    repo="github.com/mapi1/Cardio.jl",
)


text:
# ---
# title: 1003. Check If Word Is Valid After Substitutions
# id: problem1003
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/check-if-word-is-valid-after-substitutions/description/>
# hidden: true
# ---
#
# Given a string `s`, determine if it is **valid**.
#
# A string `s` is **valid** if, starting with an empty string `t = ""`, you can
# **transform**`t` **into**`s` after performing the following operation **any
# number of times** :
#
#   * Insert string `"abc"` into any position in `t`. More formally, `t` becomes `tleft + "abc" + tright`, where `t == tleft + tright`. Note that `tleft` and `tright` may be **empty**.
#
# Return `true` _if_`s` _is a **valid** string, otherwise, return_ `false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aabcbc"
#     Output: true
#     Explanation:
#     "" -> " _abc_ " -> "a _abc_ bc"
#     Thus, "aabcbc" is valid.
#
# **Example 2:**
#
#
#
#     Input: s = "abcabcababcc"
#     Output: true
#     Explanation:
#     "" -> " _abc_ " -> "abc _abc_ " -> "abcabc _abc_ " -> "abcabcab _abc_ c"
#     Thus, "abcabcababcc" is valid.
#
#
# **Example 3:**
#
#
#
#     Input: s = "abccba"
#     Output: false
#     Explanation: It is impossible to get "abccba" using the operation.
#
#
# **Example 4:**
#
#
#
#     Input: s = "cababc"
#     Output: false
#     Explanation: It is impossible to get "cababc" using the operation.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 2 * 104`
#   * `s` consists of letters `'a'`, `'b'`, and `'c'`
#
#
## @lc code=start
using LeetCode

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


text:
# TODO: Adapt python model below to new Julia format
# ################################################################################
# ##                         Gallium Arsenide (GaAs)                            ##
# ################################################################################
#
# def n_GaAs_sym():
#     """This function creates a symbolic representation (using SymPy) of the
#     a model for the temperature and wavelength dependence of GaAs's refractive
#     index. The equation form and fit parameters are based on "Improved dispersion
#     relations for GaAs and applications to nonlinear optics" by Skauli et al.,
#     JAP 94, 6447 (2003); doi: 10.1063/1.1621740
#     This model is then exported to other functions that use it and its
#     derivatives to return index, group index and GVD values as a function
#     of temperature and wavelength.
#     Variable units are lm in [um] and T in [deg C]
#     """
#
#     lm, T = sp.symbols('lm T')
#     T0 = 22 # reference temperature in [Deg C]
#     deltaT = T-T0
#     A = 0.689578
#     eps2 = 12.99386
#     G3 = 2.18176e-3
#     E0 = 1.425 - 3.7164e-4 * deltaT - 7.497e-7* deltaT**2
#     E1 = 2.400356 - 5.1458e-4 * deltaT
#     E2 = 7.691979 - 4.6545e-4 * deltaT
#     E3 = 3.4303e-2 + 1.136e-5 * deltaT
#     E_phot = (u.h * u.c).to(u.eV*u.um).magnitude / lm  #
#     n_sym = sp.sqrt( 1  + (A/np.pi)*sp.log((E1**2 - E_phot**2) / (E0**2-E_phot**2)) \
#                         + (eps2/np.pi)*sp.log((E2**2 - E_phot**2) / (E1**2-E_phot**2)) \
#                         + G3 / (E3**2-E_phot**2) )
#
#     return lm, T, n_sym
#
# ###### Temperature Dependent Sellmeier Equation for phase-matching calculations
# def n_GaAs(lm_in,T_in):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the index of refraction of 5% MgO:LiNbO3. Equation form is based on
#     "Temperature and wavelength dependent refractive index equations
#     for MgO-doped congruent and stoichiometric LiNbO3"
#     by Gayer et al., Applied Physics B 91, p.343-348 (2008)
#     Variable units are lm in [um] and T in [deg C]
#     """
# #    lm_um = Q_(lm_in).to(u.um).magnitude
# #    T_C = Q_(T_in).to(u.degC).magnitude
# #    lm, T, n_sym = n_GaAs_sym()
# #    n = sp.lambdify([lm,T],n_sym,'numpy')
# #    return n(lm_um, T_C)
#
#     lm_um = np.array([Q_(lm_in).to(u.um).magnitude]).flatten()
#     T_C = np.array([Q_(T_in).to(u.degC).magnitude]).flatten()
#     lm, T, n_sym = n_GaAs_sym()
#     n = sp.lambdify([lm,T],n_sym,'numpy')
#     output = np.zeros((T_C.size,lm_um.size))
#     for T_idx, TT in enumerate(T_C):
#         output[T_idx,:] = n(lm_um, T_C[T_idx])
#     return output
#
#
# def n_g_GaAs(lm_in,T_in):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group index of refraction of 5% MgO:LiNbO3. Equation form is based on
#     "Temperature and wavelength dependent refractive index equations
#     for MgO-doped congruent and stoichiometric LiNbO3"
#     by Gayer et al., Applied Physics B 91, p.343-348 (2008)
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_GaAs_sym()
#     n_sym_prime = sp.diff(n_sym,lm)
#     n_g_sym = n_sym - lm*n_sym_prime
#     n_g = sp.lambdify([lm,T],n_g_sym,'numpy')
#     return n_g(lm_um, T_C)
#
#
# def gvd_GaAs(lm_in,T_in):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group velocity dispersion of 5% MgO:LiNbO3. Equation form is based on
#     "Temperature and wavelength dependent refractive index equations
#     for MgO-doped congruent and stoichiometric LiNbO3"
#     by Gayer et al., Applied Physics B 91, p.343-348 (2008)
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_GaAs_sym()
#     n_sym_double_prime = sp.diff(n_sym,lm,lm)
#     c = Q_(3e8,'m/s') # unitful definition of speed of light
#     gvd_sym_no_prefactors = (lm**3)*n_sym_double_prime # symbolic gvd without unitful prefactors, to be made unitful below
#     gvd_no_prefactors = sp.lambdify([lm,T],gvd_sym_no_prefactors,'numpy') # numerical gvd without unitful prefactors
#     gvd = (1 / (2 * np.pi * (c**2))) * Q_(gvd_no_prefactors(lm_um,T_C),'um')
#     return gvd.to('fs**2 / mm')


text:
# # Granger Causality Estimation
#
#md # The notebook can be viewed here:
#md # * [![binder](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/granger_examples.ipynb)
#md # * [![nbviewer](https://img.shields.io/badge/render-nbviewer-orange.svg)](@__NBVIEWER_ROOT_URL__/generated/granger_examples.ipynb)
#
#nb # ### Acknowledgement
#nb # This work was funded by the German Federal Ministry of Education and Research ([BMBF](https://www.bmbf.de/)) in the project ALICE III under grant ref. 01IS18049B.
#
# ### Load packages

using TimeSeriesCausality
using Distributions: MvNormal
using Plots: plot
using Printf

# ### Data
# We will generate (simulate) 1e6 two-channel samples using a design (evolution) matrix A of order 3. where:
# - channel 1 is the causal (independent) time series
# - channel 2 is the effect (dependent) time series

## design (evolution) matrix
##            Channel 1->    Channel 2->
##           t-1  t-2  t-3  t-1  t-2  t-3
designer = [0.4 -0.6  0.8  0.0  0.0  0.0;  # -> Channel 1
            0.5  0.9  0.0  0.0  0.0  0.7]  # -> Channel 2

## data generation
n_samples = 128*1024  # number of time steps
segment_length = 1024  # segment length
noise_cov = MvNormal([0.25 0.0; 0.0 0.64])  # uncorrelated Cov matrix of noise
noise = rand(noise_cov, n_samples)'  # sampling from the noise distribution

## pre-allocation and initial values
signal = zeros(n_samples, 2)
signal[1:3, :] = rand(3, 2) + noise[1:3, :]

## simulation
for t in 4:n_samples
    signal[t, :] = designer * reshape(signal[t-3:t-1, :], :, 1) + noise[t, :]
end

plot(1:segment_length, [signal[1:segment_length, 1], signal[1:segment_length, 2]],
     title = "Signals",
     label = ["Channel 1" "Channel 2"],
     xlabel = "Time steps",
     ylabel = "",
     lw = 2)

# ### Granger method of causal estimation
# Granger method can estimate the causal index for a 2-channel signal. The larger the value, the stronger the causal dependence. A positive value indicates the flow of information from channel 1 to 2, and vice versa.
# We have also included JackKnife method for calculate the standard deviation of error for the estimation.

granger_idx, err_std = granger_est(signal, segment_length; order=3, method="jackknife")
@printf "Granger causality index is %.3f with std error of %.3f" granger_idx err_std

# ### Order as a hyperparameter
# Since the order (time step delay) in which the signals interact is unknown by us, we have implemented Akaike and Bayesian information criteria to estimate the most informative order.

## the range of orders to look into
order_range = 1:7

## Akaike information criterion
aic = granger_aic(signal, segment_length, order_range)

## Bayesian information criterion
bic = granger_bic(signal, segment_length, order_range)

plot(order_range, [aic, bic],
     title = "Akaike vs Bayesian information criteria",
     label = ["Akaike IC" "Bayesian IC"],
     xlabel = "Order",
     ylabel = "IC",
     lw = 2)


text:
#=
Modified from https://github.com/ryarazi/ExactRiemannProblemSolver, Roy Arazi
MIT License

Copyright (c) 2018 Roy Arazi

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

# Modified data: 2021/8/21
#Type Define:
#struct for saving the state of the system in some half
struct HydroStatus
	rho::Float64 #density
	u::Float64 #velocity
	p::Float64 #pressure
	gamma::Float64 #adiabatic index
	c::Float64 #speed of cound
	A::Float64 #Toro constant (4.8)
	B::Float64 #Toro constant (4.8)

	function HydroStatus(rho, u, p, gamma)
		c = sqrt(gamma * p / rho)
		A = 2. / (gamma + 1.) / rho
		B = (gamma - 1.) / (gamma + 1.) * p
		new(rho, u, p, gamma, c, A, B)
	end
end

#like the HydroStatus but saves only the important parameters - density, pressure and velocity
struct HydroStatusSimple
	rho::Float64 #density
	u::Float64 #velocity
	p::Float64 #pressure
end
HydroStatusSimple(status::HydroStatus) = HydroStatusSimple(status.rho, status.u, status.p)

#different types to dispatch the first guess of p_star calculation
abstract type FirstGuessScheme end
struct TwoRarefaction <: FirstGuessScheme end
struct PrimitiveValue <: FirstGuessScheme end
struct TwoShock <: FirstGuessScheme end
struct MeanPressure <: FirstGuessScheme end

#the type of the wave that will be created on some side of the contact discontinuity
abstract type WaveType end
struct Rarefaction <: WaveType end
struct Shock <: WaveType end
struct Vaccum <: WaveType end #for vaccum generated in the middle

#the side in which the calculation is happening
abstract type Side end
struct LeftSide <: Side end
struct RightSide <: Side end

#in some equations there is factor difference between the left and right side calculations
side_factor(::LeftSide) = -1
side_factor(::RightSide) = 1

#choose the some variable based on the side given
choose_side(left, right, ::LeftSide) = left
choose_side(left, right, ::RightSide) = right

# Initial guess of the p_star in iteration

# (1) Two-Rarefaction approximation, see (4.46) from Toro
function p_guess(left::HydroStatus, right::HydroStatus, ::TwoRarefaction)
		@assert left.gamma == right.gamma
		gamma = left.gamma
		gamma_power = (gamma-1) / 2. / gamma
		delta_u = right.u - left.u
		return ( (left.c + right.c - 0.5 * (gamma - 1) * delta_u) /
			   (left.c/left.p^gamma_power + right.c/right.p^gamma_power) ) ^
		       (1. / gamma_power)
end
# (2) Linearised solution based on primitive variables, see (4.47) from Toro
function p_guess(left::HydroStatus, right::HydroStatus, ::PrimitiveValue)
	delta_u = right.u - left.u
	mean_pressure = 0.5 * (left.p + right.p)
	mean_density = 0.5 * (left.rho + right.rho)
	mean_speed_of_sound = 0.5 * (left.c + right.c)
	return mean_pressure - 0.5 * delta_u * mean_density * mean_speed_of_sound
end
# (3) Two-Shock approximation, see (4.48) from Toro
function p_guess(left::HydroStatus, right::HydroStatus, ::TwoShock)
	p_hat = p_guess(left, right, PrimitiveValue())
	g_left = sqrt(left.A / (p_hat + left.B))
	g_right = sqrt(right.A / (p_hat + right.B))
	delta_u = right.u - left.u
	return (g_left * left.p + g_right * right.p - delta_u)/(g_left+g_right)
end
# (4) Arithmetic mean, see (4.49) from Toro
p_guess(left::HydroStatus, right::HydroStatus, ::MeanPressure) = 0.5*(left.p + right.p)

# (5) to provide a guess value for pressure p_star in the Star Region.
# The choice is made according to adaptive Riemann solver using
# the PrimitiveValue, TwoRarefaction and TwoShock approximate Riemann solvers.
# See Sect. 9.5 of Chapt. 9 of Ref. 1
function p_guess(left::HydroStatus, right::HydroStatus)
    Quser=2.0
    PPV=max(0.0, p_guess(left,right,PrimitiveValue()))
    Pmin=min(left.p, right.p)
    Pmax=max(left.p, right.p)
    Qmax=Pmax/Pmin
    if (Qmax≤Quser)&&(Pmin≤PPV≤Pmax)
        return PPV
    else
        if PPV<Pmin # Select Two-Rarefaction Riemann solver
            return p_guess(left,right,TwoRarefaction())
        else # Select Two-Shock Riemann solver with PVRS as estimate
            return p_guess(left,right,TwoShock())
        end
    end
end

# Riemann Solver:

# function (4.6) and (4.7) from Toro
f(p, status::HydroStatus, ::Shock) = (p - status.p) * sqrt(status.A / (p + status.B))
f(p, status::HydroStatus, ::Rarefaction) = 2. * status.c / (status.gamma - 1.) * ((p / status.p)^((status.gamma - 1.) / 2. /status.gamma) - 1.)
function f(p, status::HydroStatus)
	wave_type = p > status.p ? Shock() : Rarefaction()
	return f(p, status, wave_type)
end
# function (4.37) from Toro
f′(p, status::HydroStatus, ::Shock) = sqrt(status.A/(status.B+p))*(1-(p-status.p)/(2(status.B+p)))
f′(p, status::HydroStatus, ::Rarefaction) = 1/(status.rho*status.c)*(p/status.p)^(-(status.gamma+1)/2/status.gamma)
function f′(p, status::HydroStatus)
	wave_type = p > status.p ? Shock() : Rarefaction()
	return f′(p, status, wave_type)
end
#solver for p_star based on the first guess of the pressure
#uses Newton iteration method to solve the equation
function p_star_calc(left::HydroStatus, right::HydroStatus, TOL::Float64)
	p0 = max(TOL, p_guess(left, right))
	p=p0
	delta_u = right.u - left.u
	MAXIter=100
	I=0
	while true
		p=p0-(f(p0,left)+f(p0,right)+delta_u)/(f′(p0,left)+f′(p0,right))
		if 2.0*abs((p-p0)/(p+p0))≤TOL
			break
		end
		if p<0
			p=TOL
		end
		if I>MAXIter
			error("exceed the max num of iter")
		end
		p0=p
		I+=1
	end
	return p
end
# equation (4.9) from Toro
u_star_calc(p_star, left::HydroStatus, right::HydroStatus) = 0.5 * (left.u + right.u) +
															 0.5 * (f(p_star, right) - f(p_star, left))
# equations (4.53) and (4.60) from Toro
rho_star_calc(p_star, status::HydroStatus, ::Rarefaction) = status.rho * (p_star / status.p) ^ (1. / status.gamma)
#equations (4.50) and (4.57) from Toro
function rho_star_calc(p_star, status::HydroStatus, ::Shock)
	gamma_ratio = (status.gamma - 1.) / (status.gamma + 1.)
	pressure_ratio = p_star / status.p
	return status.rho * (pressure_ratio + gamma_ratio) / (pressure_ratio * gamma_ratio + 1.)
end
#for vaccum generation the density in the middle is zero (= vaccum)
rho_star_calc(p_star, status::HydroStatus, ::Vaccum) = 0.
# equation (4.25) and (4.34) from Toro
c_star_calc(p_star, status::HydroStatus, ::Rarefaction) = status.c * (p_star / status.p) ^ ((status.gamma - 1.) / 2. / status.gamma)
#equation (4.55) and (4.62) from Toro
head_speed_calc(p_star, u_star, status::HydroStatus,
				side::T, ::Rarefaction) where {T<:Side} = status.u + side_factor(side) * status.c

#equation (4.55) and (4.62) from Toro
function tail_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Rarefaction) where {T<:Side}
	c_star = c_star_calc(p_star, status, Rarefaction())
	return u_star + side_factor(side) * c_star
end
#equations (4.52) and (4.59) from Toro
function head_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Shock) where {T<:Side}
	gamma = status.gamma
	gamma_ratio1 = (gamma + 1.) / 2. / gamma
	gamma_ratio2 = (gamma - 1.) / 2. / gamma
	return status.u + side_factor(side) * status.c * sqrt(gamma_ratio1 * p_star / status.p + gamma_ratio2)
end
#for shock the head and the tail are the same
tail_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Shock) where {T<:Side} = head_speed_calc(p_star, u_star, status, side, Shock())

#those are the velocities of the head and the tail of the rarefaciton wave when a vaccum presents
#see (4.76), (4.77), (4.79) and (4.80) from Toro
head_speed_calc(p_star, u_star, status::HydroStatus,
				side::T, ::Vaccum) where {T<:Side} = status.u + side_factor(side) * status.c
tail_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Vaccum) where {T<:Side} = status.u - side_factor(side) * 2 * status.c / (status.gamma - 1.)

#this is the density, velocity and pressure profile of a rarefaction wave in some specific x/t
#see (4.56) and (4.63) from Toro
function rarefaction_profile(x_over_t, status::HydroStatus, side::T) where {T<:Side}
	gamma = status.gamma
	gamma_ratio = (gamma - 1.) / (gamma + 1.)
	gamma_plus = 2. / (gamma + 1.)
	gamma_minus = 2. / (gamma - 1.)

	rarefaction_factor = (gamma_plus - side_factor(side) * gamma_ratio / status.c * (status.u - x_over_t)) ^ gamma_minus

	rho = status.rho * rarefaction_factor
	u = gamma_plus * (-side_factor(side)*status.c + (1. /gamma_minus) * status.u + x_over_t)
	p = status.p * rarefaction_factor ^ gamma

	return HydroStatusSimple(rho, u, p)
end

# this is the MAIN function that samples the Riemann problem and returns the density, pressure and velocity in the space points
#specified by the array x. It gets the following parameters:
# x - an array type which saves all the points in space in which we want to sample the problem (pay attention - the contact discontinuity begins at x=0.)
# t - a floating number which tells us in what time we sample to problem
# left + right - struct from type HydroStatus (see "TypeDefine.jl") to get the initial state of the system
# TOL - the tolerance for the convergance of the p_star finding algorithm
function RiemannSolver(x, t::Float64,left::HydroStatus,right::HydroStatus;TOL::Float64 = 1.e-6)
	@assert left.gamma == right.gamma
	if left.rho == 0 #material is on the right, vaccum on the left
		return sample_riemann_side_vaccum(x, t, right, RightSide())
	elseif right.rho == 0 #material is on the left, vaccum on the right
		return sample_riemann_side_vaccum(x, t, left, LeftSide())
	else
		return sample_riemann_regular(x, t, left, right, TOL)
	end
end
function RiemannSolver(x, t::Float64,left_rho::Float64,left_u::Float64,left_p::Float64,left_gamma::Float64,
									 right_rho::Float64,right_u::Float64,right_p::Float64,right_gamma::Float64;
									 TOL::Float64 = 1.e-6)
	left=HydroStatus(left_rho,left_u,left_p,left_gamma)
	right=HydroStatus(right_rho,right_u,right_p,right_gamma)
	return RiemannSolver(x, t,left,right;TOL=TOL)
end

#sample the riemann problem in the situation where both sides are not vaccum (but vaccum can be generated through the dynamics)
function sample_riemann_regular(x, t::Float64, left::HydroStatus, right::HydroStatus, TOL::Float64)
	@assert left.rho != 0.
	@assert right.rho != 0.

	profile = similar(x, HydroStatusSimple) #we return the values of rho, u and p in each point x in space

	#the status "far" from the center = like the starting condition
	status_left = HydroStatusSimple(left)
	status_right = HydroStatusSimple(right)

	x_over_t = x ./ t

	if 2. / (left.gamma - 1.) * (left.c + right.c) <= right.u - left.u #vaccum generation
		#if vaccum is generated - the density in the middle will be zero
		p_star = 0.
		u_star = 0.
		wave_type_left = Vaccum()
		wave_type_right = Vaccum()
	else
		#if vaccum is not generated we can calculate the star profile in the middle
		p_star = p_star_calc(left, right, TOL)
		u_star = u_star_calc(p_star, left, right)
		#check what kind of waves are in every direction (shock or rarefaction)
		wave_type_left = p_star > left.p ? Shock() : Rarefaction()
		wave_type_right = p_star > right.p ? Shock() : Rarefaction()
	end

	rho_star_left = rho_star_calc(p_star, left, wave_type_left)
	status_left_star = HydroStatusSimple(rho_star_left, u_star, p_star) #the profile in the left near the contact discontinuity
	head_speed_left = head_speed_calc(p_star, u_star, left, LeftSide(), wave_type_left)
	tail_speed_left = tail_speed_calc(p_star, u_star, left, LeftSide(), wave_type_left)

	rho_star_right = rho_star_calc(p_star, right, wave_type_right)
	status_right_star = HydroStatusSimple(rho_star_right, u_star, p_star) #the profile in the right near the contact discontinuity
	head_speed_right = head_speed_calc(p_star, u_star, right, RightSide(), wave_type_right)
	tail_speed_right = tail_speed_calc(p_star, u_star, right, RightSide(), wave_type_right)

	for i = 1:length(x)
		S = x_over_t[i] #this is like the S which Toro use in Section 4.5

		#see Figure 4.14 in Toro for the flow of the following lines
		side = S < u_star ? LeftSide() : RightSide()
		status = choose_side(left, right, side)
		head_speed = choose_side(head_speed_left, head_speed_right, side)
		tail_speed = choose_side(tail_speed_left, tail_speed_right, side)

		#this is used to flip the direction of the inequallity in the branching between the left and the right side
		#the xor which will be in the following lines uses this boolean like some "controlled not"
		#when right_condition is "true" the inequallity will be flipped
		#when right_condition is "false" the inequallity will be stay the same
		right_condition = isa(side, RightSide)

		if xor(S < head_speed, right_condition)
			profile[i] = choose_side(status_left, status_right, side)
		elseif xor(S < tail_speed, right_condition) #can only happen in Rarefaction, because in Shock  head_speed == tail_speed
			profile[i] = rarefaction_profile(S, status, side)
		else
			profile[i] = choose_side(status_left_star, status_right_star, side)
		end
	end

	return profile
end

#sample the riemann problem in the situation where the gas in "side" is actually a vaccum
function sample_riemann_side_vaccum(x, t::Float64, status::HydroStatus, side::T) where {T<:Side}

	profile = similar(x, HydroStatusSimple) #we return the values of rho, u and p in each point x in space
	vaccum = HydroStatusSimple(0., 0., 0.) #vaccum
	simple = HydroStatusSimple(status)

	x_over_t = x ./ t

	tail_speed = tail_speed_calc(0., 0., status, side, Vaccum()) #the head of the rarefaciton wave
	head_speed = head_speed_calc(0., 0., status, side, Vaccum()) #the tail of the rarefaciton wave (=the boundary with the vaccum)

	#see "sample_riemann_regular" for explnation about this boolean
	right_condition = isa(side, RightSide)

	for i = 1:length(x)
		S = x_over_t[i]

		if xor(S < head_speed, right_condition)
			profile[i] = simple
		elseif xor(S < tail_speed, right_condition)
			profile[i] = rarefaction_profile(S, status, side)
		else
			profile[i] = vaccum
		end
	end

	return profile
end
#
function RiemannSolver(x, t::Float64,left_rho::Float64,left_u::Float64,left_v::Float64,left_p::Float64,left_gamma::Float64,
									 right_rho::Float64,right_u::Float64,right_v::Float64,right_p::Float64,right_gamma::Float64;
									 TOL::Float64 = 1.e-6)
	left=HydroStatus(left_rho,left_u,left_p,left_gamma)
	right=HydroStatus(right_rho,right_u,right_p,right_gamma)
	profile=RiemannSolver(x, t,left,right;TOL=TOL)

	# so far can only handel the case no vacuum appear:
	@assert (left.rho>0) && (right.rho>0) && (2. / (left.gamma - 1.) * (left.c + right.c) > right.u - left.u)
	p_star = p_star_calc(left, right, TOL)
	u_star = u_star_calc(p_star, left, right)
	x_over_t = x ./ t
	profile_v=similar(x,Float64)
	for i = 1:length(x)
		S = x_over_t[i] #this is like the S which Toro use in Section 4.5
		if S>u_star
			profile_v[i]=right_v
		else
			profile_v[i]=left_v
		end
	end
	return profile, profile_v
end