text:
const unary_operators = [:(+), :(-), :(!), :(*)]

const numeric_unary_operators = [:(+), :(-)]

const logical_unary_operators = [:(!)]

const elementary_functions = [:abs, :sign, :acos, :acosh, :asin,
                              :asinh, :atan, :atanh, :sin, :sinh,
                              :cos, :cosh, :tan, :tanh, :ceil, :floor,
                              :round, :trunc, :exp, :exp2, :expm1, :log, :log10, :log1p,
                              :log2, :exponent, :sqrt, :gamma, :lgamma, :digamma,
                              :erf, :erfc]

const two_argument_elementary_functions = [:round, :ceil, :floor, :trunc]

const special_comparison_operators = [:isless]

const comparison_operators = [:(==), :(.==), :(!=), :(.!=),
                              :(>), :(.>), :(>=), :(.>=), :(<), :(.<),
                              :(<=), :(.<=)]

const scalar_comparison_operators = [:(==), :(!=), :(>), :(>=),
                                     :(<), :(<=)]

const array_comparison_operators = [:(.==), :(.!=), :(.>), :(.>=), :(.<), :(.<=)]

const vectorized_comparison_operators = [(:(.==), :(==)), (:(.!=), :(!=)),
                                         (:(.>), :(>)), (:(.>=), :(>=)),
                                         (:(.<), :(<)), (:(.<=), :(<=))]

const binary_operators = [:(+), :(.+), :(-), :(.-), :(*), :(.*), :(/), :(./),
                          :(.^), :(div), :(mod), :(fld), :(rem)]

const induced_binary_operators = [:(^)]

const arithmetic_operators = [:(+), :(.+), :(-), :(.-), :(*), :(.*), :(/), :(./),
                              :(.^), :(div), :(mod), :(fld), :(rem)]

const induced_arithmetic_operators = [:(^)]

const biscalar_operators = [:(max), :(min)]

const scalar_arithmetic_operators = [:(+), :(-), :(*), :(/),
                                     :(div), :(mod), :(fld), :(rem)]

const induced_scalar_arithmetic_operators = [:(^)]

const array_arithmetic_operators = [:(+), :(.+), :(-), :(.-), :(.*), :(.^)]

const bit_operators = [:(&), :(|), :($)]

const unary_vector_operators = [:minimum, :maximum, :prod, :sum, :mean, :median, :std,
                                :var, :mad, :norm, :skewness, :kurtosis]

# TODO: dist, iqr, rle, inverse_rle

const pairwise_vector_operators = [:diff, :reldiff, :percent_change]

const cumulative_vector_operators = [:cumprod, :cumsum, :cumsum_kbn, :cummin, :cummax]

const ffts = [:fft]

const binary_vector_operators = [:dot, :cor, :cov, :cor_spearman]

const rowwise_operators = [:rowmins, :rowmaxs, :rowprods, :rowsums,
                           :rowmeans, :rowmedians, :rowstds, :rowvars,
                           :rowffts, :rownorms]

const columnar_operators = [:colmins, :colmaxs, :colprods, :colsums,
                            :colmeans, :colmedians, :colstds, :colvars,
                            :colffts, :colnorms]

const boolean_operators = [:any, :all]

# Swap arguments to fname() anywhere in AST. Returns the number of
# arguments swapped
function swapargs(ast::Expr, fname::Symbol)
    if ast.head == :call &&
       (ast.args[1] == fname ||
        (isa(ast.args[1], Expr) && ast.args[1].head == :curly &&
         ast.args[1].args[1] == fname)) &&
       length(ast.args) == 3

        ast.args[2], ast.args[3] = ast.args[3], ast.args[2]
        1
    else
        n = 0
        for arg in ast.args
            n += swapargs(arg, fname)
        end
        n
    end
end
function swapargs(ast, fname::Symbol)
    ast
    0
end

# Return a block consisting of both the given function and a copy of
# the function in which arguments to the function itself and any
# 2-argument calls to a function of the same name are swapped
macro swappable(func, syms...)
    if (func.head != :function && func.head != :(=)) ||
       func.args[1].head != :call || length(func.args[1].args) != 3
        error("@swappable may only be applied to functions of two arguments")
    end

    func2 = deepcopy(func)
    fname = func2.args[1].args[1]
    if isa(fname, Expr)
        if fname.head == :curly
            fname = fname.args[1]
        else
            error("Unexpected function name $fname")
        end
    end

    for s in unique([fname, syms...])
        if swapargs(func2, s) < 1
            error("No argument swapped")
        end
    end
    esc(Expr(:block, func, func2))
end

#
# Unary operator macros for DataFrames and DataArrays
#

macro dataframe_unary(f)
    esc(:($(f)(d::DataFrame) = DataFrame([$(f)(d[i]) for i=1:size(d, 2)], deepcopy(index(d)))))
end

# Apply unary operator to non-NA members of a DataArray or
# AbstractDataArray
macro dataarray_unary(f, intype, outtype, N...)
    esc(quote
        function $(f){T<:$(intype)}(d::$(isempty(N) ? :(DataArray{T}) : :(DataArray{T,$(N[1])})))
            data = similar(d.data, $(outtype))
            for i = 1:length(data)
                if !d.na[i]
                    data[i] = $(f)(d.data[i])
                end
            end
            DataArray(data, copy(d.na))
        end
        function $(f){T<:$(intype)}(adv::$(isempty(N) ? :(AbstractDataArray{T}) : :(AbstractDataArray{T,$(N[1])})))
            res = similar(adv, $(outtype))
            for i = 1:length(adv)
                res[i] = ($f)(adv[i])
            end
            res
        end
    end)
end

#
# Binary operator macros for DataFrames and DataArrays
#

# Binary operators with one scalar argument
macro dataarray_binary_scalar(vectorfunc, scalarfunc, outtype)
    esc(Expr(:block,
        # DataArray and AbstractDataArray with scalar
        # XXX It would be really nice to make this work with arbitrary
        # types, but doing so results in a bunch of method ambiguity
        # warnings
        {
            quote
                @swappable function $(vectorfunc)(a::DataArray, b::$t)
                    res = DataArray(similar(a.data, $outtype), copy(a.na))
                    for i = 1:length(a)
                        if !res.na[i]
                            res.data[i] = $(scalarfunc)(a.data[i], b)
                        end
                    end
                    res
                end $scalarfunc
                @swappable function $(vectorfunc)(a::AbstractDataArray, b::$t)
                    res = similar(a, $outtype)
                    for i = 1:length(a)
                        res[i] = $(scalarfunc)(a[i], b)
                    end
                    res
                end $scalarfunc
            end
            for t in (:String, :Number)
        }...
    ))
end

# Binary operators with two array arguments
macro dataarray_binary_array(vectorfunc, scalarfunc, outtype)
    esc(Expr(:block,
        # DataArray with other array
        {
            quote
                function $(vectorfunc)(a::$(adata ? :DataArray : :AbstractArray),
                                       b::$(bdata ? :DataArray : :AbstractArray))
                    res = DataArray(Array($outtype, promote_shape(size(a), size(b))), $narule)
                    for i = 1:length(res)
                        if !res.na[i]
                            res.data[i] = $(scalarfunc)($(adata ? :(a.data) : :a)[i],
                                                        $(bdata ? :(b.data) : :b)[i])
                        end
                    end
                    res
                end
            end
            for (adata, bdata, narule) in ((true, true, :(a.na | b.na)),
                                           (true, false, :(copy(a.na))),
                                           (false, true, :(copy(b.na))))
        }...,
        # AbstractDataArray with other array
        # Definitinons with DataArray necessary to avoid ambiguity
        {
            quote
                function $(vectorfunc)(a::$atype, b::$btype)
                    res = similar($(asim ? :a : :b), $outtype, promote_shape(size(a), size(b)))
                    for i = 1:length(a)
                        res[i] = $(scalarfunc)(a[i], b[i])
                    end
                    res
                end
            end
            for (asim, atype, btype) in ((true, :DataArray, :AbstractDataArray),
                                         (false, :AbstractDataArray, :DataArray),
                                         (true, :AbstractDataArray, :AbstractDataArray),
                                         (true, :AbstractDataArray, :AbstractArray),
                                         (false, :AbstractArray, :AbstractDataArray))
        }...,
    ))
end

macro dataframe_binary(f)
    esc(quote
        function $(f)(a::DataFrame, b::DataFrame)
            if size(a) != size(b); error("argument dimensions must match"); end
            DataFrame([$(f)(a[i], b[i]) for i=1:size(a, 2)], deepcopy(index(a)))
        end
        @swappable $(f)(a::DataFrame, b::Union(Number, String)) =
            DataFrame([$(f)(a[i], b) for i=1:size(a, 2)], deepcopy(index(a)))
        @swappable $(f)(a::DataFrame, b::NAtype) =
            DataFrame([$(f)(a[i], b) for i=1:size(a, 2)], deepcopy(index(a)))
    end)
end

# Unary operators, NA
for f in unary_operators
    @eval $(f)(d::NAtype) = NA
end

# Unary operators, DataArrays. Definitions in base should be adequate
# for AbstractDataArrays. These are just optimizations
!(d::DataArray{Bool}) = DataArray(!d.data, copy(d.na))
-(d::DataArray) = DataArray(-d.data, copy(d.na))

# Unary operators, DataFrames
@dataframe_unary !
@dataframe_unary -
# As in Base, these are identity operators
for f in (:(+), :(*))
    @eval $(f)(d::DataFrame) = d
end

# Treat ctranspose and * in a special way for now
for f in (:ctranspose, :transpose)
    @eval $(f)(d::DataArray) = DataArray($(f)(d.data), d.na')
end

# Propagates NA's
# For a dissenting view,
# http://radfordneal.wordpress.com/2011/05/21/slowing-down-matrix-multiplication-in-r/
# But we're getting 10x R while maintaining NA's
for (adata, bdata) in ((true, false), (false, true), (true, true))
    @eval begin
        function (*)(a::$(adata ? :(Union(DataVector, DataMatrix)) : :(Union(Vector, Matrix))),
                     b::$(bdata ? :(Union(DataVector, DataMatrix)) : :(Union(Vector, Matrix))))
            c = $(adata ? :(a.data) : :a) * $(bdata ? :(b.data) : :b)
            res = DataArray(c, falses(size(c)))
            # Propagation can be made more efficient by storing record of corrupt
            # rows and columns, then doing fast edits.
            $(if adata
                quote
                    n1 = size(a, 1)
                    p1 = size(a, 2)
                    corrupt_rows = falses(n1)
                    for j in 1:p1, i in 1:n1
                        # Propagate NA's
                        # Corrupt all rows based on i
                        corrupt_rows[i] |= a.na[i, j]
                    end
                    res.na[corrupt_rows, :] = true
                end
            end)
            $(if bdata
                quote
                    n2 = size(b, 1)
                    p2 = size(b, 2)
                    corrupt_cols = falses(p2)
                    for j in 1:p2, i in 1:n2
                        # Propagate NA's
                        # Corrupt all columns based on j
                        corrupt_cols[j] |= b.na[i, j]
                    end
                    res.na[:, corrupt_cols] = true
                end
            end)
            res
        end
    end
end

#
# Elementary functions
#
# XXX: The below should be revisited once we have a way to infer what
# the proper return type of an array should be.

# One-argument elementary functions that return the same type as their
# inputs
for f in (:abs, :sign)
    @eval begin
        $(f)(::NAtype) = NA
        @dataarray_unary $(f) Number T
        @dataframe_unary $(f)
    end
end

# One-argument elementary functions that always return floating points
for f in (:acos, :acosh, :asin, :asinh, :atan, :atanh, :sin, :sinh, :cos,
          :cosh, :tan, :tanh, :exp, :exp2, :expm1, :log, :log10, :log1p,
          :log2, :exponent, :sqrt, :gamma, :lgamma, :digamma, :erf, :erfc)
    @eval begin
        ($f)(::NAtype) = NA
        @dataarray_unary $(f) FloatingPoint T
        @dataarray_unary $(f) Real Float64
        @dataframe_unary $(f)
    end
end

# Elementary functions that take varargs
for f in (:round, :ceil, :floor, :trunc)
    @eval begin
        ($f)(::NAtype, args::Integer...) = NA

        # ambiguity
        @dataarray_unary $(f) Real T 1
        @dataarray_unary $(f) Real T 2
        @dataarray_unary $(f) Real T

        function $(f){T<:Real}(d::DataArray{T}, args::Integer...)
            data = similar(d.data)
            for i = 1:length(data)
                if !d.na[i]
                    data[i] = $(f)(d[i], args...)
                end
            end
            DataArray(data, copy(d.na))
        end
        function $(f){T<:Real}(adv::AbstractDataArray{T}, args::Integer...)
            res = similar(adv)
            for i = 1:length(adv)
                res[i] = ($f)(adv[i], args...)
            end
            res
        end
        $(f)(d::DataFrame, args::Integer...) =
            DataFrame([$(f)(d[i], args...) for i=1:size(d, 2)], deepcopy(index(d)))
    end
end

#
# Bit operators
#

@swappable (&)(a::NAtype, b::Bool) = b ? NA : false
@swappable (|)(a::NAtype, b::Bool) = b ? true : NA
@swappable ($)(a::NAtype, b::Bool) = NA

# To avoid ambiguity warning
@swappable (|)(a::NAtype, b::Function) = NA

for f in (:&, :|, :$)
    @eval begin
        # Scalar with NA
        ($f)(::NAtype, ::NAtype) = NA
        @swappable ($f)(::NAtype, b) = NA
        # DataFrame
        @dataframe_binary $(f)
    end
end

# DataArray with DataArray
(&)(a::DataArray{Bool}, b::DataArray{Bool}) =
    DataArray(a.data & b.data, (a.na & b.data) | (b.na & a.data))
(|)(a::DataArray{Bool}, b::DataArray{Bool}) =
    DataArray(a.data | b.data, (a.na & !b.data) | (b.na & !a.data))
($)(a::DataArray{Bool}, b::DataArray{Bool}) =
    DataArray(a.data $ b.data, a.na | b.na)

# DataArray with non-DataArray
# Need explicit definition for BitArray to avoid ambiguity
for t in (:(BitArray), :(Union(AbstractArray{Bool}, Bool)))
    @eval begin
        @swappable (&)(a::DataArray{Bool}, b::$t) = DataArray(a.data & b, a.na & b)
        @swappable (|)(a::DataArray{Bool}, b::$t) = DataArray(a.data | b, a.na & !b)
        @swappable ($)(a::DataArray{Bool}, b::$t) = DataArray(a.data $ b, a.na)
    end
end

#
# Comparison operators
#

isless(::NAtype, ::NAtype) = false
isless(::NAtype, b) = true
isless(a, ::NAtype) = false

# This is for performance only; the definition in Base is sufficient
# for AbstractDataArrays
function isequal(a::DataArray, b::DataArray)
    if size(a) != size(b) || a.na != b.na
        return false
    end
    for i = 1:length(a)
        if !a.na[i] && a.data[i] != b.data[i]
            return false
        end
    end
    return true
end
function isequal(df1::AbstractDataFrame, df2::AbstractDataFrame)
    if size(df1, 2) != size(df2, 2)
        return false
    end
    for idx in 1:size(df1, 2)
        if !isequal(df1[idx], df2[idx])
            return false
        end
    end
    return true
end

for sf in scalar_comparison_operators
    vf = symbol(".$sf")

    @eval begin
        # Array with NA
        @swappable $(vf){T,N}(::NAtype, b::AbstractArray{T,N}) =
            DataArray(Array(Bool, size(b)), trues(size(b)))

        # Scalar with NA
        $(vf)(::NAtype, ::NAtype) = NA
        $(sf)(::NAtype, ::NAtype) = NA
        @swappable $(vf)(::NAtype, b) = NA
        @swappable $(sf)(::NAtype, b) = NA

        @dataarray_binary_scalar $vf $sf Bool
        @dataarray_binary_array $vf $sf Bool

        $(sf)(a::AbstractDataArray, b::AbstractDataArray) =
            error("$sf not defined for AbstractDataArrays. Try $vf")

        # DataFrame
        @dataframe_binary $vf
    end
end

#
# Binary operators
#

# Necessary to avoid ambiguity warnings
.^(::MathConst{:e}, B::DataArray) = exp(B)
.^(::MathConst{:e}, B::AbstractDataArray) = exp(B)

for f in arithmetic_operators
    @eval begin
        # Array with NA
        @swappable $(f){T,N}(::NAtype, b::AbstractArray{T,N}) =
            DataArray(Array(T, size(b)), trues(size(b)))

        # Scalar with NA
        ($f)(::NAtype, ::NAtype) = NA
        @swappable ($f)(d::NAtype, x::Number) = NA

        # DataArray with scalar
        @dataarray_binary_scalar $f $f promote_type(eltype(a), eltype(b))

        # DataFrame
        @dataframe_binary $f
    end
end

^(::NAtype, ::NAtype) = NA
^(a, ::NAtype) = NA
^(::NAtype, ::Integer) = NA
^(::NAtype, ::Number) = NA

for (vf, sf) in ((:(+), :(+)), (:(.+), :(+)), (:(-), :(-)), (:(.-), :(-)), (:(.*), :(*)),
                 (:(.^), :(^)))
    @eval begin
        # Necessary to avoid ambiguity warnings
        @swappable ($vf)(A::BitArray, B::AbstractDataArray) = ($vf)(bitunpack(A), B)
        @swappable ($vf)(A::BitArray, B::DataArray) = ($vf)(bitunpack(A), B)

        @dataarray_binary_array $vf $sf promote_type(eltype(a), eltype(b))
    end
end

@swappable ./(A::BitArray, B::AbstractDataArray) = ./(bitunpack(A), B)
@swappable ./(A::BitArray, B::DataArray) = ./(bitunpack(A), B)

@dataarray_binary_array (./) (/) isa(a, FloatingPoint) || isa(b, FloatingPoint) ?
                                 promote_type(a, b) : Float64

for f in biscalar_operators
    @eval begin
        ($f)(::NAtype, ::NAtype) = NA
        @swappable $(f)(::Number, ::NAtype) = NA
    end
end

for f in pairwise_vector_operators
    @eval function ($f)(dv::DataVector)
        n = length(dv)
        new_data = ($f)(dv.data)
        new_na = falses(n - 1)
        new_na[1] = dv.na[1]
        for i = 2:(n - 1)
            if dv.na[i]
                new_na[i - 1] = true
                new_na[i] = true
            end
        end
        new_na[n - 1] = new_na[n - 1] || dv.na[n]
        return DataArray(new_data, new_na)
    end
end

for f in cumulative_vector_operators
    @eval function ($f)(dv::DataVector)
        new_data = ($f)(dv.data)
        new_na = falses(length(dv))
        hitna = false
        for i = 1:length(dv)
            if dv.na[i]
                hitna = true
            end
            if hitna
                new_na[i] = true
            end
        end
        return DataArray(new_data, new_na)
    end
end

for f in [unary_vector_operators; ffts]
    @eval ($f)(dv::DataVector) = any(dv.na) ? NA : ($f)(dv.data)
end

for f in binary_vector_operators
    @eval ($f)(dv1::DataVector, dv2::DataVector) =
            any(dv1.na) || any(dv2.na) ? NA : ($f)(dv1.data, dv2.data)
end

for f in (:minimum, :maximum, :prod, :sum, :mean, :median, :std, :var, :norm)
    colf = symbol("col$(f)s")
    rowf = symbol("row$(f)s")
    @eval begin
        function ($colf)(df::AbstractDataFrame)
            p = ncol(df)
            res = DataFrame()
            for j in 1:p
                res[j] = DataArray(($f)(df[j]))
            end
            colnames!(res, colnames(df))
            return res
        end
        function ($colf)(dm::AbstractDataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(p)
            for j in 1:p
                res[j] = ($f)(dm[:, j])
            end
            return res
        end
        function ($rowf)(dm::DataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(n)
            for i in 1:n
                res[i] = ($f)(DataArray(reshape(dm.data[i, :], p), reshape(dm.na[i, :], p)))
            end
            return res
        end
    end
end

#
# Boolean operators
#

function all(dv::AbstractDataArray{Bool})
    for i in 1:length(dv)
        if dv.na[i]
            return NA
        end
        if !dv.data[i]
            return false
        end
    end
    true
end

function all(dv::AbstractDataArray{Bool})
    for i in 1:length(dv)
        if isna(dv[i])
            return NA
        end
        if !dv[i]
            return false
        end
    end
    true
end

function any(dv::DataArray{Bool})
    has_na = false
    for i in 1:length(dv)
        if !dv.na[i]
            if dv.data[i]
                return true
            end
        else
            has_na = true
        end
    end
    has_na ? NA : false
end

function any(dv::AbstractDataArray{Bool})
    has_na = false
    for i in 1:length(dv)
        if !isna(dv[i])
            if dv[i]
                return true
            end
        else
            has_na = true
        end
    end
    has_na ? NA : false
end

function all(df::AbstractDataFrame)
    for i in 1:size(df, 2)
        x = all(df[i])
        if isna(x)
            return NA
        end
        if !x
            return false
        end
    end
    true
end

function any(df::AbstractDataFrame)
    has_na = false
    for i in 1:size(df, 2)
        x = any(df[i])
        if !isna(x)
            if x
                return true
            end
        else
            has_na = true
        end
    end
    has_na ? NA : false
end

function Stats.range{T}(dv::AbstractDataVector{T})
    return DataVector[minimum(dv), maximum(dv)]
end

function rle{T}(v::AbstractVector{T})
    n = length(v)
    current_value = v[1]
    current_length = 1
    values = similar(v, n)
    total_values = 1
    lengths = Array(Int16, n)
    total_lengths = 1
    for i in 2:n
        if v[i] == current_value
            current_length += 1
        else
            values[total_values] = current_value
            total_values += 1
            lengths[total_lengths] = current_length
            total_lengths += 1
            current_value = v[i]
            current_length = 1
        end
    end
    values[total_values] = current_value
    lengths[total_lengths] = current_length
    return (values[1:total_values], lengths[1:total_lengths])
end

function rle{T}(v::AbstractDataVector{T})
    n = length(v)
    current_value = v[1]
    current_length = 1
    values = DataArray(T, n)
    total_values = 1
    lengths = Array(Int16, n)
    total_lengths = 1
    for i in 2:n
        if isna(v[i]) || isna(current_value)
            if isna(v[i]) && isna(current_value)
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        else
            if v[i] == current_value
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        end
    end
    values[total_values] = current_value
    lengths[total_lengths] = current_length
    return (values[1:total_values], lengths[1:total_lengths])
end

## inverse run-length encoding
function inverse_rle{T, I <: Integer}(values::AbstractVector{T}, lengths::Vector{I})
    total_n = sum(lengths)
    pos = 0
    res = similar(values, total_n)
    n = length(values)
    for i in 1:n
        v = values[i]
        l = lengths[i]
        for j in 1:l
            pos += 1
            res[pos] = v
        end
    end
    return res
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to ECI to ECEF transformations using satellite state vector.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Get the current EOP Data.
#
# TODO: The EOP data obtained from IERS website does not match the values in the
# examples in [1]. However, this should be enough, because 1) the individual
# functions at the low level are tested using the same values of [1], and 2) the
# difference is smaller than 30 cm.

eop_iau1980  = read_iers_eop("./eop_IAU1980.txt",  Val(:IAU1980))
eop_iau2000a = read_iers_eop("./eop_IAU2000A.txt", Val(:IAU2000A))

# File: ./src/transformations/sv_eci_to_ecef.jl
# =============================================

# Functions: sv_eci_to_ecef
# -------------------------

# The rotations functions were already heavily tested in `ecef_to_eci.jl`,
# `fk5/fk5.jl`, and `iau2006/iau2006.jl`. Hence, here we will do only some minor
# testing involving the following transformations:
#
#   GCRF  => ITRF (FK5)
#   J2000 => PEF  (FK5, no corrections)
#   GCRF  => ITRF (IAU-2010)
#   GCRF  => TIRS (IAU-2010, no corrections)

################################################################################
#                                  IAU-76/FK5
################################################################################

# GCRF to ICRF
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction [1, p. 230].
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf =  5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#   v_gcrf =    -4.7432201570 i +    0.7905364970 j +    5.533755270  k [km/s]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#   v_itrf =    -3.225636520  i -    2.872451450  j +    5.531924446  k [km/s]
#
################################################################################

@testset "Function sv_eci_to_ecef GCRF => ITRF" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]
    v_gcrf = [-4.7432201570; 0.7905364970; 5.5337557270]
    sv_gcrf = orbsv(JD_UTC, r_gcrf, v_gcrf)
    sv_itrf = sv_eci_to_ecef(sv_gcrf, GCRF(), ITRF(), JD_UTC, eop_iau1980)

    @test sv_itrf.t === JD_UTC

    @test sv_itrf.r[1] ≈ -1033.4793830 atol=3e-4
    @test sv_itrf.r[2] ≈ +7901.2952754 atol=3e-4
    @test sv_itrf.r[3] ≈ +6380.3565958 atol=3e-4

    @test sv_itrf.v[1] ≈ -3.225636520  atol=8e-7
    @test sv_itrf.v[2] ≈ -2.872451450  atol=8e-7
    @test sv_itrf.v[3] ≈ +5.531924446  atol=8e-7
end

# J2000 to PEF
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction [1, p. 230].
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   UT1    = UTC - 0.4399619/86400
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#   v_j2000 =  -4.7432196000 i + 0.7905366000    j + 5.5337561900    k [km/s]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_pef  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function sv_eci_to_ecef J2000 => PEF" begin
    JD_UT1 = date_to_jd(2004,4,6,7,51,28.386009) - 0.4399619/86400

    r_j2000  = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000  = [-4.7432196000; 0.7905366000; 5.5337561900]
    sv_j2000 = orbsv(JD_UT1, r_j2000, v_j2000)
    sv_pef   = sv_eci_to_ecef(sv_j2000, J2000(), PEF(), JD_UT1)

    @test sv_pef.t === JD_UT1

    @test sv_pef.r[1] ≈ -1033.47503130 atol=1e-7
    @test sv_pef.r[2] ≈ +7901.30558560 atol=1e-7
    @test sv_pef.r[3] ≈ +6380.34453270 atol=1e-7

    @test sv_pef.v[1] ≈ -3.2256327470  atol=1e-7
    @test sv_pef.v[2] ≈ -2.8724425110  atol=1e-7
    @test sv_pef.v[3] ≈ +5.5319312880  atol=1e-7
end

################################################################################
#                                IAU-2006/2010
################################################################################

# GCRF to ITRF
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-2010 reduction [1, p. 220].
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895290   i + 6123.01139910   j + 6378.13693380   k [km]
#   v_gcrf =   -4.7432201610 i +    0.7905364950 j +    5.533755240  k [km/s]
#
# one gets:
#
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#   v_itrf =    -3.225636520  i -    2.872451450  j +    5.531924446  k [km/s]
#
################################################################################

@testset "Function sv_eci_to_ecef GCRF => ITRF" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    r_gcrf = [5102.50895290; 6123.01139910; 6378.13693380]
    v_gcrf = [-4.7432201610; 0.7905364950; 5.5337557240]
    sv_gcrf = orbsv(JD_UTC, r_gcrf, v_gcrf)
    sv_itrf = sv_eci_to_ecef(sv_gcrf, GCRF(), ITRF(), JD_UTC, eop_iau2000a)

    @test sv_itrf.t === JD_UTC

    @test sv_itrf.r[1] ≈ -1033.4793830 atol=3e-4
    @test sv_itrf.r[2] ≈ +7901.2952754 atol=3e-4
    @test sv_itrf.r[3] ≈ +6380.3565958 atol=3e-4

    @test sv_itrf.v[1] ≈ -3.225636520  atol=8e-7
    @test sv_itrf.v[2] ≈ -2.872451450  atol=8e-7
    @test sv_itrf.v[3] ≈ +5.531924446  atol=8e-7
end

# GCRF to TIRS
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-2010 reduction [1, p. 220].
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   UT1    = UTC - 0.4399619/86400
#   r_gcrf = 5102.50895290   i + 6123.01139910   j + 6378.13693380   k [km]
#   v_gcrf =   -4.7432201610 i +    0.7905364950 j +    5.533755240  k [km/s]
#
# one gets:
#
#   r_tirs  = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_tirs  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function sv_eci_to_ecef TIRS => GCRF" begin
    JD_UT1 = date_to_jd(2004,4,6,7,51,28.386009) - 0.4399619/86400

    r_gcrf = [5102.50895290; 6123.01139910; 6378.13693380]
    v_gcrf = [-4.7432201610; 0.7905364950; 5.5337557240]
    sv_gcrf = orbsv(JD_UT1, r_gcrf, v_gcrf)
    sv_tirs = sv_eci_to_ecef(sv_gcrf, GCRF(), TIRS(), JD_UT1)

    @test sv_tirs.t === JD_UT1

    @test sv_tirs.r[1] ≈ -1033.47503120 atol=3e-4
    @test sv_tirs.r[2] ≈ +7901.30558560 atol=3e-4
    @test sv_tirs.r[3] ≈ +6380.34453270 atol=3e-4

    @test sv_tirs.v[1] ≈ -3.2256327470  atol=8e-7
    @test sv_tirs.v[2] ≈ -2.8724425110  atol=8e-7
    @test sv_tirs.v[3] ≈ +5.5319312880  atol=8e-7
end

ext: jl


text:
module TransformReduce

push!(LOAD_PATH, joinpath(@__DIR__, "../", "core"))

using Base.Threads
using GenericBlas

@inline function transform(policy::Any, range::Int64, transformer::Any)

@threads for i in 1:range
           transformer(i)
         end
end #transform

@inline function transform_reduce(policy::Any, range::Int64, init::Any, reducer::Any, transformer::Any)
         local res = init
@threads for i in 1:range
           local val = transformer(i)
           res       = reducer(res, val)
         end
end #transform_reduce

end #TrnasformReduce

max_stars_repo_name: alexstrel/solver-tuner

max_issues_repo_name: alexstrel/solver-tuner

max_forks_repo_name: alexstrel/solver-tuner


text:
#Basic Example of IRLS in a noiseless situation

using PyPlot
using CompressedSensing

srand(0)

#Signal Dimensions
n=1000
#measurement dimension
m=400
#Approximate Sparsity
k=100

#Create a random sparse signal
Signal = zeros(n)
Signal[int(rand(k)*(n-1)+1)] = randn(k)

#Create a random measurement matrix
MM = randn(m,n)

#Measure the signal using the measurement matrix
Measurement = MM*Signal

#Reconstruct using IRLS
Reconstruction = IRLS(MM,Measurement,verbose=true)

#plot the results
bar(1:length(Signal),Signal,width=.01) #easiest method for plotting thin lines
SigHandle = plt.Line2D([], [],label="Signal",color="black") #by using bar() plot above, the legend
						#gets all funky, so to fix it make a handle that contains the correct
						#information
ReconPlot=scatter(1:length(Reconstruction),Reconstruction,label="IRLS Reconstruction")
xlim(0,length(Signal))
legend(loc=1,handles=[SigHandle,ReconPlot])
savefig("Example_1_Fig_1.png")

max_stars_repo_path: examples/Example 1.jl

max_stars_repo_name: dahlend/CompressedSensing

max_issues_repo_path: examples/Example 1.jl

max_forks_repo_path: examples/Example 1.jl


text:
md"""
Лекция 8

## Проектираование на основе построения иерархии типов, параметрические типы

В лекции 7 была рассмотрена следующая иерархию типов, позволившая нам писать обобщенный код:

AbstractRobot:
|                        move!(::AbstractCoordRobot, ::HorizonSide),
|                        isborder(::AbstractCoordRobot, ::HorizonSide),
|                        putmarker!(::AbstractCoordRobot),
|                        ismarker(::AbstractCoordRobot),
|                        temperature(::AbstractCoordRobot)
|
|                        snake(::AbstractCoordRobot, NTuple(2,HorizonSide))
|
|--AbsctractBorderRobot: movements!(::AbsctractBorderRobot, ::HorisonSide),
|  |                     movements!(::AbsctractBorderRobot, ::HorisonSide, ::Integer),
|  |                          (предполагается наличие функции try_move!(::Robot, ::HorizonSide))
|  |
|  |--PutmarkersBorderRobot: try_move!(::PutmarkersBorderRobot, ::HorizonSide), get(::PutmarkerRobot)
|  |--CountmarkerBorderRobot: try_move!(::PutmarkersBorderRobot, ::HorizonSide), get(::CountmarkerRobot)
|  |--....
|
|-- AbstractCoordRobot:
|  |                    get_coords(::AbstractCoordRobot),
|  |                    move!(::AbstractCoordRobot, ::HorizonSide)
|  |
|  |- CoordRobot:
|  |              get(::AbstractCoordRobot)
|  |              get_coords(::AbstractCoordRobot), move!(::AbstractCoordRobot, ::HorizonSide)
|
|--....

Допустим теперь, что при перемещении робота по полю с обходом перегородок понадобилось ещё и вычислять текущие
координаты робота. Например, допустим, что требуется не постчитать число маркеров на поле, а получить их
координаты (например, записать в массив).

Как можно было бы поступить в таком случае?

Например, можно было бы спроектировать, следующий тип данных.

"""
struct CoordmarkersBorderRobot <: AbstractBorderRobot
    robot::Robot
    coord::Coord
    marker_coord::Vector{NTuple{2,Int}}
end

md"""
Пример создания объекта типа CoormarkersBorderRobot:
"""
r=Robot(animate=true)

r = CoormarkersBorderRobot(r, Coord(), Int[])

md"""

При этом, однако, пришлось бы фактически повторить определения функций
    get_coords(::AbstractCoordRobot), move!(::AbstractCoordRobot, ::HorizonSide),
(уже имеющихся для типа AbstractCoordRobot) т.е.
"""
get_coord(robot::CoordmarkersBorderRobot) = get_coord(robot.coord) #...get_coord -???
function move!(robot::CoordmarkersBorderRobot, side::HorizonSide)
    move!(robot,side)
    move!(robot.coord, side)
end

md"""
И, конечно, при этом потребуется ещё переопределить и функцию try_move! следующим образом.
"""

function try_move!(robot::CoordmarkersBorderRobot, side::HorizonSide)
    try_move!(robot.robot, side)
    if ismarker(robot)
        push!(robot.marker_coord, get_coord(robot))
    end
end

md"""

Однако необходимости писать функции
    get_coords(::CoordmarkersBorderRobot), move!(::CoordmarkersBorderRobot, ::HorizonSide),
(повторяя фактически опрделения одноименных функций для типа AbstractCoordRobot)
можно было бы избежать, если бы вместо типа AbsctractBorderRobot нами был бы спроектирован соответствующий
параметрический тип.

Вот как это могло бы выглядеть.

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

AbstractRobot:
|                        move!(::AbstractRobot, ::HorizonSide),
|                        isborder(::AbstractRobot, ::HorizonSide),
|                        putmarker!(::AbstractRobot),
|                        ismarker(::AbstractRobot),
|                        temperature(::AbstractRobot)
|
|--AbsctractBorderRobot{TypeRobot}: # - параметрический тип, где TypeRobot - это параметр этого типа
|  |                     movements!(::AbsctractBorderRobot{TypeRobot},::HorisonSide),
|  |                     movements!(::AbsctractBorderRobot{TypeRobot},::HorisonSide,Integer),
|  |                     try_move!(::AbstractBorderRobot{TypeRobot}, ::HorizonSide)
|  |                     get(robot::AbsctractBorderRobot{TypeRobot}) = robot.robot # get(robot::Robot) = robot
|  |
|  |--PutmarkersBorderRobot{TypeRobot}: # - теперь это особая разновидность абстрактных типов (все параметрические типы - абстрактные, хотя и определяются как структура)
|  |                    try_move!(::PutmarkersBorderRobot{TypeRobot}, ::HorizonSide),
|  |
|  |--CountmarkerBorderRobot:
|  |                    try_move!(::CountmarkerBorderRobot, ::HorizonSide),
|  |
|  |--....
|
|--AbstractCoordRobot
|  |
|  |-CoordRobot{TypeRobot}:
|  |                    get(robot::CoordRobot{TypeRobot}) = robot.robot  #::Robot
|
|--PutmarkerRobot
|--CountmarkersRobot
|--....

get(robot::Robot) = robot # - теперь это необходимо доопределить
#------------------------------------------------------

"""

abstract type AbstracrBorderRobot{TypeRobot} <: AbstractRobot end

function movements!(robot::AbsctractBorderRobot{TypeRobot},side::HorisonSide)
    n=0
    while try_move!(get(robot), side)
        n+=1
    end
    return n
end

function movements!(robot::AbstractBorderRobot{TypeRobot}, side::HorizonSide, num_steps::Integer)
    for _ in 1:num_steps
        try_move!(get(robot), side)
    end
end

#----------------------------------------------
get(robot::Robot) = robot

#--------------------------------------------------------
struct CoordmarkersBorderRobot{TypeRobot} <: AbstractBorderRobot{TypeRobot}
    robot::TypeRobot
    markers_coord::Vector{NTuple{2,Int}}
end

get(robot::CoordmarkersBorderRobot{CoordRobot})=get(robot.robot)

function try_move!(robot::CoordmarkersBorderRobot{CoordRobot}, side::HorizonSide)
    try_move!(get(robot), side)
    if ismarker(robot)
        push!(robot.marker_coord, get_coord(robot))
    end
end

md"""
ЗАМЕЧАНИЕ

Если определен тип CoordRobot
то тип CoordmarkersBorderRobot{CoordRobot} будет конкретными

Но тип CoordmarkersBorderRobot{TypeRobot} будет абстрактным параметрическим типом (если TypeRobot - это неопределенный параметр).

В данном случае функции get и try_move! были определены для конкретного типа CoordmarkersBorderRobot{CoordRobot}
(т.е. эти функции являются узко специализированными, не являются обобщенными)
#-----------------------------------------------------

ПРИМЕРЫ СОЗДАНИЯ ОБЪЕКТОВ
"""
r = Robot()

r1 = CoordmarkersBorderRobot{Robot}(r, Int[])
# - r1 - на самом деле координаты не отслеживает
#        и поэтому для типа CoordmarkersBorderRobot{Robot} нами не была определена функция try_move
#        (как и функция get)
#   Но тем не менее создать объект r1 можно, но использовать его для наших целей не получится

r2 = CoordRobot(r) # по умолчанию устанавливаются нулевые начальные координаты

get(r2) #--> r2.robot (::Robot)

r3 = CoordmarkersBorderRobot{CoordRobot}(r2, Int[]) # - отслеживает координаты

get(r3) #--> get(r3.robot) --> r3.robot.robot (::Robot)

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

md"""
В виду бесполезности типа CoordmarkersBorderRobot{Robot}, а также бесполезности всех остальных возможных подтипов
типа CoordmarkersBorderRobot{TypeRobot} за исключением типа
CoordmarkersBorderRobot{CoordRobot}
тип CoordmarkersBorderRobot не следовало делать параметрическим.

Правильно было бы сделать так.
"""

struct CoordmarkersBorderRobot <: AbstractBorderRobot{CoordRobot}
    robot::CoordRobot
    markers_coord::Vector{NTuple{2,Int}}
end

get(robot::CoordmarkersBorderRobot) = get(robot.robot)

function try_move!(robot::CoordmarkersBorderRobot, side::HorizonSide)
    try_move!(get(robot), side)
    if ismarker(robot)
        push!(robot.marker_coord, get_coord(robot))
    end
end

md"""
ЗАДАЧА
ДАНО: Робот - в юго-западном углу пояля, на котором имеются внутренние прямоугольные или прямолинейные перегородки
РЕЗУЛЬТАТ: Функция вернула массив с координатами замаркированных клеток
"""

function markers_coord(r::Robot)
    coord_roobot = CoordRobot(r)
    robot = CoordmarkersBorderRobot(coord_roobot, Int[])
    snake!(robot)
    return robot.markers_coord
end

md"""
Если определить ещё следующую функцию
"""

get_markers_coord(robot::CoordmarkersBorderRobot) = robot.markers_coord

md"""

и если при этом функция snake! возвращает ссылку на робота (принимает её и её же возвращает),
то определение функции markers_coord(::Robot) могло бы быть совсем кратким:
"""

markers_coord(r::Robot) = get_markers_coord(snake!(CoordmarkersBorderRobot(CoordRobot(r), Int[])))

md"""
Для удобства можно было бы также определить ещё один конструктор:
CoordmarkersBorderRobot(r::Robot) = CoordmarkersBorderRobot(CoordRobot(r), Int[]))

тогда определение markers_coord(::Robot) могло бы стать еще более лаконичным:
"""
markers_coord(r::Robot) = get_markers_coord(snake!(CoordmarkersBorderRobot(r))

md"""
Пример проектирования ещё одного типа
"""

struct PutmarkersBorderRobot{TypeRobot} <: AbstractBorderRobot{TypeRobot}
    robot::TypeRobot
end

md"""
Здесь, в отличие от CoormarkersBorderRobot тип PutmarkersBorderRobot{TypeRobot} определен как параметрический.

Это не лишено смысла, т.к. робот-постановщик маркеров можно быть реализован как на базе типа Robot, так и
на базе какого-нибудь другого типа, например, - CoordRobot

При этом, например, функции get и try_move могут быть определены как обобщенные функции, так и как узкоспециализированные
"""
get(robot::PutmarkersBorderRobot{TypeRobot}) where TypeRobot = get(robot.robot)

function try_move!(robot::PutmarkersBorderRobot{TypeRobot}, side::HorizonSide) where TypeRobot
    if try_move!(get(robot), side)
        putmarker!(robot)
    end
end

Так определенные функции будут обобщенными, поскольку тип аргумента PutmarkersBorderRobot{TypeRobot} является
 абстрактным: зависит от параметра TypeRobot. То, что TypeRobot - это параметр, а не имя какого-то ранее опрделенного
 типа (абстрактного или конкретного) указывается с помощью ключевого слова where.

Следующие определния являются определениями узкоспециализированных функций (не обобщенных)

get(robot::PutmarkersBorderRobot{Robot}) = get(robot.robot)

function try_move!(robot::PutmarkersBorderRobot{Robot}, side::HorizonSide)
    if try_move!(get(robot), side)
        putmarker!(robot)
    end
end

В отличие от сделанных выше определений обобщенных функций get, try_move!, такие определения потребовались
для каждого фактического значения параметра TypeRobot.

md"""
# Элементы функционального программирования

## Функции - как объекты 1-го класса
"""
x = 0

s = sin  # <: Function

md"""
Пусть, например, мы хотим перемещать робота в заданном направлении до тех пор, пока не наступит некоторое событие,
выраженное соответствующим условием.

Это реализуется следующей функцией:
"""

function movements!(condition::Function, r, side)
    n=0
    while !condition(r)
        move!(r,side)
        n+=1
    end
    return n
end

#Теперь эту функцию можно вызывать, например, так:

movements!(isborder, r, side)

md"""
## Анонимные функции
Анонимные функции представляют собой функциональные значения (<: Function) без имени
 (подобно тому, как бывают значения типа Int без имени)

Например:

condition = robot -> !isborder(robot)

Здесь переменной condition присвоено некоторое функциональное значение (значение анонимной функции)
"""


text:
module dispersion

using Interpolations

_fraunhofer_lines = Dict(
  "i" => ("Hg", 0.36501, "blue"),
  "h" => ("Hg", 0.40466, "blue"),
  "g" => ("Hg", 0.43585, "blue"),
  "F'" => ("Cd", 0.47999, "green"),
  "F" => ("Hβ", 0.486134, "infrared"),
  "e" => ("Hg", 0.546073, "infrared"),
  "d" => ("He", 0.58756, "infrared"),
  "D" => ("Na", 0.5893, "red"),
  "C'" => ("Cd", 0.64385, "red"),
  "C" => ("Hα", 0.656281, "red"),
  "r" => ("He", 0.70652, "ultraviolet"),
  "A'" => ("K", 0.7682, "violet"),
  "s" => ("Cs", 0.85211, "yellow"),
  "t" => ("Hg", 1.01398, "yellow"),
  "t2" => ("Ni", 0.299444, ""),
  "T" => ("Fe", 0.302108, ""),
  "P" => ("Ti+", 0.336112, ""),
  "N" => ("Fe", 0.358121, ""),
  "L" => ("Fe", 0.382044, ""),
  "K" => ("Ca+", 0.3933666, ""),
  "H" => ("Ca+", 0.396847, ""),
  "h2" => ("Hδ", 0.410175, ""),
  "G2" => ("Ca", 0.430774, ""),
  "G" => ("Fe", 0.43079, ""),
  "G'" => ("Hγ", 0.434047, ""),
  "e2" => ("Fe", 0.438355, ""),
  "d2" => ("Fe", 0.466814, ""),
  "c" => ("Fe", 0.495761, ""),
  "b4" => ("Mg", 0.516733, ""),
  "b3" => ("Fe", 0.516891, ""),
  "b2" => ("Mg", 0.51727, ""),
  "b1" => ("Mg", 0.518362, ""),
  "E2" => ("Fe", 0.527039, ""),
  "D3" => ("He", 0.5875618, ""),
  "D2" => ("Na", 0.588995, ""),
  "D1" => ("Na", 0.589592, ""),
  "a" => ("O2", 0.627661, ""),
  "B" => ("O2", 0.686719, ""),
  "A" => ("O2", 0.75937, ""),
  "Z" => ("O2", 0.822696, ""),
  "y" => ("O2", 0.898765, "")
)

line_λum(line::AbstractString) = get(_fraunhofer_lines,line,NaN)[2]
line_source(line::AbstractString) = get(_fraunhofer_lines,line,"")[1]
line_color(line::AbstractString) = get(_fraunhofer_lines,line,"")[3]
line_standards() = keys(filter((k,v)->!isempty(v[3]), _fraunhofer_lines))

# standard dispersion calculation wavelengths
const λum_D = line_λum("D")
const λum_C = line_λum("C")
const λum_F = line_λum("F")

# forumulas used by the database
@enum(INDEXFN, Sellmeier =1, Sellmeier2 =2, Polynomial =3, RFInfo =4,
  Cauchy =5, Gases =6, Herzberger =7, Retro =8, Exotic =9)

# Sellmeier (formula 1)
function refindex{T<:Real}(::Type{Val{Sellmeier}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local λsq::T = λum^2
  local nsq::T = one(T) + C[1]
  for i = 2:2:N-1
    nsq += C[i] * λsq / (λsq - C[i+1]^2)
  end
  return sqrt(nsq)
end

# Sellmeier2 (formula 2)
function refindex{T<:Real}(::Type{Val{Sellmeier2}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local λsq::T = λum^2
  local nsq::T = one(T) + C[1]
  for i = 2:2:N-1
    nsq += C[i] * λsq / (λsq - C[i+1])
  end
  sqrt(nsq)
end

# Polynomial (formula 3)
function refindex{T<:Real}(::Type{Val{Polynomial}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local nsq::T = C[1]
  for i = 2:2:N-1
    nsq += C[i] * λum^C[i+1]
  end
  sqrt(nsq)
end

# RFInfo (formula 4)
function refindex{T<:Real}(::Type{Val{RFInfo}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=11 "index equation requires an odd number of coefficients ≥ 11"
  local λsq::T = λum^2
  local nsq::T = C[1]
  for i = 2:4:9
    nsq += C[i] * λum^C[i+2] / (λsq - C[i+3]^C[i+4])
  end
  for i = 10:2:N-1
    nsq += C[i]*λum^C[i+1]
  end
  sqrt(nsq)
end

# Cauchy (formula 5)
function refindex{T<:Real}(::Type{Val{Cauchy}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local n::T = C[1]
  for i = 2:2:N-1
    n += C[i]*λum^C[i+1]
  end
  n
end

# Gasses (formula 6)
function refindex{T<:Real}(::Type{Val{Gases}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local n::T = one(T) + C[1]
  local invsqrtλ::T = λum^-2
  for i = 2:2:N-1
    n += C[i] / (C[i+1] - invsqrtλ)
  end
  n
end

# Herzberger (formula 7)
function refindex{T<:Real}(::Type{Val{Herzberger}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert N>=6 "index equation requires 6 coefficients"
  local n::T = C[1]
  local λsq::T = λum^2
  n += C[2] / (λsq - 0.028)
  n += C[3]*(1/(λsq - 0.028))^2
  for i = 4::N
    n += C[i] * λum^(2*(i-3))
  end
  n
end

# Retro (formula 8)
function refindex{T<:Real}(::Type{Val{Retro}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert N>=4 "index equation requires 4 coefficients"
  local λsq::T = λum^2
  local r::T = C[1] + C[2]*λsq/(λsq-C[3]) + C[4]*λsq
  sqrt(-2r - 1) / sqrt(r - 1)
end

# Exotic (formula 9)
function refindex{T<:Real}(::Type{Val{Exotic}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert N>=6 "index equation requires 6 coefficients"
  local λsq::T = λum^2
  local nsq::T = C[1] + C[2]/(λsq-C[3]) + C[4]*(λum - C[5])/((λum - C[5])^2 + C[6])
  sqrt(nsq)
end

function refindex{T<:Real}(m::Material, λum::T)
  if 1 <= m.fnum <= 9 && length(m.C) > 0
    refindex(Val{INDEXFN(m.fnum)}, m.C, λum)
  elseif !isnull(m.Ninterp) && m.range[1] <= λum <= m.range[2]
    get(m.Ninterp)[λum]
  else
    NaN
  end
end

immutable SeriesInterp{T<:Real}
  extrema::Tuple{T, T}
  interp::Interpolations.GriddedInterpolation
  function SeriesInterp(xvals::Array{T,1}, yvals::Array{T,1})
    new(extrema(xvals), interpolate((xvals,), vals, Gridded(Linear())))
  end
  SeriesInterp(table::Array{T, 2}) = SeriesInterp(vec(table[:,1]), vec(table[:,2]))
end

interpolate{T<:Real}(intp::SeriesInterp, x::T) = extrema[1] <= x <= extrema[2] ? intp.interp[x] : NaN

type Material{T<:Real}
  name::AbstractString
  nD::T
  VD::T
  fnum::Int
  extrema::Tuple{T, T}
  C::Array{T,1}
  nintp::Nullable{SeriesInterp{T}}
  kintp::Nullable{SeriesInterp{T}}
  info::Dict{AbstractString,Any}
  comments::AbstractString
  references::AbstractString
  Material(name) = new(name)
end

end # module


text:
# test projections onto intersection for julia for single image desaturation
# as set theoretical estimation (a feasibility problem)
# We will use a VERY simple learning appraoch to obtain 'good' constraints. This
# learing works with just a few even <10 training examples.
using Distributed
@everywhere using SetIntersectionProjection
using MAT
using Interpolations
using Statistics
using LinearAlgebra
using SparseArrays
using Random
using StatsBase

@everywhere mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

#data directory for loading and writing results

#select working precision
FL=32
if     FL==64
  TF = Float64
  TI = Int64
elseif FL==32
  TF = Float32
  TI = Int32
end

#load a very small data set (12 images only) (Mablab files for compatibility with matlab only solvers for comparison...)
file  = matopen(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Data/Ternate_patch.mat"))
mtrue = read(file, "Ternate_patch")
mtrue = convert(Array{TF,3},mtrue)

#split data tensor in a training and evaluation data
m_train      = mtrue[[1:3 ; 5:6 ; 8 ; 10:13 ; 15:20],:,:]
m_evaluation = mtrue[[4, 7, 9, 14],:,:]
m_est        = zeros(TF,size(m_evaluation)) #allocate results

#computational grid for the training images (all images assumed to be on the same grid here)
comp_grid = compgrid((1, 1),(size(m_evaluation,2), size(m_evaluation,3)))
#convert model and computational grid parameters
comp_grid=compgrid( ( convert(TF,comp_grid.d[1]),convert(TF,comp_grid.d[2]) ), comp_grid.n )

#create 'observed data' by creating artificial saturation of images
d_obs = deepcopy(m_evaluation)
d_obs[d_obs.>125.0f0] .= 125.0f0
d_obs[d_obs.<60.0f0]  .= 60.0f0

# train by observing constraints on the data images in training data set
observations = constraint_learning_by_obseration(comp_grid,m_train)

#define a few constraints and what to do with the observations
constraint = Vector{SetIntersectionProjection.set_definitions}()

#bounds:
m_min     = 0.0
m_max     = 255.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#relaxed histogram constraint:
m_min     = observations["hist_min"]
m_max     = observations["hist_max"]
set_type  = "histogram"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#relaxed histogram constraint on discrete derivative of the image:
# m_min     = observations["hist_TV_min"]
# m_max     = observations["hist_TV_max"]
# set_type  = "histogram"
# TD_OP     = "TV"
# app_mode  = ("matrix","")
# custom_TD_OP = ([],false)
# push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

# #rank that preserves 95% of the training images:
# m_min     = 0
# observations["rank_095"]=sort(vec(observations["rank_095"]))
# m_max     = convert(TI,round(quantile(observations["rank_095"],0.25)))
# set_type  = "rank"
# TD_OP     = "identity"
# app_mode  = ("matrix","")
# custom_TD_OP = ([],false)
# push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#nuclear norm constraint
m_min     = 0
m_max     = convert(TF,quantile(vec(observations["nuclear_norm"]),0.25))
set_type  = "nuclear"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#nuclear norm constraint on the x-derivative of the image
m_min     = 0
m_max     = convert(TF,quantile(vec(observations["nuclear_Dx"]),0.25))
set_type  = "nuclear"
TD_OP     = "D_x"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#nuclear norm constraint on the z-derivative of the image
m_min     = 0
m_max     = convert(TF,quantile(vec(observations["nuclear_Dz"]),0.25))
set_type  = "nuclear"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#anisotropic total-variation constraint:
m_min     = 0
m_max     = convert(TF,quantile(vec(observations["TV"]),0.25))
set_type  = "l1"
TD_OP     = "TV"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#l2 constraint on discrete derivatives of the image:
m_min     = 0
m_max     = convert(TF,quantile(vec(observations["D_l2"]),0.25))
set_type  = "l2"
TD_OP     = "TV"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#l1 constraint on DFT coefficients
m_min     = 0
m_max     = convert(TF,quantile(vec(observations["DFT_l1"]),0.50))
set_type  = "l1"
TD_OP     = "DFT"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#bound constraints on x-derivative
m_min     = convert(TF,quantile(vec(observations["D_x_min"]),0.15))
m_max     = convert(TF,quantile(vec(observations["D_x_max"]),0.85))
set_type  = "bounds"
TD_OP     = "D_x"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#bound constraints on z-derivative
m_min     = convert(TF,quantile(vec(observations["D_z_min"]),0.15))
m_max     = convert(TF,quantile(vec(observations["D_z_max"]),0.85))
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#annulus constraint
m_min     = quantile(observations["annulus"],0.15)
m_max     = quantile(observations["annulus"],0.85)
set_type  = "annulus"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#annulus constraint on discrete gradients
m_min     = quantile(observations["TV_annulus"],0.15)
m_max     = quantile(observations["TV_annulus"],0.85)
set_type  = "annulus"
TD_OP     = "TV"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#PARSDMM options:
options=PARSDMM_options()
options=default_PARSDMM_options(options,options.FL)
options.evol_rel_tol = 1f-6
options.feas_tol     = 0.001f0
options.obj_tol      = 0.0002f0
options.adjust_gamma           = true
options.adjust_rho             = true
options.adjust_feasibility_rho = true
options.Blas_active            = true
options.maxit                  = 300
set_zero_subnormals(true)

options.feasibility_only     = false #compute projection of initial guess
options.parallel             = true
options.zero_ini_guess       = false
BLAS.set_num_threads(2)
#FFTW.set_num_threads(2) #can't use this in Julia >0.6

(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL) #obtain projector and transform-domain operator pairs

#add identity matrix as linear operator for the desaturation data-fit constraint
push!(TD_OP,SparseMatrixCSC{TF}(LinearAlgebra.I,comp_grid.n[1]*comp_grid.n[2],comp_grid.n[1]*comp_grid.n[2]))
push!(set_Prop.AtA_offsets,[0]) #these are dummy values, actual ofsetts are automatically detected
push!(set_Prop.ncvx,false)
push!(set_Prop.banded,true)
push!(set_Prop.AtA_diag,true)
push!(set_Prop.dense,false)
push!(set_Prop.tag,("bounds","identity","matrix",""))

#also add a projector onto the data constraint:
#i.e. , or l<=(A*x-m)<=u (or, a norm ||A*x-m||=< sigma)
data = vec(d_obs[1,:,:])
LBD  = data .- 2.0;  LBD=convert(Vector{TF},LBD);
UBD  = data .+ 2.0;  UBD=convert(Vector{TF},UBD);
ind_min_clip = findall(data.==60.0f0)
ind_max_clip = findall(data.==125.0f0)
LBD[ind_min_clip] .= 0.0f0
UBD[ind_max_clip] .= 255.0f0

push!(P_sub,input -> project_bounds!(input,LBD,UBD))

dummy=zeros(TF,prod(comp_grid.n))
(TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

x_ini               = vec(d_obs[1,:,:])
x_ini[ind_max_clip] .= 225.0f0
x_ini[ind_min_clip] .= 0.0f0

options.rho_ini      = ones(TF,length(TD_OP))*1000f0
for i=1:length(options.rho_ini)
  if set_Prop.ncvx[i]==true
    options.rho_ini[i] = 10f0
  end
end

for i=1:size(d_obs,1)
  SNR(in1,in2)=20*log10(norm(in1)/norm(in1-in2))
 [@spawnat pid y[:L][1]=TD_OP[:L][1]*x_ini for pid in y.pids]
 #[@spawnat pid l[:L][1]=randn(TF,length(l[:L][1])) for pid in l.pids]

  p2proj = deepcopy(x_ini) #don't couple initial guess and point to project.
  @time (x,log_PARSDMM) = PARSDMM(p2proj,AtA,TD_OP,set_Prop,P_sub,comp_grid,options,x_ini,[],y)
  m_est[i,:,:]=reshape(x,comp_grid.n)
  println("SNR:", round(SNR(vec(m_evaluation[i,:,:]),vec(m_est[i,:,:])),digits=2))
  println("PSNR:", round(psnr(vec(m_evaluation[i,:,:]),vec(m_est[i,:,:]),maximum(m_evaluation[i,:,:])),digits=2))

  if i+1<=size(d_obs,1)
    data = vec(d_obs[i+1,:,:])
    LBD = data.-2.0;  LBD=convert(Vector{TF},LBD);
    UBD = data.+2.0;  UBD=convert(Vector{TF},UBD);
    ind_min_clip = findall(data.==60.0f0)
    ind_max_clip = findall(data.==125.0f0)
    LBD[ind_min_clip] .= 0.0f0
    UBD[ind_max_clip] .= 255.0f0
    P_sub[end] = x -> project_bounds!(x,LBD,UBD)

    global x_ini = vec(d_obs[i+1,:,:])
    x_ini[ind_max_clip] .= 225.0f0
    x_ini[ind_min_clip] .= 0.0f0
  end

end

SNR(in1,in2)=20*log10(norm(in1)/norm(in1-in2))

ENV["MPLBACKEND"]="qt5agg"
using PyPlot

#all results in one figure
figure(figsize=(4.8, 3.15))
FS  = 5
LS  = 3
TML = 1
PD  = 1
for i=1:size(m_est,1)
  subplot(3,4,i);imshow(d_obs[i,:,:],cmap="gray",vmin=0.0,vmax=255.0); title("Observed",FontSize=FS);tick_params(labelsize=LS,length=TML,pad=PD)
end
for i=1:size(m_est,1)
  subplot(3,4,i+8);imshow(m_est[i,:,:],cmap="gray",vmin=0.0,vmax=255.0); title(string("PARSDMM, PSNR=", round(psnr(vec(m_evaluation[i,:,:]),vec(m_est[i,:,:]),maximum(m_evaluation[i,:,:])),digits=2)),FontSize=FS);;tick_params(labelsize=LS,length=TML,pad=PD)
end
for i=1:size(m_est,1)
  subplot(3,4,i+4);imshow(m_evaluation[i,:,:],cmap="gray",vmin=0.0,vmax=255.0); title("True",FontSize=FS);tick_params(labelsize=LS,length=TML,pad=PD)
end
tight_layout()
PyPlot.subplots_adjust(wspace=-0.5, hspace=0.3)
savefig("desaturation_results.png",bbox_inches="tight",dpi=600)
close("all")

#First 8 in 1 figure
figure(figsize=(4.8, 2.4))
for i=1:8
  subplot(2,4,i);imshow(m_train[i,:,:],cmap="gray",vmin=0.0,vmax=255.0);axis("off") #title("training image", fontsize=10)
end
tight_layout()
PyPlot.subplots_adjust(wspace=-0.5, hspace=0.1)
savefig("training_data_first8.png",bbox_inches="tight",dpi=600)
close()

#plot training images
figure();title("training image", fontsize=10)
for i=1:16
  subplot(4,4,i);imshow(m_train[i,:,:],cmap="gray",vmin=0.0,vmax=255.0);axis("off") #title("training image", fontsize=10)
end
savefig("training_data_all.eps",bbox_inches="tight",dpi=600)
close()

for i=1:16
  figure();title(string("training image", i), fontsize=10)
  imshow(m_train[i,:,:],cmap="gray",vmin=0.0,vmax=255.0);axis("off") #title("training image", fontsize=10)
  savefig(string("training_data_", i,".eps"),bbox_inches="tight",dpi=600)
  close()
end

#plot results individually
for i=1:size(m_est,1)
    figure();imshow(d_obs[i,:,:],cmap="gray",vmin=0.0,vmax=255.0); title("observed");
    savefig(string("saturized_observed",i,".eps"),bbox_inches="tight",dpi=600)
    figure();imshow(m_est[i,:,:],cmap="gray",vmin=0.0,vmax=255.0); title(string("PARSDMM, SNR=", round(SNR(vec(m_evaluation[i,:,:]),vec(m_est[i,:,:])),digits=2)))
    savefig(string("PARSDMM_desaturation",i,".eps"),bbox_inches="tight",dpi=600)
    figure();imshow(m_evaluation[i,:,:],cmap="gray",vmin=0.0,vmax=255.0); title("True")
    savefig(string("desaturation_evaluation",i,".eps"),bbox_inches="tight",dpi=600)
    close("all")
end

#save stuff
file = matopen("m_est.mat", "w")
write(file, "m_est", convert(Array{Float64,3},m_est))
close(file)

file = matopen("m_evaluation.mat", "w")
write(file, "m_evaluation", convert(Array{Float64,3},m_evaluation))
close(file)

file = matopen("m_train.mat", "w")
write(file, "m_train", convert(Array{Float64,3},m_train))
close(file)

file = matopen("d_obs.mat", "w")
write(file, "d_obs", convert(Array{Float64,3},d_obs))
close(file)

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
include("cfd.jl")
include("structs.jl")
include("thermo.jl")
include("dump.jl")
include("plot.jl")
include("test.jl")
include("SimulationParameters.jl")

## https://stackoverflow.com/questions/37200025/how-to-import-custom-module-in-julia
using DelimitedFiles, ..SimulationParameters
sp = SimulationParameters

function update!(syswork::System, sysnext::System, p::Parameters,
    t::Integer)
# -------------------- Predictor ------------------------- #
## ρ̃ᵢⱼnext = ρᵢⱼwork + F⋅Δt
    sysnext_pred = copy(sysnext)

    continuity_equation!(syswork, sysnext_pred, p)
    find_pressure!(sysnext_pred, p)
    apply_bc_pres_sat_den!(sysnext_pred, p, t)
    darcy!(sysnext_pred, p)
    apply_bc_velocity!(sysnext_pred, p)
# -------------------------------------------------------- #

# -------------------- Corrector ------------------------- #
## ρᵢⱼnext = ρᵢⱼwork + (F + F̃) / 2 ⋅Δt
    continuity_equation!(syswork, sysnext_pred, sysnext, p)
    find_pressure!(sysnext, p)
    apply_bc_pres_sat_den!(sysnext, p, t)
    darcy!(sysnext, p)
    apply_bc_velocity!(sysnext, p)
# -------------------------------------------------------- #
end

function init(p::Parameters)
    check_bc(p.bc)
    P = fill(p.P₀, p.nx, p.ny)
    s = zeros(p.nx, p.ny, 2)
    u = zeros(p.nx, p.ny, 2)
    v = zeros(p.nx, p.ny, 2)
    ρ = zeros(p.nx, p.ny, 2)
    F = zeros(p.nx, p.ny, 2)
    bc_RHS = initial_RHS(p)

    sys = System(u, v, ρ, F, s, P, bc_RHS)

    initial_conditions!(sys, p)

    return sys
end

function filtration(p::Parameters)
    syswork = init(p)
    sysnext = copy(syswork)

    for t = 1 : p.nsteps
        check_saturation(syswork, p)
        check_density(syswork, p)
        check_pressure(syswork)
        check_continuity(syswork, sysnext, p, t)
        println("Step i = ", t)
        update!(syswork, sysnext, p, t)
        sysnext, syswork = syswork, sysnext
        # dump(syswork, p)
        dump_total_liquid_flux(syswork, p, t)
    end
    return syswork
end

function main(p)
    sys = filtration(p)
    plot(sys, p)
    dump(sys, p)
end

## -------------------------- BC -------------------------- ##

a_u = zeros(max(nx, ny), 4, 2)
a_v = zeros(max(nx, ny), 4, 2)
a_ρ = zeros(max(nx, ny), 4, 2)
a_s = zeros(max(nx, ny), 4, 2)
a_P = zeros(max(nx, ny), 4, 2)

b_u = zeros(max(nx, ny), 4, 2)
b_v = zeros(max(nx, ny), 4, 2)
b_ρ = zeros(max(nx, ny), 4, 2)
b_s = zeros(max(nx, ny), 4, 2)
b_P = zeros(max(nx, ny), 4, 2)

## h[:][1][k] -> left
##
## h[:][2][k] -> right
##
## h[:][3][k] -> top
##
## h[:][4][k] -> bottom

k = 1 : 2

## top & bottom
@views a_v[:, :, k]     .= 1

## left & right
@views a_u[:, :, k]     .= 1

## left & right & top & bottom
@views a_ρ[:, :, k]     .= 1

## inlet at bottom [1, nx ÷ 2]
@views a_s[inlet, 4, k] .= 1

## walls & outlet
@views b_s[:, 1:3,  k]  .= 1
@views b_s[wall, 4, k]  .= 1

## outlet
@views a_P[:, 3, k]     .= 1

## inlet
@views a_P[inlet, 4, k] .= 1

## walls left & right
@views b_P[:, 1:2, k]   .= 1

## bottom wall
@views b_P[wall, 4, k]  .= 1

u_bc = BoundaryCondition(a_u, b_u)
v_bc = BoundaryCondition(a_v, b_v)
ρ_bc = BoundaryCondition(a_ρ, b_ρ)
s_bc = BoundaryCondition(a_s, b_s)
P_bc = BoundaryCondition(a_P, b_P)

bc = BoundaryConditions(u_bc, v_bc, ρ_bc, s_bc, P_bc)

## Problem!!
## https://discourse.julialang.org/t/writing-functions-for-types-defined-in-another-module/31895/4
## https://discourse.julialang.org/t/structs-in-modules/34317
println("TODO: fix")
println("typeof(sp.bc): ", typeof(sp.bc))
println("typeof(bc): ", typeof(bc))

p = Parameters(sp.nsteps, sp.nx, sp.ny, sp.L, sp.Δx,
               sp.Δy, sp.Δt, sp.φ, sp.K, sp.Pin,
               sp.Pout, sp.P₀, sp.ψ, sp.ψ₀, sp.M, sp.μ,
               bc, sp.wall, sp.inlet, sp.relax_steps)

main(p)

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
using LoewnerJohnEllipsoids
using Base.Test

# write your own tests here
function test_unit_circle()
	polyhedron = box([-1.,-1], [1., 1])
	e = inner_ellipsoid(polyhedron)
	@show e
	@test isapprox(e.d, [0,0], atol=1e-4)
	@test isapprox(e.C, [1 0; 0 1], atol=1e-4)
end
test_unit_circle()

function test_shifted_scaled_box()
	polyhedron = box([1,1], [2,2])
	e = inner_ellipsoid(polyhedron)
	@show e
	@test isapprox(e.d, [1.5, 1.5], atol=1e-4)
	@test isapprox(e.C, [0.5 0; 0 0.5], atol=1e-4)
end
test_shifted_scaled_box()

function test_outer_unit_box()
	points = Any[[-1,-1], [-1,1], [1,-1], [1,1]]
	e = outer_ellipsoid(points)
	@show e
	@test isapprox(e.d, [0, 0], atol=1e-4)
	@test isapprox(e.C, [sqrt(2) 0; 0 sqrt(2)], atol=1e-4)
end
test_outer_unit_box()

function test_barrier()
	polyhedron = box([1,1], [2,2])
	e = inner_ellipsoid(polyhedron)
	@test isapprox(barrier_value(e, [1,1.5]), 0, atol=1e-4)
	@test isapprox(barrier_value(e, [2,1.5]), 0, atol=1e-4)
	@test barrier_value(e, [1.5, 1.5]) < 0
	@test barrier_value(e, [2.5, 2.5]) > 0
end
test_barrier()


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# 2-arg show: show schema and list # of metadata entries if non-zero
function Base.show(io::IO, table::Table)
    ncols = length(Tables.columnnames(table))
    print(io, "$(typeof(table)) with $(Tables.rowcount(table)) rows, $(ncols) columns,")
    meta = getmetadata(table)
    if meta !== nothing && !isempty(meta)
        print(io, " ", length(meta), " metadata entries,")
    end
    sch = Tables.schema(table)
    print(io, " and schema:\n")
    show(IOContext(io, :print_schema_header => false), sch)
    return nothing
end

# 3-arg show: show schema and show metadata entries adaptively according to `displaysize`
function Base.show(io::IO, mime::MIME"text/plain", table::Table)
    display_rows, display_cols = displaysize(io)
    ncols = length(Tables.columnnames(table))
    meta = getmetadata(table)
    if meta !== nothing
        display_rows -= 1 # decrement for metadata header line
        display_rows -= min(length(meta), 2) # decrement so we can show at least 2 lines of metadata
    end
    print(io, "$(typeof(table)) with $(Tables.rowcount(table)) rows, $(ncols) columns, and ")
    sch = Tables.schema(table)
    print(io, "schema:\n")
    schema_context = IOContext(io, :print_schema_header => false, :displaysize => (max(display_rows, 3), display_cols))
    schema_str = sprint(show, mime, sch; context=schema_context)
    print(io, schema_str)
    display_rows -= (count("\n", schema_str) + 1) # decrement for number of lines printed
    if meta !== nothing
        print(io, "\n\nwith metadata given by a ")
        show(IOContext(io, :displaysize => (max(display_rows, 5), display_cols)), mime, meta)
    end
    return nothing
end

max_stars_repo_name: jonalm/arrow-julia

max_issues_repo_name: jonalm/arrow-julia

max_forks_repo_name: jonalm/arrow-julia


text:
#using GaussianProcesses
import GaussianProcesses: Base.size
function useme()
	include("import.jl")
end
#useme()
println("Hello")

max_stars_repo_name: krystophny/profit

max_issues_repo_name: krystophny/profit

max_forks_repo_name: krystophny/profit


text:
using PeaksAndValleys
using Base.Test

using Compat

function test_initial_pivot()
    # Strictly increasing.
    data = linspace(1.0, 2.0, 10)
    @test identify_initial_pivot(data, 0.1, -0.1) == VALLEY

    # Strictly decreasing.
    data = linspace(1, 0.5, 10)
    @test identify_initial_pivot(data, 0.1, -0.1) == PEAK

    # Increasing but kinked.
    data = [1.0, 0.99, 1.1]
    @test identify_initial_pivot(data, 0.1, -0.1) == PEAK

    # Decreasing but kinked.
    data = [1.0, 1.01, 0.9]
    @test identify_initial_pivot(data, 0.1, -0.1) == VALLEY

    # Strictly increasing but under threshold.
    data = linspace(1, 1.01, 10)
    @test identify_initial_pivot(data, 0.1, -0.1) == VALLEY

    # Strictly decreasing but under threshold.
    data = linspace(1, 0.99, 10)
    @test identify_initial_pivot(data, 0.1, -0.1) == PEAK

    # Increasing but under threshold and kinked.
    data = [1.0, 0.99, 1.02]
    @test identify_initial_pivot(data, 0.1, -0.1) == VALLEY

    # Decreasing but under threshold and kinked.
    data = [1.0, 1.01, 0.99]
    @test identify_initial_pivot(data, 0.1, -0.1) == PEAK
end

function test_pivots()
    # The relative change should be a negative number. When coding,
    # I kept mistakenly using positive numbers for both. Having
    # a positive and a negative makes things more clear, without
    # having to use named parameters.
    data = [1.0, 2.0, 3.0]
    @test_throws AssertionError peak_valley_pivots(data, 0.1, 0.1)

    # Strictly increasing.
    data = linspace(1, 10, 10)
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = zeros(data)
    expected_result[1], expected_result[end] = VALLEY, PEAK
    @test result == expected_result

    # Strictly decreasing.
    data = linspace(10, 0, 10)
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = zeros(data)
    expected_result[1], expected_result[end] = PEAK, VALLEY
    @test result == expected_result

    # Strictly increasing but under threshold.
    data = linspace(1.0, 1.05, 10)
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = zeros(data)
    expected_result[1], expected_result[end] = VALLEY, PEAK
    @test data[1] < data[end]
    @test result == expected_result

    # Strictly decreasing but under threshold.
    data = linspace(1.05, 1.0, 10)
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = zeros(data)
    expected_result[1], expected_result[end] = PEAK, VALLEY
    @test result == expected_result

    # Increasing but kinked.
    data = [1.0, 0.99, 1.1]
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = [PEAK, VALLEY, PEAK]
    @test result == expected_result

    # Decreasing but kinked.
    data = [1.0, 1.01, 0.9]
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = [VALLEY, PEAK, VALLEY]
    @test result == expected_result

    # Single peaked.
    data = [1.0, 1.2, 1.05]
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = [VALLEY, PEAK, VALLEY]
    @test result == expected_result

    # Single valleyed.
    data = [1.0, 0.9, 1.2]
    result = peak_valley_pivots(data, 0.1, -0.1)
    expected_result = [PEAK, VALLEY, PEAK]
    @test result == expected_result
end

function test_segment_return()
    # Strictly increasing.
    data = linspace(1.0, 100.0, 10)
    pivots = peak_valley_pivots(data, 0.1, -0.1)
    @test compute_segment_returns(data, pivots) == [99.0]

    # Strictly decreasing.
    data = linspace(100.0, 1.0, 10)
    pivots = peak_valley_pivots(data, 0.1, -0.1)
    @test compute_segment_returns(data, pivots) == [-0.99]

    # Rise-Fall-Rise.
    data = [1.0, 1.05, 1.1, 1.0, 0.9, 1.5]
    pivots = peak_valley_pivots(data, 0.1, -0.1)
    result = compute_segment_returns(data, pivots)
    expected = [0.1, -0.181818, 0.6666666]
    for (a,b) in zip(result, expected)
        @test_approx_eq_eps a b 1e-2
    end
end

function test_max_drawdown()
    # Strictly increasing.
    data = linspace(1.0, 100.0, 10)
    @test max_drawdown(data) == 0.0

    # Strictly decreasing.
    data = linspace(100.0, 1.0, 10)
    @test_approx_eq max_drawdown(data) 0.99

    # Rise-Fall-Rise.
    data = [1.0, 1.05, 1.1, 1.0, 0.9, 1.5]
    @test_approx_eq max_drawdown(data) 0.18181818181818188
end

function test_modes_to_pivots()
    data = [1, 0, 0, 0, -1, 0, 0, 1, -1, 0, 1]
    result = pivots_to_modes(data)
    expected = [1, -1, -1, -1, -1, 1, 1, 1, -1, 1, 1]

    @test result == expected
end

test_initial_pivot()
test_pivots()
test_segment_return()
test_max_drawdown()
test_modes_to_pivots()


text:
using POMDPModelTools
using POMDPs
using POMDPModels
using Random
using Test
using Pkg
using POMDPSimulators
using BeliefUpdaters
using POMDPPolicies
import Distributions.Categorical
using SparseArrays
using DiscreteValueIteration
import CommonRLInterface

using Documenter: doctest

doctest(POMDPModelTools)

@testset "POMDPModelTools" begin
    @testset "ordered" begin
        include("test_ordered_spaces.jl")
    end

    # require POMDPModels
    @testset "genbeliefmdp" begin
        include("test_generative_belief_mdp.jl")
    end
    @testset "implement" begin
        include("test_implementations.jl")
    end
    @testset "weightediter" begin
        include("test_weighted_iteration.jl")
    end
    @testset "sparsecat" begin
        include("test_sparse_cat.jl")
    end
    @testset "bool" begin
        include("test_bool.jl")
    end
    @testset "deterministic" begin
        include("test_deterministic.jl")
    end
    @testset "uniform" begin
        include("test_uniform.jl")
    end
    @testset "implicit" begin
        include("test_implicit.jl")
    end
    @testset "terminalstate" begin
        include("test_terminal_state.jl")
    end

    # require POMDPModels
    @testset "info" begin
        include("test_info.jl")
    end
    @testset "obsweight" begin
        include("test_obs_weight.jl")
    end

    @testset "visolve" begin
        include("test_fully_observable_pomdp.jl")
        include("test_underlying_mdp.jl")
    end

    @testset "vis" begin
        include("test_visualization.jl")
    end

    @testset "evaluation" begin
        include("test_evaluation.jl")
    end

    @testset "pretty printing" begin
        include("test_pretty_printing.jl")
    end

    @testset "sparse tabular" begin
        include("test_tabular.jl")
    end

    @testset "reward model" begin
        include("test_reward_model.jl")
    end

    @testset "matrices" begin
        include("test_matrices.jl")
    end

    @testset "Common RL" begin
        include("test_common_rl.jl")
    end
end


text:
# Quadratic Costs
n,m = 3,2
Q = Diagonal(@SVector fill(0.1, n))
R = Diagonal(@SVector fill(0.01, m))
Qf = Diagonal(@SVector fill(10., n))
xf = @SVector ones(n)

quadcost = QuadraticCost(Q,R)
@test quadcost.Q == Q
@test quadcost.q == zeros(n)
@test quadcost.R == R
@test quadcost.r == zeros(m)
@test quadcost.c == 0
@test state_dim(quadcost) == n
@test control_dim(quadcost) == m

quadcost = LQRCost(Q,R,xf)
@test quadcost.Q == Q
@test quadcost.q == -Q*xf
@test quadcost.R == R
@test quadcost.r == zeros(m)
@test quadcost.c == 0.5*xf'Q*xf

quadcost = TO.LQRCostTerminal(Q,xf)
@test quadcost.Q == Q
@test quadcost.q == -Q*xf
@test isempty(quadcost.R)
@test isempty(quadcost.r)
@test quadcost.c == 0.5*xf'Q*xf

q = @SVector rand(n)
quadcost = QuadraticCost(Q,q,0.5)
@test quadcost.Q == Q
@test quadcost.q == q
@test isempty(quadcost.R)
@test isempty(quadcost.r)
@test quadcost.c == 0.5

Q2 = Diagonal(@SVector [-1, 0.1, 0.2])
R2 = Diagonal(@SVector [0, 0.1])
@test_throws ArgumentError QuadraticCost(Q2,R)
@test_logs (:warn, "R is not positive definite") QuadraticCost(Q,R2)


text:
#################################
"""
    converttodatetime!(fgmax::VisClaw.FGmax, t0::Dates.DateTime)
    converttodatetime!(gauge::VisClaw.Gauge, unit::Dates.DateTime)
    converttodatetime!(amrs::VisClaw.AMR, unit::Dates.DateTime)
    converttodatetime!(track::VisClaw.Track, unit::Dates.DateTime)

Time unit converter to Dates.DateTime
"""
function converttodatetime!(fgmax::VisClaw.FGmax, t0::Dates.DateTime)
    ## return if already converted
    isa(fgmax.unittime, Dates.DateTime) && (return fgmax)

    ## factor
    ratio = timedict[fgmax.unittime]

    ## (temporal) to avoid NaN convert error
    fgmax.tD[isnan.(fgmax.tD)] .= 0.0
    fgmax.tarrival[isnan.(fgmax.tarrival)] .= 0.0
    fgmax.tv[isnan.(fgmax.tv)] .= 0.0
    fgmax.tM[isnan.(fgmax.tM)] .= 0.0
    fgmax.tMflux[isnan.(fgmax.tMflux)] .= 0.0
    fgmax.tDmin[isnan.(fgmax.tDmin)] .= 0.0

    ## convert
    fgmax.unittime = :DateTime
    fgmax.tD = @. t0 + Dates.Second(round(ratio*fgmax.tD))
    fgmax.tarrival = @. t0 + Dates.Second(round(ratio*fgmax.tarrival))
    if !isempty(fgmax.tv)
        fgmax.tv = @. t0 + Dates.Second(round(ratio*fgmax.tv))
    end
    if !isempty(fgmax.tM)
        fgmax.tM = @. t0 + Dates.Second(round(ratio*fgmax.tM))
        fgmax.tMflux = @. t0 + Dates.Second(round(ratio*fgmax.tMflux))
        fgmax.tDmin = @. t0 + Dates.Second(round(ratio*fgmax.tDmin))
    end
    return fgmax
end
#################################

#################################
function converttodatetime!(gauge::VisClaw.Gauge, t0::Dates.DateTime)
    ## return if already converted
    isa(gauge.unittime, Dates.DateTime) && (return gauge)

    ## factor
    ratio = timedict[gauge.unittime]

    ## convert
    gauge.unittime = :DateTime
    gauge.time = @. t0 + Dates.Second(round(ratio*gauge.time))
    # return value
    return gauge
end
#################################

#################################
function converttodatetime!(amrs::VisClaw.AMR, t0::Dates.DateTime)
    ## return if already converted
    isa(amrs.unittime, Dates.DateTime) && (return amrs)

    ## factor
    ratio = timedict[amrs.unittime]

    ## convert
    amrs.unittime = :DateTime
    amrs.timelap = @. t0 + Dates.Second(round(ratio*amrs.timelap))
    # return value
    return amrs
end
#################################

#################################
function converttodatetime!(track::VisClaw.Track, t0::Dates.DateTime)
    ## return if already converted
    isa(track.unittime, Dates.DateTime) && (return track)

    ## factor
    ratio = timedict[track.unittime]

    ## convert
    track.unittime = :DateTime
    track.timelap = @. t0 + Dates.Second(round(ratio*track.timelap))
    # return value
    return track
end
#################################


text:
using NeuralPDE
using Quadrature, Cubature, Cuba
using Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots
using DelimitedFiles
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum

function hamilton_jacobi(strategy, minimizer, maxIters)

    ##  DECLARATIONS
    @parameters  t x1 x2 x3 x4
    @variables   u(..)

    Dt = Differential(t)

    Dx1 = Differential(x1)
    Dx2 = Differential(x2)
    Dx3 = Differential(x3)
    Dx4 = Differential(x4)

    Dxx1 = Differential(x1)^2
    Dxx2 = Differential(x2)^2
    Dxx3 = Differential(x3)^2
    Dxx4 = Differential(x4)^2

    # Discretization
    tmax         = 1.0
    x1width      = 1.0
    x2width      = 1.0
    x3width      = 1.0
    x4width      = 1.0

    tMeshNum     = 10
    x1MeshNum    = 10
    x2MeshNum    = 10
    x3MeshNum    = 10
    x4MeshNum    = 10

    dt   = tmax/tMeshNum
    dx1  = x1width/x1MeshNum
    dx2  = x2width/x2MeshNum
    dx3  = x3width/x3MeshNum
    dx4  = x4width/x4MeshNum

    domains = [t ∈ Interval(0.0,tmax),
               x1 ∈ Interval(0.0,x1width),
               x2 ∈ Interval(0.0,x2width),
               x3 ∈ Interval(0.0,x3width),
               x4 ∈ Interval(0.0,x4width)]

    ts  = 0.0 : dt : tmax
    x1s = 0.0 : dx1 : x1width
    x2s = 0.0 : dx2 : x2width
    x3s = 0.0 : dx3 : x3width
    x4s = 0.0 : dx4 : x4width

    λ = 1.0f0

    # Operators
    Δu = Dxx1(u(t,x1,x2,x3,x4)) + Dxx2(u(t,x1,x2,x3,x4)) + Dxx3(u(t,x1,x2,x3,x4)) + Dxx4(u(t,x1,x2,x3,x4)) # Laplacian
    ∇u = [Dx1(u(t,x1,x2,x3,x4)), Dx2(u(t,x1,x2,x3,x4)),Dx3(u(t,x1,x2,x3,x4)),Dx4(u(t,x1,x2,x3,x4))]

    # Equation
    eq = Dt(u(t,x1,x2,x3,x4)) + Δu - λ*sum(∇u.^2) ~ 0  #HAMILTON-JACOBI-BELLMAN EQUATION

    terminalCondition =  log((1 + x1*x1 + x2*x2 + x3*x3 + x4*x4)/2) # see PNAS paper

    bcs = [u(tmax,x1,x2,x3,x4) ~ terminalCondition]  #PNAS paper again

    ## NEURAL NETWORK
    n = 10   #neuron number

    chain = FastChain(FastDense(5,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))   #Neural network from Flux library

    indvars = [t,x1,x2,x3,x4]   #phisically independent variables
    depvars = [u(t,x1,x2,x3,x4)]       #dependent (target) variable

    dim = length(domains)

    losses = []
    error  = []
    times = []

    dx_err = 0.2

    error_strategy = GridTraining(dx_err)

    initθ = Float64.(DiffEqFlux.initial_params(chain))
    eltypeθ = eltype(initθ)
    parameterless_type_θ = DiffEqBase.parameterless_type(initθ)

    phi = NeuralPDE.get_phi(chain,parameterless_type_θ)
    derivative = NeuralPDE.get_numeric_derivative()

    _pde_loss_function = NeuralPDE.build_loss_function(eq,indvars,depvars,
                                             phi,derivative,nothing,chain,initθ,error_strategy)

    bc_indvars = NeuralPDE.get_variables(bcs,indvars,depvars)
    _bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars,
                                              phi,derivative,nothing,chain,initθ,error_strategy,
                                              bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]

    train_sets = NeuralPDE.generate_training_sets(domains,dx_err,[eq],bcs,eltypeθ,indvars,depvars)
    train_domain_set, train_bound_set = train_sets

    pde_loss_functions = [NeuralPDE.get_loss_function(_pde_loss_function,
                                          train_domain_set[1],eltypeθ,
                                          parameterless_type_θ,error_strategy)]

    bc_loss_functions = [NeuralPDE.get_loss_function(bc,
                                          train_set,eltypeθ,
                                          parameterless_type_θ,error_strategy) for (train_set,bc) in zip(train_bound_set, _bc_loss_functions)]
    loss_functions =  [pde_loss_functions;bc_loss_functions]
    function loss_function_(θ,p)
        return sum(map(l->l(θ) ,loss_functions))
    end

    cb_ = function (p,l)
        deltaT_s = time_ns() #Start a clock when the callback begins, this will evaluate questo misurerà anche il calcolo degli uniform error

        ctime = time_ns() - startTime - timeCounter #This variable is the time to use for the time benchmark plot
        append!(times, ctime/10^9) #Conversion nanosec to seconds
        append!(losses, l)
        loss_ = loss_function_(p,nothing)
        append!(error, loss_)

        timeCounter = timeCounter + time_ns() - deltaT_s #timeCounter sums all delays due to the callback functions of the previous iterations

        return false
    end

    @named pde_system = PDESystem(eq, bcs, domains, indvars, depvars)

    discretization = NeuralPDE.PhysicsInformedNN(chain, strategy; init_params =initθ)
    prob = NeuralPDE.discretize(pde_system,discretization)

    timeCounter = 0.0
    startTime = time_ns() #Fix initial time (t=0) before starting the training
    res = GalacticOptim.solve(prob, minimizer, cb = cb_, maxiters=maxIters)

    phi = discretization.phi

    params = res.minimizer

    # Model prediction
    domain = [ts,x1s,x2s,x3s,x4s]

    u_predict = [reshape([first(phi([t,x1,x2,x3,x4],res.minimizer)) for x1 in x1s for x2 in x2s for x3 in x3s for x4 in x4s], (length(x1s),length(x2s), length(x3s),length(x4s))) for t in ts]  #matrix of model's prediction

    return [error, params, domain, times, losses]
end

maxIters = [(1,1,1,1000,1000,1000,1000),(1,1,1,300,300,300,300)] #iters for ADAM/LBFGS
# maxIters = [(1,1,1,1,1,2,2),(1,1,1,3,3,3,3)] #iters for ADAM/LBFGS

strategies = [NeuralPDE.QuadratureTraining(quadrature_alg = CubaCuhre(), reltol = 1, abstol = 1e-4, maxiters = 100),
              NeuralPDE.QuadratureTraining(quadrature_alg = HCubatureJL(), reltol = 1, abstol = 1e-4, maxiters = 100, batch = 0),
              NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh(), reltol = 1, abstol = 1e-4, maxiters = 100),
              NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLp(), reltol = 1, abstol = 1e-4, maxiters = 100),
              NeuralPDE.GridTraining(0.2),
              NeuralPDE.StochasticTraining(400 ; bcs_points= 50),
              NeuralPDE.QuasiRandomTraining(400 ; bcs_points= 50)]

strategies_short_name = ["CubaCuhre",
                        "HCubatureJL",
                        "CubatureJLh",
                        "CubatureJLp",
                        "GridTraining",
                        "StochasticTraining",
                        "QuasiRandomTraining"]

minimizers = [ADAM(0.005),
              #BFGS()]
              LBFGS()]

minimizers_short_name = ["ADAM",
                         "LBFGS"]
                        #"BFGS"]

# Run models
error_res =  Dict()
domains = Dict()
params_res = Dict()  #to use same params for the next run
times = Dict()
losses_res = Dict()

print("Starting run")
## Convergence

for min =1:length(minimizers) # minimizer
      for strat=1:length(strategies) # strategy
            # println(string(strategies_short_name[strat], "  ", minimizers_short_name[min]))
            res = hamilton_jacobi(strategies[strat], minimizers[min], maxIters[min][strat])
            push!(error_res, string(strat,min)     => res[1])
            push!(params_res, string(strat,min) => res[2])
            push!(domains, string(strat,min)        => res[3])
            push!(times, string(strat,min)        => res[4])
            push!(losses_res, string(strat,min)   => res[5])
      end
end

#Plotting the first strategy with the first minimizer out from the loop to initialize the canvas
current_label = string(strategies_short_name[1], " + " , minimizers_short_name[1])
error = Plots.plot(times["11"], error_res["11"], yaxis=:log10, label = current_label)#, xlims = (0,10))#legend = true)#, size=(1200,700))
plot!(error, times["21"], error_res["21"], yaxis=:log10, label = string(strategies_short_name[2], " + " , minimizers_short_name[1]))
plot!(error, times["31"], error_res["31"], yaxis=:log10, label = string(strategies_short_name[3], " + " , minimizers_short_name[1]))
plot!(error, times["41"], error_res["41"], yaxis=:log10, label = string(strategies_short_name[4], " + " , minimizers_short_name[1]))
plot!(error, times["51"], error_res["51"], yaxis=:log10, label = string(strategies_short_name[5], " + " , minimizers_short_name[1]))
plot!(error, times["61"], error_res["61"], yaxis=:log10, label = string(strategies_short_name[6], " + " , minimizers_short_name[1]))
plot!(error, times["71"], error_res["71"], yaxis=:log10, label = string(strategies_short_name[7], " + " , minimizers_short_name[1]))

plot!(error, times["12"], error_res["12"], yaxis=:log10, label = string(strategies_short_name[1], " + " , minimizers_short_name[2]))
plot!(error, times["22"], error_res["22"], yaxis=:log10, label = string(strategies_short_name[2], " + " , minimizers_short_name[2]))
plot!(error, times["32"], error_res["32"], yaxis=:log10, label = string(strategies_short_name[3], " + " , minimizers_short_name[2]))
plot!(error, times["42"], error_res["42"], yaxis=:log10, label = string(strategies_short_name[4], " + " , minimizers_short_name[2]))
plot!(error, times["52"], error_res["52"], yaxis=:log10, label = string(strategies_short_name[5], " + " , minimizers_short_name[2]))
plot!(error, times["62"], error_res["62"], yaxis=:log10, label = string(strategies_short_name[6], " + " , minimizers_short_name[2]))
plot!(error, times["72"], error_res["72"], yaxis=:log10, title = string("Hamilton Jacobi convergence ADAM/LBFGS"), ylabel = "log(error)",xlabel = "t", label = string(strategies_short_name[7], " + " , minimizers_short_name[2]))

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput