text:
#
# oiformat1.jl --
#
# Define 1st revision of OI-FITS format.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015: Éric Thiébaut.
#
#------------------------------------------------------------------------------

# OI-FITS FORMAT DESCRIPTION TABLES
#
# The format of the OI-FITS data block is described by a vector of strings
# like:
#
#   ["KEYWORD FORMAT UNITS DESCR",
#     ...,
#     ...,
#    "---------------------------",
#    "COLUMN  FORMAT UNITS DESCR",
#     ...,
#     ...,
#     ...]
#
# where:
#
#   KEYWORD = keyword for HDU header.
#   COLUMN = column name for table (TTYPE).
#   FORMAT = nL where n is an integer and L a letter;
#       for keywords, 0 means optional and 1 means required;
#       for columns, a negative number means abs(n)*NWAVE.
#   UNITS = default units.
#   DESCR = description/comment.
#
# There may be any number of keyword definitions and any number of column
# definitions, the two parts are separated by a dash line like
# "--------------".
#

# OI_TARGET definition (1st revision):
add_def("OI_TARGET", 1,
        ["OI_REVN    1I -      revision number of the table definition",
         "------------------------------------------------------------",
         "TARGET_ID  1I -      index number",
         "TARGET    16A -      target name",
         "RAEP0      1D deg    RA at mean equinox",
         "DECEP0     1D deg    DEC at mean equinox",
         "EQUINOX    1E yr     equinox",
         "RA_ERR     1D deg    error in RA at mean equinox",
         "DEC_ERR    1D deg    error in DEC at mean equino",
         "SYSVEL     1D m/s    systemic radial velocity",
         "VELTYP     8A -      reference for radial velocity",
         "VELDEF     8A -      definition of radial velocity",
         "PMRA       1D deg/yr proper motion in RA",
         "PMDEC      1D deg/yr proper motion in DEC",
         "PMRA_ERR   1D deg/yr error of proper motion in RA",
         "PMDEC_ERR  1D deg/yr error of proper motion in DEC",
         "PARALLAX   1E deg    parallax",
         "PARA_ERR   1E deg    error in parallax",
         "SPECTYP   16A -      spectral type"])

# OI_ARRAY definition (1st revision):
add_def("OI_ARRAY", 1,
        ["OI_REVN    1I - revision number of the table definition",
         "ARRNAME    1A - array name for cross-referencing",
         "FRAME      1A - coordinate frame",
         "ARRAYX     1D m array center X-coordinate",
         "ARRAYY     1D m array center Y-coordinate",
         "ARRAYZ     1D m array center Z-coordinate",
         "------------------------------------------------------------",
         "TEL_NAME  16A - telescope name",
         "STA_NAME  16A - station name",
         "STA_INDEX  1I - station index",
         "DIAMETER   1E m element diameter",
         "STAXYZ     3D m station coordinates relative to array center"])

# OI_WAVELENGTH definition (1st revision):
add_def("OI_WAVELENGTH", 1,
        ["OI_REVN    1I - revision number of the table definition",
         "INSNAME    1A - name of detector for cross-referencing",
         "------------------------------------------------------------",
         "EFF_WAVE   1E m effective wavelength of channel",
         "EFF_BAND   1E m effective bandpass of channel"])

# OI_VIS definition (1st revision):
add_def("OI_VIS", 1,
        ["OI_REVN    1I -   revision number of the table definition",
         "DATE-OBS   1A -   UTC start date of observations",
         "ARRNAME    0A -   name of corresponding array",
         "INSNAME    1A -   name of corresponding detector",
         "------------------------------------------------------------",
         "TARGET_ID  1I -   target number as index into OI_TARGET table",
         "TIME       1D s   UTC time of observation",
         "MJD        1D day modified Julian Day",
         "INT_TIME   1D s   integration time",
         "VISAMP    -1D -   visibility amplitude",
         "VISAMPERR -1D -   error in visibility amplitude",
         "VISPHI    -1D deg visibility phase",
         "VISPHIERR -1D deg error in visibility phase",
         "UCOORD     1D m   U coordinate of the data",
         "VCOORD     1D m   V coordinate of the data",
         "STA_INDEX  2I -   station numbers contributing to the data",
         "FLAG      -1L -   flag"])

# OI_VIS2 definition (1st revision):
add_def("OI_VIS2", 1,
        ["OI_REVN    1I -   revision number of the table definition",
         "DATE-OBS   1A -   UTC start date of observations",
         "ARRNAME    0A -   name of corresponding array",
         "INSNAME    1A -   name of corresponding detector",
         "------------------------------------------------------------",
         "TARGET_ID  1I -   target number as index into OI_TARGET table",
         "TIME       1D s   UTC time of observation",
         "MJD        1D day modified Julian Day",
         "INT_TIME   1D s   integration time",
         "VIS2DATA  -1D -   squared visibility",
         "VIS2ERR   -1D -   error in squared visibility",
         "UCOORD     1D m   U coordinate of the data",
         "VCOORD     1D m   V coordinate of the data",
         "STA_INDEX  2I -   station numbers contributing to the data",
         "FLAG      -1L -   flag"])

# OI_T3 definition (1st revision):
add_def("OI_T3", 1,
        ["OI_REVN    1I -   revision number of the table definition",
         "DATE-OBS   1A -   UTC start date of observations",
         "ARRNAME    0A -   name of corresponding array",
         "INSNAME    1A -   name of corresponding detector",
         "------------------------------------------------------------",
         "TARGET_ID  1I -   target number as index into OI_TARGET table",
         "TIME       1D s   UTC time of observation",
         "MJD        1D day modified Julian Day",
         "INT_TIME   1D s   integration time",
         "T3AMP     -1D -   triple product amplitude",
         "T3AMPERR  -1D -   error in triple product amplitude",
         "T3PHI     -1D deg triple product phase",
         "T3PHIERR  -1D deg error in triple product phase",
         "U1COORD    1D m   U coordinate of baseline AB of the triangle",
         "V1COORD    1D m   V coordinate of baseline AB of the triangle",
         "U2COORD    1D m   U coordinate of baseline BC of the triangle",
         "V2COORD    1D m   V coordinate of baseline BC of the triangle",
         "STA_INDEX  3I -   station numbers contributing to the data",
         "FLAG      -1L -   flag"])

# OI_SPECTRUM definition (1st revision):
add_def("OI_SPECTRUM", 1,
        ["OI_REVN    1I -   revision number of the table definition",
         "DATE-OBS   1A -   UTC start date of observations",
         "INSNAME    1A -   name of corresponding detector",
        #"FOV        1D -   area of sky over which flux is integrated",
         "------------------------------------------------------------",
         "TARGET_ID  1I -   target number as index into OI_TARGET table",
         "MJD        1D day modified Julian Day",
         "INT_TIME   1D s   integration time",
         "FLUXDATA  -1D -   flux",
         "FLUXERR   -1D -   flux error"])


text:
# test dbr

function dummy(λ::AbstractArray{T,N}, x::S, y::S) where {T<:Number, N, S<:Number}
    Nc = ones.(length(λ)) .* (x + im*y)
end # EOF dummy(...)

# Define beam
λi = 400 # intial wavelength [nm]
λf = 1000 # final wavelength [nm]
λ = LinRange(λi, λf, λf-λi+1) # wavelength range [nm]
# λ = λi:λf
λ0 = 700. # reference wavelength
θ = [0.] # angle of incidence [degrees]
beam = PlaneWave(λ, λ0, θ)

# Define layers
l0 = Geometrical(dummy(beam.λ, 1., 0.), 0.)
l1 = Geometrical(dummy(beam.λ, 1.5, 0.1), 70.)
l2 = Optical(dummy(beam.λ, 1.1, 0.1), 1/4)
l3 = Geometrical(dummy(beam.λ, 2.5, 0.6), 0.)

# Refractive index of materials: elements represent the index position below in the structure (surface down to substrate)
nseq = [l0 l1 l2 l3]

# calculation of the electromagnetic field profile: yes (1) or no (0)
emfflag = 1
# subdivision of each layer for the calculation of the EMF
h = 10
# calculation of the photonic dispersion: yes (1) or no (0)
pbgflag = 1

# call main script
results1 = thinfilmoptics(beam, nseq, emfflag, h, pbgflag)

ext: jl


text:
"""
    Haversine(radius)

The haversine distance between two locations on a sphere of given `radius`.

Locations are described with longitude and latitude in degrees.
The computed distance has the same units as that of the radius.
"""
struct Haversine{T<:Real} <: Metric
    radius::T
end

const VecOrLengthTwoTuple{T} = Union{AbstractVector{T}, NTuple{2, T}}

function (dist::Haversine)(x::VecOrLengthTwoTuple, y::VecOrLengthTwoTuple)
    length(x) == length(y) == 2 || haversine_error()

    @inbounds begin
        # longitudes
        Δλ = deg2rad(y[1] - x[1])

        # latitudes
        φ₁ = deg2rad(x[2])
        φ₂ = deg2rad(y[2])
    end

    Δφ = φ₂ - φ₁

    # haversine formula
    a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2

    # distance on the sphere
    2 * dist.radius * asin( min(√a, one(a)) ) # take care of floating point errors
end

haversine(x::VecOrLengthTwoTuple, y::VecOrLengthTwoTuple, radius::Real) = Haversine(radius)(x, y)

@noinline haversine_error() = throw(ArgumentError("expected both inputs to have length 2 in Haversine distance"))


text:
function ALparse(userinput, logicset; verbose = true)
    if (nfeasible(logicset)==0) && !occursin(r"∈", userinput)
      replthrow("You are working with an empty set! Try inputing: a,b,c in 1:3")
      return
    end

    try
      logicset = logicalparse(userinput, logicset = logicset, verbose = verbose)
    catch ex
      replthrow("\ncommand Failed! $ex")
      return
    end

    push!(activehistory, logicset)
    update!(sessionhistory, logicset)
    setreplset!(logicset)

end


text:
export wfst2dot, create_pdf, read_wfst

# Returns a string representation of the supplied WFST in the DOT language for
# use with Graphviz.
function wfst2dot(wfst::Wfst)
    s = "digraph FST {\n"
    for node in wfst.states
        nodetext = replace(string(node), r"\"", "\\\"")
        if node in wfst.final_states
            s = string("$s\t\"$nodetext\" [shape=doublecircle, color=purple, ",
                "label=\"$nodetext$(haskey(wfst.final_weights, node) ?
                string("/", wfst.final_weights[node]) : "")\"];\n")
        elseif node in wfst.initial_states
            s = string("$s\t\"$nodetext\" [shape=circle, color=green, ",
                "label=\"$nodetext$(haskey(wfst.initial_weights, node) ?
                string("/", wfst.initial_weights[node]) : "")\"];\n")
        else
            s = "$s\t\"$nodetext\" [shape=circle];\n"
        end
    end
    for arc in wfst.arcs
        # fromtext and totext used to allow us to prepend quotes in the
        # actual nodes names with backslashes so that bash pipes the desired
        # text to graphviz.
        fromtext = replace(string(arc.from), r"\"", "\\\"")
        totext = replace(string(arc.to), r"\"", "\\\"")
        s = string("$s\t\"$fromtext\" -> \"$totext\" ",
             "[label=\"$(string(arc.input)):$(string(arc.output))",
             "/$(arc.weight)\"];\n")
    end
    s = "$s}"
    return s
end

function create_pdf(wfst::Wfst, filename::String)
    dotstring = wfst2dot(wfst)
    run(`echo $dotstring` |> `dot -Tpdf -o $filename`)
end

# Reads a text representation of a WFST in AT&T format and returns the WFST.
function read_wfst(text, semiring)
    wfst = Wfst(semiring)
    lines = split(text, "\n")
    for line in lines
        line_items = split(line)
        if length(line_items) == 5
            # Then it's an arc
            from = int(line_items[1])
            to = int(line_items[2])
            input = string(line_items[3])
            output = string(line_items[4])
            weight = float64(line_items[5])
            add_arc!(wfst, from, to, input, output, weight)
        elseif length(line_items) == 2
            # Then it specifies a final weight.
            final = int(line_items[1])
            weight = float64(line_items[2])
            add_final_state!(wfst, final, weight)
        elseif length(line_items) != 0
            error("Invalid line length")
        end
    end
    # This is going to break when we introduce other semirings because of the
    # weight.
    add_initial_state!(wfst, int(split(lines[1])[1]), 1.0)
    return wfst
end


text:
# ippf  variable to hold the the IPP function name
# sr    single rate
# mr    multirate
# Tt    datatype of filter taps (coefficients)
# Tx    datatype of signal

const WIN_BARTLETT = 0
const WIN_BLACKMAN = 1
const WIN_HAMMING  = 2
const WIN_HANN     = 3
const WIN_RECT     = 4
const LOWPASS      = 0
const HIGHPASS     = 1
const BANDPASS     = 2
const BANDSTOP     = 3

export  conv,
        xcorr,
        autocorr,
        autocorrb,
        autocorru,
        FIRFilter,
        filt,
        FIRInit

################################################################################
#               ____ ____ _  _ _  _ ____ _    _  _ ___ _ ____ _  _             #
#               |    |  | |\ | |  | |  | |    |  |  |  | |  | |\ |             #
#               |___ |__| | \|  \/  |__| |___ |__|  |  | |__| | \|             #
################################################################################

conv_types = [ :IPP32f,
               :IPP64f,
               :IPP16s ]

for ( julia_fun, ippf_prefix ) in [ ( :conv, "ippsConv" ) ]

    julia_fun! = symbol(string(julia_fun, '!'))

    for T in conv_types

        ippf = string( ippf_prefix, IPPSuffix( T ) )

        if eval(T) == IPP16s # integer versions require a scaling factor
            ippf *= "_Sfs"
            @eval begin
                function $(julia_fun!)( y::Array{$T}, x1::Array{$T}, x2::Array{$T}; scale = 0  )
                    ny  = length( y )
                    nx1 = length( x1 )
                    nx2 = length( x2 )
                    ny == nx1+nx2-1 || error("Length(y) must equal length(x1)+length(x2)-1" )
                    if ny > 3
                        @ippscall( $ippf, ( Ptr{$T},    IPPInt,     Ptr{$T},    IPPInt,     Ptr{$T},    IPPInt  ),
                                            x1,         nx1,        x2,         nx2,        y,          scale   )
                    end
                    return y
                end

                $(julia_fun)(  x1::Array{$T}, x2::Array{$T}; args... ) = $(julia_fun!)( similar(x1, length(x1) + length(x2) - 1), x1, x2; args... )

            end
        else
            @eval begin
                function $(julia_fun!)( y::Array{$T}, x1::Array{$T}, x2::Array{$T}  )
                    ny  = length( y )
                    nx1 = length( x1 )
                    nx2 = length( x2 )
                    ny == nx1+nx2-1 || error( "length(y) must be length(x1)+length(x2)-1" )
                    if ny > 3
                        @ippscall( $ippf, ( Ptr{$T},    IPPInt,     Ptr{$T},    IPPInt,     Ptr{$T} ),
                                            x1,         nx1,        x2,         nx2,        y       )
                    end
                    return y
                end

                $(julia_fun)(  x1::Array{$T}, x2::Array{$T} ) = $(julia_fun!)( similar(x1, length(x1) + length(x2) - 1), x1, x2 )
            end
        end
    end
end

################################################################################
#               ____ ____ ____ ____ ____ _    ____ ___ _ ____ _  _             #
#               |    |  | |__/ |__/ |___ |    |__|  |  | |  | |\ |             #
#               |___ |__| |  \ |  \ |___ |___ |  |  |  | |__| | \|             #
################################################################################

xcorr_types = [ :IPP32f,
                :IPP64f,
                :IPP32fc,
                :IPP64fc,
                :IPP16s  ]

for ( julia_fun, ippf_prefix ) in [ ( :xcorr, "ippsCrossCorr" ) ]

    julia_fun! = symbol(string(julia_fun, '!'))

    for T in xcorr_types

        ippf = string( ippf_prefix, IPPSuffix( T ) )

        if eval(T) == IPP16s # scaled version
            ippf *= "_Sfs"
            @eval begin
                function $(julia_fun!)( y::Array{$T}, x1::Array{$T}, x2::Array{$T}; scale = 0, lowlag = 0)
                    ny  = length( y )
                    nx1 = length( x1 )
                    nx2 = length( x2 )
                    ny > 1 && nx1 > 1 && nx2 > 1 || error()
                    if ny > 3
                        @ippscall( $ippf, ( Ptr{$T},    IPPInt,     Ptr{$T},    IPPInt,     Ptr{$T},    IPPInt, IPPInt,  IPPInt ),
                                            x1,         nx1,        x2,         nx2,        y,          ny,     lowlag,  scale  )
                    end
                    return y
                end

                $(julia_fun)(  x1::Array{$T}, x2::Array{$T}; args... ) = $(julia_fun!)( similar(x1, length(x1) + length(x2) - 1), x1, x2; args... )

            end
        else
            @eval begin
                function $(julia_fun!)( y::Array{$T}, x1::Array{$T}, x2::Array{$T}; lowlag::Integer = 0)
                    ny  = length( y )
                    nx1 = length( x1 )
                    nx2 = length( x2 )
                    ny > 1 && nx1 > 1 && nx2 > 1 || error() # TODO: Check constrainTx
                    if ny > 3
                        @ippscall( $ippf, ( Ptr{$T},    IPPInt,     Ptr{$T},    IPPInt,     Ptr{$T},    IPPInt, IPPInt  ),
                                            x1,         nx1,        x2,         nx2,        y,          ny,     lowlag  )
                    end
                    return y
                end

                $(julia_fun)(  x1::Array{$T}, x2::Array{$T}; args... ) = $(julia_fun!)( similar(x1, length(x1) + length(x2) - 1), x1, x2; args... )
            end
        end
    end
end

autocorr_types = [ :IPP16s,
                   :IPP32f,
                   :IPP64f,
                   :IPP16fc,
                   :IPP32fc,
                   :IPP64fc  ]

for ( julia_fun, ippf_prefix )  in  [   ( :autocorr,  "ippsAutoCorr"       ),
                                        ( :autocorrb, "ippsAutoCorr_NormA" ),
                                        ( :autocorru, "ippsAutoCorr_NormB" )   ]

    julia_fun! = symbol(string(julia_fun, '!'))

    for T = autocorr_types

        ippf = string( ippf_prefix, '_', IPPSuffix( T ) )

        if eval(T) == IPP16s
            ippf *= "_Sfs"
            @eval begin
                function $(julia_fun!)( y::Array{$T}, x::Array{$T}; scale = 0 )
                    ny = length( y )
                    nx = length( x )
                    ny == nx || throw( DimensionMismatch("length(y) != length(x) ") )
                    if ny > 0
                        @ippscall( $ippf, ( Ptr{$T},    IPPInt, Ptr{$T},    IPPInt, IPPInt  ),
                                            x,          nx,     y,          ny,     scale   )
                    end

                    return y
                end

                $(julia_fun)(  x::Array{$T}; args... ) = $(julia_fun!)( similar(x, length(x)*2-1), x; args... )

            end
        else
            @eval begin
                function $(julia_fun!)( y::Array{$T}, x::Array{$T}  )
                    ny = length( y )
                    nx = length( x )
                    ny == nx|| throw( DimensionMismatch("length(y) != length(x) ") ) #todo: check this
                    if ny > 0
                        @ippscall( $ippf, ( Ptr{$T},    IPPInt, Ptr{$T},    IPPInt  ),
                                            x,          nx,     y,          ny      )
                    end

                    return y
                end

                $(julia_fun)(  x::Array{$T} ) = $(julia_fun!)( similar(x), x )

            end
        end
    end
end

################################################################################
#      ____ _ ____    ____ _ _    ___ ____ ____    ___ _   _ ___  ____         #
#      |___ | |__/    |___ | |     |  |___ |__/     |   \_/  |__] |___         #
#      |    | |  \    |    | |___  |  |___ |  \     |    |   |    |___         #
################################################################################
# These are common to most of the ipps FIR Filter functions
# I haven't included functions that require integer scaling.
# They take difference parameters. Hope I haven't painted
# myself into a corner.
FIRFilterTypes =    [   ( :IPP32f,      :IPP32f  ),
                        ( :IPP32fc,     :IPP32fc ),
                        ( :IPP64f,      :IPP64f  ),
                        ( :IPP64fc,     :IPP64fc ),
                        ( :IPP32f,      :IPP16s  ),
                        ( :IPP64f,      :IPP16s  ),
                        ( :IPP64f,      :IPP32f  ),
                        ( :IPP64f,      :Int32   ),
                        ( :IPP64fc,     :IPP16sc ),
                        ( :IPP64fc,     :IPP32sc ),
                        ( :IPP64fc,     :IPP32fc )   ]

# IPP assings a pointer to state
# This points to an invisuble-to-us struct in buffer
# Intel's documentation is weak in this area, but we allocate the buffer, they
# instantsiat a stuct inside that buffer, and use state to point the location of
# that state

# Both the singe-rate and multirate states are the same, but are in different
# types to help with dispatching

# Single rate state
type FIRSRState
    pointer::Ptr{Void}
    buffer::Vector{Uint8}
end

# Multirate state
type FIRMRState
    pointer::Ptr{Void}
    buffer::Vector{Uint8}
end

# Polyphase state
type FIRPPState
    pointer::Ptr{Void}
    buffer::Vector{Uint8}
end

# FIRFilter type, its type signature includes the FIRState type, which is how functions taking
# FIRFilter as a parapter get diuspatched to the correct single-rate or multirate IPP functions
type FIRFilter{Tt, Tx, Ts}
    taps::Vector{Tt}
    delayLine::Vector{Tx}
    state::Ts
    upFactor::Int
    downFactor::Int
    upPhase::Int
    downPhase::Int
end

# Single rate constructor
function FIRFilter{Tx, Tt}( ::Type{Tx}, taps::Vector{Tt} )
    tapsLen   = length( taps )
    pointer   = convert( Ptr{Uint8}, 0 )
    bufLen    = FIRRequiredStateSize( Tt, Tx, tapsLen )
    delayLine = zeros( Tx, tapsLen )
    buffer    = zeros( Uint8, bufLen )
    state     = FIRSRState( pointer, buffer )
    FIRInit!( state, taps, delayLine )
    FIRFilter( taps, delayLine, state, 1, 1, 0, 0 )
end

# Multirate constructor
function FIRFilter{ Tx, Tt }( ::Type{Tx}, taps::Array{Tt, 1}, upFactor, downFactor, upPhase = 0, downPhase = 0 )
    tapsLen   = length( taps )
    pointer   = convert( Ptr{Uint8}, 0 )
    bufLen    = FIRRequiredStateSize( Tt, Tx, tapsLen, upFactor, downFactor )
    delayLine = zeros( Tx, tapsLen )
    buffer    = zeros( Uint8, bufLen )
    state     = FIRMRState( pointer, buffer )
    FIRInit!( state, taps, delayLine, upFactor, downFactor, upPhase, downPhase )
    FIRFilter( taps, delayLine, state, upFactor, downFactor, upPhase, downPhase )
end

FIRFilter{ Tx, Tt }( ::Type{Tx}, taps::Array{Tt, 1}, resampleRatio::Rational, upPhase = 0, downPhase = 0 ) = FIRFilter( Tx, taps, num(resampleRatio), den(resampleRatio), upPhase, downPhase )

################################################################################
#           ____ _ ____    ____ _ _    ___ ____ ____    _ _  _ _ ___           #
#           |___ | |__/    |___ | |     |  |___ |__/    | |\ | |  |            #
#           |    | |  \    |    | |___  |  |___ |  \    | | \| |  |            #
################################################################################
#   Name:         ippsFIRGetStateSize, ippsFIRMRGetStateSize,
#                 ippsFIRInit, ippsFIRMRInit
#   Purpose:      ippsFIRGetStateSize   - calculates the size of the FIR Statestructure;
#                 ippsFIRInit           - initialize FIR state - set taps and delay line
#                 using external memory buffer;
#   Parameters:
#       pTaps       - pointer to the filter coefficients;
#       tapsLen     - number of coefficients;
#       pDlyLine    - pointer to the delay line values, can be NULL;
#       ppState     - pointer to the FIR state created or NULL;
#       upFactor    - multi-rate up factor;
#       upPhase     - multi-rate up phase;
#       downFactor  - multi-rate down factor;
#       downPhase   - multi-rate down phase;
#       pStateSize  - pointer where to store the calculated FIR State structure
#                                                              size (in bytes);
#    Return:
#       status      - status value returned, its value are
#          ippStsNullPtrErr       - pointer(s) to the data is NULL
#          ippStsFIRLenErr        - tapsLen <= 0
#          ippStsFIRMRFactorErr   - factor <= 0
#          ippStsFIRMRPhaseErr    - phase < 0 || factor <= phase
#          ippStsNoErr            - otherwise

# IppStatus ippsFIRMRRequiredStateSize_64fc(int tapsLen, int upFactor, int downFactor, int* pBufferSize);

for ( julia_fun, ippf_prefix1, ippf_prefix2 )  in  [ (   :FIRRequiredStateSize,  "ippsFIR",  "GetStateSize"  ) ]
    for ( Tt, Tx ) in FIRFilterTypes

        ippf_suffix = Tt == Tx ? IPPSuffix( Tt ) : IPPSuffix( (Tt, Tx) )

        # Single rate version
        ippfsr = string( ippf_prefix1, ippf_prefix2, ippf_suffix )
        @eval begin
            function $(julia_fun)( ::Type{$Tt}, ::Type{$Tx}, tapslen::Integer  )
                buffersize = IPPInt[0]  # the function will fill this value
                tapslen > 0 || error()
                @ippscall( $ippfsr, ( IPPInt,     Ptr{IPPInt} ),
                                      tapslen,    buffersize    )

                return buffersize[1]
            end
        end

        # Multirate version
        ippfmr = string( ippf_prefix1, "MR", ippf_prefix2, ippf_suffix )
        @eval begin
            function $(julia_fun)( ::Type{$Tt}, ::Type{$Tx}, tapslen::Integer, upFactor::Integer, downFactor::Integer )
                buffersize = IPPInt[0]  # the function will fill this value
                tapslen > 0 || error() # todo: check this
                @ippscall( $ippfmr, ( IPPInt,   IPPInt,    IPPInt,      Ptr{IPPInt} ),
                                      tapslen,  upFactor,  downFactor,  buffersize    )

                return buffersize[1]
            end
            $(julia_fun)( ::Type{$Tt}, ::Type{$Tx}, tapslen::Integer, resamp_ratio::Rational ) = $(julia_fun)( $Tt, $Tx, tapslen, num(resamp_ratio), den(resamp_ratio) )
        end
    end
end

for ( julia_fun, ippf_prefix1, ippf_prefix2 )  in  [ (   :FIRInit,  "ippsFIR",  "Init"  ) ]
    for ( Tt, Tx ) in FIRFilterTypes

        ippf_suffix = Tt == Tx ? IPPSuffix( Tt ) : IPPSuffix( (Tt, Tx) )
        julia_fun! = symbol(string(julia_fun, '!'))

        # Single rate version
        ippfsr     = string( ippf_prefix1, ippf_prefix2, ippf_suffix )
        @eval begin

            function $(julia_fun!)( state::FIRSRState, taps::Vector{ $Tt }, delayLine::Vector{ $Tx } )
                ntaps          = length( taps )
                ndelayline     = length( delayLine )
                statePtr       = Array( Ptr{Void}, 1 )
                @ippscall( $ippfsr, (   Ptr{Void},         Ptr{$Tt},    IPPInt,     Ptr{$Tx},      Ptr{Uint8}   ),
                                        pointer(statePtr), taps,        ntaps,      delayLine,     state.buffer )
                state.pointer = statePtr[1]
                nothing
            end

        end

        # Multirate version
        ippfmr  = string( ippf_prefix1, "MR", ippf_prefix2, ippf_suffix )
        @eval begin

            function $(julia_fun!)(  state::FIRMRState, taps::Vector{ $Tt }, delayLine::Vector{ $Tx }, upFactor, downFactor, upPhase, downPhase )
                ntaps          = length( taps )
                ndelayline     = length( delayLine )
                statePtr       = Array( Ptr{Void}, 1 )
                0 <= upPhase < upFactor && 0 <= downPhase < downFactor || error()
                @ippscall( $ippfmr, (   Ptr{Void},         Ptr{$Tt}, IPPInt, IPPInt,    IPPInt,  IPPInt,     IPPInt,    Ptr{$Tx},  Ptr{Uint8}   ),
                                        pointer(statePtr), taps,     ntaps,  upFactor,  upPhase, downFactor, downPhase, delayLine, state.buffer )
                state.pointer = statePtr[1]
                nothing
            end

        end
    end
end

################################################################################
#      ____ _ ____    ____ _ _    ___ ____ ____    ____ _  _ ____ ____         #
#      |___ | |__/    |___ | |     |  |___ |__/    |___  \/  |___ |            #
#      |    | |  \    |    | |___  |  |___ |  \    |___ _/\_ |___ |___         #
################################################################################
#
#  Names:         ippsFIR
#  Purpose:       FIR filter. Vector filtering
#  Parameters:
#      pSrcDst     - pointer to the input/output vector in in-place operation
#      pSrc        - pointer to the input vector
#      pDst        - pointer to the output vector
#      numIters    - number iterations (for single-rate equal length data vector)
#      pState      - pointer to the filter state
#      scaleFactor - scale factor value
#  Return:
#      ippSTxContextMatchErr  - wrong state identifier
#      ippSTxNullPtrErr       - pointer(s) to the data is NULL
#      ippSTxSizeErr          - numIters is less or equal zero
#      ippSTxNoErr            - otherwise
#  Note: for Multi-Rate filtering
#          length pSrc = numIters*downFactor
#          length pDst = numIters*upFactor
#               for inplace functions max this values
#
for ( julia_fun, ippf_prefix )  in  [ (   :filt,  "ippsFIR"  ) ]
    for ( Tt, Tx ) in FIRFilterTypes

        ippf_suffix = Tt == Tx ? IPPSuffix( Tt ) : IPPSuffix( (Tt, Tx) )
        julia_fun!  = symbol(string(julia_fun, '!'))

        # single rate
        ippfsr = string( ippf_prefix, ippf_suffix )
        @eval begin

            # TODO: make a lower level version of filt so direct form filters dont need to create a fir filter object
            #       they should only need to call getsize and init
            # using a pre-instanciated FIRFilter object, and a preallocated buffer
            function $(julia_fun!)( self::FIRFilter{$Tt, $Tx, FIRSRState}, buffer::Vector{ $Tx }, signal::Vector{ $Tx } )
                sigLen = length( signal )
                bufLen = length( buffer )
                bufLen >= sigLen || error( "output buffer is too small" )
                @ippscall( $ippfsr,  (  Ptr{$Tx},       Ptr{$Tx},       IPPInt,    Ptr{Void}                    ),
                                        signal,         buffer,         sigLen,    self.state.pointer           )
                buffer
            end

            # using a pre-instanciated FIRFilter object
            $(julia_fun)( self::FIRFilter{$Tt, $Tx, FIRSRState}, signal::Vector{$Tx} ) = $(julia_fun!)( self, similar( signal ), signal )

            # direct form, with taps and signal, and a preallocated buffer
            function $(julia_fun!)( buffer::Vector{ $Tx }, taps::Vector{ $Tt }, signal::Vector{ $Tx })
                self = FIRFilter( $Tx, taps )
                filt!( self, buffer, signal )
            end

            # direct form, with taps and signal
            $(julia_fun)( taps::Vector{ $Tt }, signal::Vector{ $Tx } ) = $(julia_fun!)( similar(signal), taps, signal )

        end

        # multirate
        ippfmr = string( ippf_prefix, "MR", ippf_suffix )
        @eval begin

            # using a pre-instanciated FIRFilter object, and a preallocated buffer
            function $(julia_fun!)( self::FIRFilter{ $Tt, $Tx, FIRMRState }, buffer::Vector{ $Tx }, signal::Array{$Tx} )
                sigLen = length( signal )
                outLen = length( buffer )
                sigLen % self.downFactor == 0 || error( "signal length must be an integer multiple of decimation")
                sigLen * self.upFactor == outLen * self.downFactor || error( "buffer length needs to be >= sigal len * interpolation/decimation" )
                iterations  = int( sigLen/self.downFactor )
                @ippscall( $ippfsr,  (  Ptr{$Tx},  Ptr{$Tx},  IPPInt,     Ptr{Void}            ),
                                        signal,    buffer,    iterations, self.state.pointer   )
                buffer
            end

            # using a pre-instanciated FIRFilter object
            function $(julia_fun)( self::FIRFilter{$Tt, $Tx, FIRMRState}, signal::Vector{ $Tx } )
                sigLen = length( signal )
                bufLen = int( sigLen * self.upFactor / self.downFactor )
                buffer = similar( signal, bufLen )
                $(julia_fun!)( self, buffer, signal )
            end

            # direct form, with taps and signal, and a preallocated buffer
            function $(julia_fun!)( buffer::Vector{ $Tx }, taps::Vector{ $Tt }, signal::Vector{ $Tx }, upFactor::Integer, downFactor::Integer, upPhase = 0, downPhase = 0)
                self = FIRFilter( $Tx, taps, upFactor, downFactor, upPhase, downPhase )
                $(julia_fun!)( self, buffer, signal )
            end

            # direct form, with taps and signal
            function $(julia_fun)( taps::Vector{ $Tt }, signal::Vector{ $Tx }, upFactor::Integer, downFactor::Integer, upPhase = 0, downPhase = 0)
                sigLen = length( signal )
                bufLen = int( sigLen * upFactor / downFactor )
                buffer = similar( signal, bufLen )
                $(julia_fun!)( buffer, taps, signal, upFactor, downFactor, upPhase, downPhase )
            end
        end
    end
end

######################################################################################
# ____ _ ____    ___ ____ ___  ____    ____ ____ _  _ ____ ____ ____ ___ _ ____ _  _ #
# |___ | |__/     |  |__| |__] [__     | __ |___ |\ | |___ |__/ |__|  |  | |  | |\ | #
# |    | |  \     |  |  | |    ___]    |__] |___ | \| |___ |  \ |  |  |  | |__| | \| #
######################################################################################
#  Names:      ippsFIRGenLowpass_64f, ippsFIRGenHighpass_64f, ippsFIRGenBandpass_64f
#              ippsFIRGenBandstop_64f
#
#  Purpose:    This function computes the lowpass FIR filter coefficients
#              by windowing of ideal (infinite) filter coefficients segment
#
#  Parameters:
#      rfreq             cut off frequency (0 < rfreq < 0.5)
#
#      taps              pointer to the array which specifies
#                        the filter coefficients;
#
#      tapsLen           the number of taps in taps[] array (tapsLen>=5);
#
#      winType           the ippwindow switch variable,
#                        which specifies the smoothing window type;
#
#      doNormal          if doNormal=0 the functions calculates
#                        non-normalized sequence of filter coefficients,
#                        in other cases the sequence of coefficients
#                        will be normalized.
#  Return:
#   ippStsNullPtrErr     the null pointer to taps[] array pass to function
#   ippStsSizeErr        the length of coefficient's array is less than five
#   ippStsSizeErr        the low or high frequency isn't satisfy
#                                    the condition 0 < rLowFreq < 0.5
#   ippStsNoErr          otherwise
#
#
#  IPPAPI(IppStatus, ippsFIRGenLowpass_64f( Ipp64f rfreq,    Ipp64f* taps,      int tapsLen, IppWinType winType, IppBool doNormal))
#  IPPAPI(IppStatus, ippsFIRGenHighpass_64f(Ipp64f rfreq,    Ipp64f* taps,      int tapsLen, IppWinType winType, IppBool doNormal))
#  IPPAPI(IppStatus, ippsFIRGenBandpass_64f(Ipp64f rLowFreq, Ipp64f  rHighFreq, Ipp64f* taps, int tapsLen, IppWinType winType, IppBool doNormal))
#  IPPAPI(IppStatus, ippsFIRGenBandstop_64f(Ipp64f rLowFreq, Ipp64f  rHighFreq, Ipp64f* taps, int tapsLen, IppWinType winType, IppBool doNormal))

for ( julia_fun, ippf_prefix )  in  [   ( :highpass,  "ippsFIRGenHighpass" ),
                                        ( :lowpass,   "ippsFIRGenLowpass"  )  ]
    julia_fun! = symbol(string(julia_fun, '!'))

    for T in [ IPP64f ]

        ippf_suffix = IPPSuffix( T )
        ippf        = string( ippf_prefix, ippf_suffix )

        @eval begin
            function $(julia_fun!){Tt}( buffer::Vector{Tt}, ω::FloatingPoint, window = WIN_HAMMING, normalize = true )
                Ntaps = length ( buffer )
                0.0 < ω < 0.5   || error( "cutoff frequency ω must satisfy:  0.0 < ω < 0.5" )
                Ntaps >= 5      || error( "number of taps must satisfy: Ntaps >= 5")
                innerBuffer = Tt == $T ? buffer : Array( $T, Ntaps )
                @ippscall( $ippf,  (  IPP64f, Ptr{$T},      IPPInt, IPPInt,     Bool        ),
                                      ω,      innerBuffer,  Ntaps,  window, normalize   )
                return Tt == $T ? buffer : [ buffer[i] = convert( Tt, innerBuffer[i] ) for i = 1:Ntaps ]
            end # function

            $(julia_fun){Tt}( ::Type{Tt}, Ntaps::Integer, ω::FloatingPoint, window = WIN_HAMMING, normalize = true ) = $(julia_fun!)( Array(Tt, Ntaps), ω, window, normalize )
        end # eval
    end # type loop
end # function loop

for ( julia_fun, ippf_prefix )  in  [   ( :bandpass,  "ippsFIRGenBandpass"  ),
                                        ( :bandstop,  "ippsFIRGenBandstop"  )  ]

    julia_fun! = symbol(string(julia_fun, '!'))

    for T in [ IPP64f ]

        ippf_suffix = IPPSuffix( T )
        ippf        = string( ippf_prefix, ippf_suffix )

        @eval begin
            function $(julia_fun!){Tt}( buffer::Vector{Tt}, ω1::FloatingPoint, ω2::FloatingPoint, window = WIN_HAMMING, normalize = true )
                Ntaps = length ( buffer )
                0.0 < ω1 < ω2 < 0.5 || error( "cutoff frequencies ω1, ω2  must satisfy:  0.0 < ω1 < ω2 < 0.5" )
                Ntaps >= 5          || error( "number of taps must satisfy: Ntaps >= 5" )
                innerBuffer = Tt == $T ? buffer : Array( $T, Ntaps )
                @ippscall( $ippf,  (  IPP64f, IPP64f, Ptr{$T},      IPPInt, IPPInt,     Bool        ),
                                      ω1,     ω2,     innerBuffer,  Ntaps,  window, normalize   )
                return Tt == $T ? buffer : [ buffer[i] = innerBuffer[i] for i = 1:Ntaps ]
            end # function

            $(julia_fun){Tt}( ::Type{Tt}, Ntaps::Integer, ω1::FloatingPoint, ω2::FloatingPoint, window = WIN_HAMMING, normalize = true ) = $(julia_fun!)( Array(Tt, Ntaps), ω1, ω2, window, normalize )
        end # eval
    end # type loop
end # function loop


text:
using ModelingToolkit, Test, SparseArrays
@variables a b c

# Auxiliary Functions and Constants
get_sparsity_pattern(h::Array{Expression}) = sparse(Int64.(map(~, h .=== ModelingToolkit.Constant(0))))
get_sparsity_pattern(h::SparseMatrixCSC{Expression,Int64}) = sparse(Int64.(map(~, h .=== ModelingToolkit.Constant(0))))
get_sparsity_pattern(h::SparseVector{Expression,Int64}) = sparse(Int64.(map(~, h .=== ModelingToolkit.Constant(0))))

input = [1, 2, 3]

# ===== Dense tests =====
# Arrays of Matrices
h_dense_arraymat = [[a 1; b 0], [0 0; 0 0], [a c; 1 0]] # empty array support required
function h_dense_arraymat_julia!(out, x)
    a, b, c = x
    out[1] .= [a[1] 1; b[1] 0]
    out[2] .= [0 0; 0 0]
    out[3] .= [a[1] c[1]; 1 0]
end

h_dense_arraymat_str = ModelingToolkit.build_function(h_dense_arraymat, [a, b, c])
h_dense_arraymat_ip! = eval(h_dense_arraymat_str[2])
out_1_arraymat = [Array{Int64}(undef, 2, 2) for i in 1:3]
out_2_arraymat = [similar(x) for x in out_1_arraymat]
h_dense_arraymat_julia!(out_1_arraymat, input)
h_dense_arraymat_ip!(out_2_arraymat, input)
@test out_1_arraymat == out_2_arraymat

# Arrays of Matrices, Heterogeneous element types
test_exp = [exp(a) * exp(b), a]
h_dense_arraymat_het = [ModelingToolkit.hessian(t, [a, b]) for t in test_exp]
function h_dense_arraymat_het_julia!(out, x)
    a, b, c = x
    out[1] .= [exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1]); exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1])]
    out[2] .= [0 0; 0 0]
end

h_dense_arraymat_het_str = ModelingToolkit.build_function(h_dense_arraymat_het, [a, b, c])
h_dense_arraymat_het_ip! = eval(h_dense_arraymat_het_str[2])
out_1_arraymat_het = [Array{Float64}(undef, 2, 2) for i in 1:2]
out_2_arraymat_het = [similar(x) for x in out_1_arraymat_het]
h_dense_arraymat_het_julia!(out_1_arraymat_het, input)
h_dense_arraymat_het_ip!(out_2_arraymat_het, input)
@test out_1_arraymat_het == out_2_arraymat_het

# Arrays of 1D Vectors
h_dense_arrayvec = [[a, 0, c], [0, 0, 0], [1, a, b]] # same for empty vectors, etc.
function h_dense_arrayvec_julia!(out, x)
    a, b, c = x
    out[1] .= [a[1], 0, c[1]]
    out[2] .= [0, 0, 0]
    out[3] .= [1, a[1], b[1]]
end

h_dense_arrayvec_str = ModelingToolkit.build_function(h_dense_arrayvec, [a, b, c])
h_dense_arrayvec_ip! = eval(h_dense_arrayvec_str[2])
out_1_arrayvec = [Vector{Int64}(undef, 3) for i in 1:3]
out_2_arrayvec = [Vector{Int64}(undef, 3) for i in 1:3]
h_dense_arrayvec_julia!(out_1_arrayvec, input)
h_dense_arrayvec_ip!(out_2_arrayvec, input)
@test out_1_arrayvec == out_2_arrayvec

# Arrays of Arrays of Matrices
h_dense_arrayNestedMat = [[[a 1; b 0], [0 0; 0 0]], [[b 1; a 0], [b c; 0 1]]]
function h_dense_arrayNestedMat_julia!(out, x)
    a, b, c = x
    out[1][1] .= [a[1] 1; b[1] 0]
    out[1][2] .= [0 0; 0 0]
    out[2][1] .= [b[1] 1; a[1] 0]
    out[2][2] .= [b[1] c[1]; 0 1]
end

h_dense_arrayNestedMat_str = ModelingToolkit.build_function(h_dense_arrayNestedMat, [a, b, c])
h_dense_arrayNestedMat_ip! = eval(h_dense_arrayNestedMat_str[2])
out_1_arrayNestedMat = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]] # avoid undef broadcasting issue
out_2_arrayNestedMat = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]]
h_dense_arrayNestedMat_julia!(out_1_arrayNestedMat, input)
h_dense_arrayNestedMat_ip!(out_2_arrayNestedMat, input)
@test out_1_arrayNestedMat == out_2_arrayNestedMat

# Arrays of Arrays of Matrices, Heterogeneous element types
test_exp = [exp(a) * exp(b), a]
h_dense_arrayNestedMat_het = [[ModelingToolkit.hessian(t, [a, b]) for t in test_exp], [[ModelingToolkit.Constant(0) ModelingToolkit.Constant(0); ModelingToolkit.Constant(0) ModelingToolkit.Constant(0)], [ModelingToolkit.Constant(0) ModelingToolkit.Constant(0); ModelingToolkit.Constant(0) ModelingToolkit.Constant(0)]]]
function h_dense_arrayNestedMat_het_julia!(out, x)
    a, b, c = x
    out[1][1] .= [exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1]); exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1])]
    out[1][2] .= [0 0; 0 0]
    out[2][1] .= [0 0; 0 0]
    out[2][2] .= [0 0; 0 0]
end
h_dense_arrayNestedMat_het_str = ModelingToolkit.build_function(h_dense_arrayNestedMat_het, [a, b, c])
h_dense_arrayNestedMat_het_ip! = eval(h_dense_arrayNestedMat_het_str[2])
out_1_arrayNestedMat_het = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]] # avoid undef broadcasting issue
out_2_arrayNestedMat_het = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]]
h_dense_arrayNestedMat_julia!(out_1_arrayNestedMat_het, input)
h_dense_arrayNestedMat_ip!(out_2_arrayNestedMat_het, input)
@test out_1_arrayNestedMat_het == out_2_arrayNestedMat_het

# ===== Sparse tests =====
# Array of Matrices
h_sparse_arraymat = sparse.([[a 1; b 0], [0 0; 0 0], [a c; 1 0]])
function h_sparse_arraymat_julia!(out, x)
    a, b, c = x
    out[1][1, 1] = a[1]
    out[1][1, 2] = 1
    out[1][2, 1] = b[1]
    out[2] = sparse([0 0; 0 0]) # no undef constructor for SparseMatrixCSC
    out[3][1, 1] = a[1]
    out[3][1, 2] = c[1]
    out[3][2, 1] = 1
end

h_sparse_arraymat_str = ModelingToolkit.build_function(h_sparse_arraymat, [a, b, c])
h_sparse_arraymat_ip! = eval(h_sparse_arraymat_str[2])
h_sparse_arraymat_sparsity_patterns = map(get_sparsity_pattern, h_sparse_arraymat)
out_1_arraymat = [similar(h) for h in h_sparse_arraymat_sparsity_patterns]
out_2_arraymat = [similar(h) for h in h_sparse_arraymat_sparsity_patterns] # can't do similar() because it will just be #undef, with the wrong sparsity pattern
h_sparse_arraymat_julia!(out_1_arraymat, input)
h_sparse_arraymat_ip!(out_2_arraymat, input)
@test out_1_arraymat == out_2_arraymat

# Array of 1D Vectors
h_sparse_arrayvec = sparse.([[a, 0, c], [0, 0, 0], [1, a, b]])
function h_sparse_arrayvec_julia!(out, x)
    a, b, c = x
    out[1][1] = a[1]
    out[1][3] = c[1]
    out[2] = sparse([0, 0, 0]) # necessary because sparsity pattern is 3 elements with 0 stored, not 0 elements
    out[3][1] = 1
    out[3][2] = a[1]
    out[3][3] = b[1]
end

h_sparse_arrayvec_str = ModelingToolkit.build_function(h_sparse_arrayvec, [a, b, c])
h_sparse_arrayvec_ip! = eval(h_sparse_arrayvec_str[2])
h_sparse_arrayvec_sparsity_patterns = map(get_sparsity_pattern, h_sparse_arrayvec)
out_1_arrayvec = [similar(h) for h in h_sparse_arrayvec_sparsity_patterns]
out_2_arrayvec = [similar(h) for h in h_sparse_arrayvec_sparsity_patterns]
h_sparse_arrayvec_julia!(out_1_arrayvec, input)
h_sparse_arrayvec_ip!(out_2_arrayvec, input)
@test out_1_arrayvec == out_2_arrayvec

# Arrays of Arrays of Matrices
h_sparse_arrayNestedMat = [sparse.([[a 1; b 0], [0 0; 0 0]]), sparse.([[b 1; a 0], [b c; 0 1]])]
function h_sparse_arrayNestedMat_julia!(out, x)
    a, b, c = x
    out[1][1][1, 1] = a[1]
    out[1][1][1, 2] = 1
    out[1][1][2, 1] = b[1]
    out[1][2] = sparse([0 0; 0 0])
    out[2][1][1, 1] = b[1]
    out[2][1][1, 2] = 1
    out[2][1][2, 1] = a[1]
    out[2][2][1, 1] = b[1]
    out[2][2][1, 2] = c[1]
    out[2][2][2, 2] = 1
end

h_sparse_arrayNestedMat_str = ModelingToolkit.build_function(h_sparse_arrayNestedMat, [a, b, c])
h_sparse_arrayNestedMat_ip! = eval(h_sparse_arrayNestedMat_str[2])
h_sparse_arrayNestedMat_sparsity_patterns = [map(get_sparsity_pattern, h) for h in h_sparse_arrayNestedMat]
out_1_arrayNestedMat = [[similar(h_sub) for h_sub in h] for h in h_sparse_arrayNestedMat_sparsity_patterns]
out_2_arrayNestedMat = [[similar(h_sub) for h_sub in h] for h in h_sparse_arrayNestedMat_sparsity_patterns]
h_sparse_arrayNestedMat_julia!(out_1_arrayNestedMat, input)
h_sparse_arrayNestedMat_ip!(out_2_arrayNestedMat, input)
@test out_1_arrayNestedMat == out_2_arrayNestedMat

# Additional Tests
# Returning 0-element structures (corresponding to empty Jacobians)
# Arrays of Matrices
h_empty = [[a b; c 0], Array{Expression,2}(undef, 0,0)]
h_empty_str = ModelingToolkit.build_function(h_empty, [a, b, c])
h_empty_ip! = eval(h_empty_str[2])
out = [Matrix{Int64}(undef, 2, 2), Matrix{Int64}(undef, 0, 0)]
h_empty_ip!(out, input) # should just not fail

# Array of Vectors
h_empty_vec = [[a, b, c, 0], Vector{Expression}(undef,0)]
h_empty_vec_str = ModelingToolkit.build_function(h_empty_vec, [a, b, c])
h_empty_vec_ip! = eval(h_empty_vec_str[2])
out = [Vector{Int64}(undef, 4), Vector{Int64}(undef, 0)]
h_empty_vec_ip!(out, input) # should just not fail

# Arrays of Arrays of Matrices
h_emptyNested = [[[a b; c 0]], Array{Array{Expression, 2}}(undef, 0)] # emptyNested array of arrays
h_emptyNested_str = ModelingToolkit.build_function(h_emptyNested, [a, b, c])
h_emptyNested_ip! = eval(h_emptyNested_str[2])
out = [[[1 2;3 4]], Array{Array{Int64,2},1}(undef, 0)]
h_emptyNested_ip!(out, input) # should just not fail


text:
using Gnuplot
function hline_at(y; color = "black", lw = 2)
    "set arrow from graph 0, first $(y) to graph 1, first $(y) nohead lc rgb '$(color)' lw $(lw)"
end
@gp [2,3] [3,3] "w l notit lc 'red'" # horizontal line from 2 to 3 at 1
@gp :- xrange =(0,4) yrange=(0,4)
# full horizontal line at y = 1
@gp :- hline_at(1; lw = 3)
@gp :- hline_at(2; color = "blue")
save(term ="pngcairo size 600,400", output ="hlines.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using Test
using GeoParams

@testset "MaterialParameters" begin

# This tests the MaterialParameters structure
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);

# Define a struct for a first phase
Phase1 = SetMaterialParams(Name="test1", Phase=22,
                CreepLaws= (PowerlawViscous(), LinearViscous(η=1e21Pa*s)),
                Gravity  = ConstantGravity(g=11.0m/s^2),
                Density  = ConstantDensity(),
                CharDim  = CharUnits_GEO);

@test   Phase1.Density[1].ρ.val ≈ 2.9e-16 rtol=1e-6
@test   Phase1.Gravity[1].g*1.0 ≈ 1.1e19  rtol=1e-6
@test   Phase1.CreepLaws[1].η0*1.0==0.1
@test   Phase1.CreepLaws[2].η*1.0==100.0

Phase2 = SetMaterialParams(Name="test1",Phase=2,
                CreepLaws= (LinearViscous(), LinearViscous(η=1e21Pa*s)),
                Density  = ConstantDensity());

# Non-dimensionalize all values:
Phase2 = nondimensionalize(Phase2, CharUnits_GEO)
@test Phase2.CreepLaws[2].η.val ≈ 100.0
@test Phase2.Density[1].ρ.val ≈ 2.9e-16

# Dimensionalize all values again:
Phase2 = dimensionalize(Phase2, CharUnits_GEO)
@test Phase2.Density[1].ρ.val ≈ 2900
@test Phase2.CreepLaws[2].η.val ≈ 1e21

# Create array with several phases
MatParam        =   Vector{MaterialParams}(undef, 2);
Phase           =   1;
MatParam[Phase] =   SetMaterialParams(Name="Upper Crust", Phase=Phase,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density  = ConstantDensity(ρ=2900kg/m^3));
Phase           =   2;
MatParam[Phase] =   SetMaterialParams(Name="Lower Crust", Phase=Phase,
                        CreepLaws= (PowerlawViscous(n=5.0), LinearViscous(η=1e21Pa*s)),
                        Density  = PT_Density(ρ0=3000kg/m^3));

@test MatParam[2].Density[1].α.val ≈  3.0e-5
@test MatParam[2].CreepLaws[1].n.val == 5.0

# test adding phase Diagrams
MatParam        =   Array{MaterialParams, 1}(undef, 1);
Phase           =   1;
MatParam[Phase] =   SetMaterialParams(Name="Mantle", Phase=Phase,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = PerpleX_LaMEM_Diagram("test_data/Peridotite.in"), CharDim=GEO_units());

@test MatParam[1].Density[1].Rho(1,10) ≈ 3.1688441499707393e-18

end


text:
struct DivConformity <: Conformity end

struct RaviartThomas <: ReferenceFEName end

const raviart_thomas = RaviartThomas()

"""
    RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et

The `order` argument has the following meaning: the divergence of the  functions in this basis
is in the Q space of degree `order`.

"""
function RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et

  D = num_dims(p)

  if is_n_cube(p)
    prebasis = QCurlGradMonomialBasis{D}(et,order)
  elseif is_simplex(p)
    prebasis = PCurlGradMonomialBasis{D}(et,order)
  else
    @notimplemented "H(div) Reference FE only available for cubes and simplices"
  end

  nf_nodes, nf_moments = _RT_nodes_and_moments(et,p,order,GenericField(identity))

  face_own_dofs = _face_own_dofs_from_moments(nf_moments)

  face_dofs = face_own_dofs

  dof_basis = MomentBasedDofBasis(nf_nodes, nf_moments)

  ndofs = num_dofs(dof_basis)

  metadata = nothing

  reffe = GenericRefFE{RaviartThomas}(
    ndofs,
    p,
    prebasis,
    dof_basis,
    DivConformity(),
    metadata,
    face_dofs)

  reffe
end

function ReferenceFE(p::Polytope,::RaviartThomas, order)
  RaviartThomasRefFE(Float64,p,order)
end

function ReferenceFE(p::Polytope,::RaviartThomas,::Type{T}, order) where T
  RaviartThomasRefFE(T,p,order)
end

function Conformity(reffe::GenericRefFE{RaviartThomas},sym::Symbol)
  hdiv = (:Hdiv,:HDiv)
  if sym == :L2
    L2Conformity()
  elseif sym in hdiv
    DivConformity()
  else
    @unreachable """\n
    It is not possible to use conformity = $sym on a Raviart Thomas reference FE.

    Possible values of conformity for this reference fe are $((:L2, hdiv...)).
    """
  end
end

function get_face_own_dofs(reffe::GenericRefFE{RaviartThomas}, conf::DivConformity)
  get_face_dofs(reffe)
end

function _RT_nodes_and_moments(::Type{et}, p::Polytope, order::Integer, phi::Field) where et

  D = num_dims(p)
  ft = VectorValue{D,et}
  pt = Point{D,et}

  nf_nodes = [ zeros(pt,0) for face in 1:num_faces(p)]
  nf_moments = [ zeros(ft,0,0) for face in 1:num_faces(p)]

  fcips, fmoments = _RT_face_values(p,et,order,phi)
  frange = get_dimrange(p,D-1)
  nf_nodes[frange] = fcips
  nf_moments[frange] = fmoments

  if (order > 0)
    ccips, cmoments = _RT_cell_values(p,et,order,phi)
    crange = get_dimrange(p,D)
    nf_nodes[crange] = ccips
    nf_moments[crange] = cmoments
  end

  nf_nodes, nf_moments
end

# Ref FE to faces geomaps
function _ref_face_to_faces_geomap(p,fp)
  cfvs = get_face_coordinates(p,num_dims(fp))
  nc = length(cfvs)
  freffe = LagrangianRefFE(Float64,fp,1)
  fshfs = get_shapefuns(freffe)
  cfshfs = fill(fshfs, nc)
  fgeomap = lazy_map(linear_combination,cfvs,cfshfs)
end

function _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)
  nc = length(fgeomap)
  c_fips = fill(fips,nc)
  c_wips = fill(wips,nc)
  pquad = lazy_map(evaluate,fgeomap,c_fips)

  if is_n_cube(p)
    c_detwips = c_wips
  elseif is_simplex(p)
    # Must account for diagonals in simplex discretizations to get the correct
    # scaling
    Jt1 = lazy_map(∇,fgeomap)
    Jt1_ips = lazy_map(evaluate,Jt1,c_fips)
    det_J = lazy_map(Broadcasting(meas),Jt1_ips)

    c_detwips = collect(lazy_map(Broadcasting(*),c_wips,det_J))
  end

  c_fips, pquad, c_detwips
end

function _broadcast(::Type{T},n,b) where T
  c = Array{T}(undef,size(b))
  for (ii, i) in enumerate(b)
    c[ii] = i⋅n
  end
  return c
end

function _RT_face_moments(p, fshfs, c_fips, fcips, fwips,phi)
  nc = length(c_fips)
  cfshfs = fill(fshfs, nc)
  cvals = lazy_map(evaluate,cfshfs,c_fips)
  cvals = [fwips[i].*cvals[i] for i in 1:nc]
  fns = get_facet_normal(p)

  # Must express the normal in terms of the real/reference system of
  # coordinates (depending if phi≡I or phi is a mapping, resp.)
  # Hence, J = transpose(grad(phi))

  Jt = fill(∇(phi),nc)
  Jt_inv = lazy_map(Operation(pinvJt),Jt)
  det_Jt = lazy_map(Operation(meas),Jt)
  change = lazy_map(*,det_Jt,Jt_inv)
  change_ips = lazy_map(evaluate,change,fcips)

  cvals = [ _broadcast(typeof(n),n,J.*b) for (n,b,J) in zip(fns,cvals,change_ips)]

  return cvals
end

# It provides for every face the nodes and the moments arrays
function _RT_face_values(p,et,order,phi)

  # Reference facet
  @assert is_simplex(p) || is_n_cube(p) "We are assuming that all n-faces of the same n-dim are the same."
  fp = Polytope{num_dims(p)-1}(p,1)

  # geomap from ref face to polytope faces
  fgeomap = _ref_face_to_faces_geomap(p,fp)

  # Nodes are integration points (for exact integration)
  # Thus, we define the integration points in the reference
  # face polytope (fips and wips). Next, we consider the
  # n-face-wise arrays of nodes in fp (constant cell array c_fips)
  # the one of the points in the polytope after applying the geopmap
  # (fcips), and the weights for these nodes (fwips, a constant cell array)
  # Nodes (fcips)
  degree = (order+1)*2
  fquad = Quadrature(fp,degree)
  fips = get_coordinates(fquad)
  wips = get_weights(fquad)

  c_fips, fcips, fwips = _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)

  # Moments (fmoments)
  # The RT prebasis is expressed in terms of shape function
  fshfs = MonomialBasis(et,fp,order)

  # Face moments, i.e., M(Fi)_{ab} = q_RF^a(xgp_RFi^b) w_Fi^b n_Fi ⋅ ()
  fmoments = _RT_face_moments(p, fshfs, c_fips, fcips, fwips, phi)

  return fcips, fmoments

end

function _RT_cell_moments(p, cbasis, ccips, cwips)
  # Interior DOFs-related basis evaluated at interior integration points
  ishfs_iips = evaluate(cbasis,ccips)
  return cwips.⋅ishfs_iips
end

_p_filter(e,order) = (sum(e) <= order)

# It provides for every cell the nodes and the moments arrays
function _RT_cell_values(p,et,order,phi)
  # Compute integration points at interior
  degree = 2*(order+1)
  iquad = Quadrature(p,degree)
  ccips = get_coordinates(iquad)
  cwips = get_weights(iquad)

  # Cell moments, i.e., M(C)_{ab} = q_C^a(xgp_C^b) w_C^b ⋅ ()
  if is_n_cube(p)
    cbasis = QGradMonomialBasis{num_dims(p)}(et,order-1)
  elseif is_simplex(p)
    T = VectorValue{num_dims(p),et}
    cbasis = MonomialBasis{num_dims(p)}(T,order-1, _p_filter)
  else
    @notimplemented
  end
  cell_moments = _RT_cell_moments(p, cbasis, ccips, cwips )

  # Must scale weights using phi map to get the correct integrals
  # scaling = meas(grad(phi))
  Jt = ∇(phi)
  Jt_inv = pinvJt(Jt)
  det_Jt = meas(Jt)
  change = det_Jt*Jt_inv
  change_ips = evaluate(change,ccips)

  cmoments = change_ips.⋅cell_moments

  return [ccips], [cmoments]

end

function _face_own_dofs_from_moments(f_moments)
  face_dofs = Vector{Int}[]
  o = 1
  for moments in f_moments
    ndofs = size(moments,2)
    dofs = collect(o:(o+ndofs-1))
    push!(face_dofs,dofs)
    o += ndofs
  end
  face_dofs
end

struct Moment <: Dof end

struct MomentBasedDofBasis{P,V} <: AbstractVector{Moment}
  nodes::Vector{P}
  face_moments::Vector{Array{V}}
  face_nodes::Vector{UnitRange{Int}}

  function MomentBasedDofBasis(nodes,f_moments,f_nodes)
    P = eltype(nodes)
    V = eltype(eltype(f_moments))
    new{P,V}(nodes,f_moments,f_nodes)
  end

  function MomentBasedDofBasis(f_nodes,f_moments)
    P = eltype(eltype(f_nodes))
    V = eltype(eltype(f_moments))
    nodes = P[]
    face_nodes = UnitRange{Int}[]
    nfaces = length(f_nodes)
    n = 1
    for fi in 1:nfaces
      nodes_fi = f_nodes[fi]
      nini = n
      for node_fi in nodes_fi
        push!(nodes,node_fi)
        n += 1
      end
      nend = n-1
      push!(face_nodes,nini:nend)
    end
    new{P,V}(nodes,f_moments,face_nodes)
  end
end

Base.size(a::MomentBasedDofBasis) = (length(a.nodes),)
Base.axes(a::MomentBasedDofBasis) = (axes(a.nodes,1),)
# @santiagobadia : Not sure we want to create the moment dofs
Base.getindex(a::MomentBasedDofBasis,i::Integer) = Moment()
Base.IndexStyle(::MomentBasedDofBasis) = IndexLinear()

get_nodes(b::MomentBasedDofBasis) = b.nodes
get_face_moments(b::MomentBasedDofBasis) = b.face_moments
get_face_nodes_dofs(b::MomentBasedDofBasis) = b.face_nodes

function num_dofs(b::MomentBasedDofBasis)
  n = 0
  for m in b.face_moments
    n += size(m,2)
  end
  n
end

function return_cache(b::MomentBasedDofBasis,field)
  cf = return_cache(field,b.nodes)
  vals = evaluate!(cf,field,b.nodes)
  ndofs = num_dofs(b)
  r = _moment_dof_basis_cache(vals,ndofs)
  c = CachedArray(r)
  (c, cf)
end

function _moment_dof_basis_cache(vals::AbstractVector,ndofs)
  T = eltype(vals)
  r = fill(zero(eltype(T))*0.,ndofs)
end

function _moment_dof_basis_cache(vals::AbstractMatrix,ndofs)
  _, npdofs = size(vals)
  T = eltype(vals)
  r = fill(zero(eltype(T))*0.,ndofs,npdofs)
end

function evaluate!(cache,b::MomentBasedDofBasis,field)
  c, cf = cache
  vals = evaluate!(cf,field,b.nodes)
  dofs = c.array
  _eval_moment_dof_basis!(dofs,vals,b)
  dofs
end

function _eval_moment_dof_basis!(dofs,vals::AbstractVector,b)
  o = 1
  z = zero(eltype(dofs))
  face_nodes = b.face_nodes
  face_moments = b.face_moments
  for face in 1:length(face_moments)
    moments = face_moments[face]
    if length(moments) != 0
      nodes = face_nodes[face]
      ni,nj = size(moments)
      for j in 1:nj
        dofs[o] = z
        for i in 1:ni
          dofs[o] += moments[i,j]⋅vals[nodes[i]]
        end
        o += 1
      end
    end
  end
end

function _eval_moment_dof_basis!(dofs,vals::AbstractMatrix,b)
  o = 1
  na = size(vals,2)
  z = zero(eltype(dofs))
  face_nodes = b.face_nodes
  face_moments = b.face_moments
  for face in 1:length(face_moments)
    moments = face_moments[face]
    if length(moments) != 0
      nodes = face_nodes[face]
      ni,nj = size(moments)
      for j in 1:nj
        for a in 1:na
          dofs[o,a] = z
          for i in 1:ni
            dofs[o,a] += moments[i,j]⋅vals[nodes[i],a]
          end
        end
        o += 1
      end
    end
  end
end

struct ContraVariantPiolaMap <: Map end

function evaluate!(
  cache,
  ::Broadcasting{typeof(∇)},
  a::Fields.BroadcastOpFieldArray{ContraVariantPiolaMap})
  v, Jt, detJ,sign_flip = a.args
  # Assuming J comes from an affine map
  ∇v = Broadcasting(∇)(v)
  k = ContraVariantPiolaMap()
  Broadcasting(Operation(k))(∇v,Jt,detJ,sign_flip)
end

function lazy_map(
  ::Broadcasting{typeof(gradient)},
  a::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}})
  v, Jt, detJ,sign_flip = a.args
  ∇v = lazy_map(Broadcasting(∇),v)
  k = ContraVariantPiolaMap()
  lazy_map(Broadcasting(Operation(k)),∇v,Jt,detJ,sign_flip)
end

function evaluate!(cache,::ContraVariantPiolaMap,
                   v::Number,
                   Jt::Number,
                   detJ::Number,
                   sign_flip::Bool)
  ((-1)^sign_flip*v)⋅((1/detJ)*Jt)
end

function evaluate!(cache,
                   k::ContraVariantPiolaMap,
                   v::AbstractVector{<:Field},
                   phi::Field,
                   sign_flip::AbstractVector{<:Field})
  Jt = ∇(phi)
  detJ = Operation(meas)(Jt)
  Broadcasting(Operation(k))(v,Jt,detJ,sign_flip)
end

function lazy_map(
  k::ContraVariantPiolaMap,
  cell_ref_shapefuns::AbstractArray{<:AbstractArray{<:Field}},
  cell_map::AbstractArray{<:Field},
  sign_flip::AbstractArray{<:AbstractArray{<:Field}})

  cell_Jt = lazy_map(∇,cell_map)
  cell_detJ = lazy_map(Operation(meas),cell_Jt)

  lazy_map(Broadcasting(Operation(k)),cell_ref_shapefuns,cell_Jt,cell_detJ,sign_flip)
end


text:
""" Creates a body-centered lattice with a single site """
bcc(T::Type=Float64; unit=u"nm") = Crystal(
    T[-0.5 0.5 0.5; 0.5 -0.5 0.5; 0.5 0.5 -0.5] * unit,
    position=T[0, 0, 0] * unit
)

""" Creates an face centered lattice with a single site """
fcc(T::Type=Float64; unit=u"nm") = Crystal(
    T[0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0] * unit,
    position=T[0, 0, 0] * unit
)


text:
"""
    plotabm(model::ABM{<: GraphSpace}; ac, as, am, kwargs...)
This function is the same as `plotabm` for `ContinuousSpace`, but here the three key
functions `ac, as, am` do not get an agent as an input but a vector of agents at
each node of the graph. Their output is the same.

Here `as` defaults to `length`. Internally, the `graphplot` recipe is used, and
all other `kwargs...` are propagated there.
"""
function plotabm(
    model::ABM{<:GraphSpace};
    ac = x -> "#765db4",
    as = length,
    am = x -> :circle,
    kwargs...,
)

    N = nodes(model)
    ncolor = Vector(undef, length(N))
    weights = zeros(length(N))
    markers = Vector(undef, length(N))
    for (i, n) in enumerate(N)
        a = get_node_agents(n, model)
        ncolor[i] = ac(a)
        weights[i] = as(a)
        markers[i] = am(a)
    end

    graphplot(
        model.space.graph,
        node_weights = weights,
        nodeshape = markers,
        nodecolor = ncolor,
        color = "black",
        markerstrokecolor = "black",
        markerstrokewidth = 1.5;
        kwargs...,
    )
end


text:
# This opens 14*3 = 42 websockets which close from the client side after
# timeouts up to 8192 seconds.
# Not included in runtests.jl because this takes two and a half hours to run.
# Checking against erros equires inspecting REPL.

using Test
using WebSockets
import Dates.now
include("../../benchmark/functions_open_browsers.jl")
include("timeout_functions.jl")

const WSIDLIST = Dict{typeof(time_ns()) , String}()
const WSLIST = Dict{typeof(time_ns()) , WebSocket}()
const CLIENTTASKS = Dict{typeof(time_ns()) , Task}()
const CLOSINGMESSAGES = Dict{typeof(time_ns()) , String}()
const CURRENT = Vector{Pair{WebSocket, String}}()
const T0 = now()
# At Http 0.7 / WebSockets 1.3, there is no obviously nice interface to limiting
# the maximum simultaneous connections. We settle for a naughty as well as ugly way, and state our
# intention to fix this up (sometime). It should not concern timeouts, just
# our ability to open many sockets at a time.
@eval WebSockets.HTTP.ConnectionPool default_connection_limit = 32

const SERVER = WebSockets.ServerWS(handle, gatekeeper, rate_limit = 1000//1)
const OLDLOGGER = WebSockets.global_logger()

WebSockets.global_logger(WebSocketLogger())
# Uncomment to include logging messages from HTTP
#WebSockets.global_logger(WebSocketLogger(shouldlog= (_, _, _, _, _) -> true))

const SERTASK = @async WebSockets.serve(SERVER, 8000)
open_a_browser()
# We'll also open a second type of browser. If all are available, the default sequence is
# ["chrome", "firefox", "iexplore", "safari", "phantomjs"]
open_a_browser()

# The browsers are working on opening client side websockets. Meanwhile,
# Julia will open some, too:
for i = 0:13
    sec = 2^i
    wsh = clientwsh(sec)
    push!(CLIENTTASKS, time_ns() => @async WebSockets.open(wsh, "ws://127.0.0.1:8000"))
    # We want to yield for asyncronous functions started by incoming
    # requests from browsers. Otherwise, the browsers could perhaps become bored.
    yield()
end

# The time to open a websocket may depend on things like the browser updating,
# or, for Julia, if compilation is necessary.
sleep(24)
if checktasks() == 14
    @wslog "All client tasks are running without error"
end
@test length(WSLIST) == 42
for (key, ws) in WSLIST
    push!(CURRENT, ws => get(WSIDLIST, key, "See WSIDLIST directly"))
end
@wslog CURRENT
@async begin
    # Wait for all timeouts (the longest is 8192s)
    @wslog "\e[32m --- The final tests will run in 8250s---\n
    \t For more viewing pleasure, now inspect CURRENT\n
    \t\e[39m or alternatively WSIDLIST WSLIST  CLIENTTASKS  CLOSINGMESSAGES and SERVER"
    sleep(8250)
    #put!(SERVER.in("Job well done!"))
    for (key, ws) in WSLIST
        @wslog ws, WSIDLIST[key]
        @test !WebSockets.isopen(ws)
    end

    close(SERVER)

    #Check that all the sockets were closed for the right reason
    for clmsg in values(CLOSINGMESSAGES)
        @test occursin("seconds are up", clmsg)
    end
    WebSockets.global_logger(OLDLOGGER)
end
nothing


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   API functions for all the orbit representations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export get_epoch, get_a, get_e, get_i, get_Ω, get_ω, get_f, get_raan, get_argp,
       get_M, get_r, get_v, get_rv

"""
    get_epoch(orb)

Return the epoch of the orbit representation `orb` [m].

"""
get_epoch(orb)

"""
    get_a(orb)

Return the semi-major axis of the orbit representation `orb` [m].

"""
get_a

"""
    get_e(orb)

Return the eccentricity of the orbit representation `orb`.

"""
get_e

"""
    get_i(orb)

Return the inclination of the orbit representation `orb` [rad].

"""
get_i

"""
    get_Ω(orb)

Return the right ascention of the ascending node (RAAN) of the orbit
representation `orb` [rad].

"""
get_Ω

"""
    get_ω(orb)

Return the argument of periapsis of the representation `orb` [rad].

"""
get_ω

"""
    get_f(orb)

Return the true anomaly of the representation `orb` [rad].

"""
get_f

"""
    get_raan(orb)

Return the right ascention of the ascending node (RAAN) of the orbit
representation `orb` [rad].

"""
@inline get_raan(orb) = get_Ω(orb)

"""
    get_argp(orb)

Return the argument of periapsis of the representation `orb` [rad].

"""
@inline get_argp(orb) = get_ω(orb)

"""
    get_M(orb)

Return the mean anomaly of the representation `orb` [rad].

"""
get_M

"""
    get_r(orb)

Return the position vector of the representation `orb` [m].

"""
get_r

"""
    get_r(orb)

Return the velocity vector of the representation `orb` [m].

"""
get_v

"""
    get_rv(orb)

Return the position and velocity vector of the representation `orb` [m].

"""
get_rv


text:
"""
    format_version(v::Union{VersionNumber, AbstractString}) -> String

Strip everything but the major and minor release from a `VersionNumber`.
Strings are left in their original form.
"""
format_version(v::VersionNumber) = "$(v.major).$(v.minor)"
format_version(v::AbstractString) = string(v)

const ALLOWED_FAILURES = ["nightly"]  # TODO: Update this list with new RCs.
const DEFAULT_CI_VERSIONS = map(format_version, [default_version(), VERSION, "nightly"])
const DEFAULT_CI_VERSIONS_NO_NIGHTLY = map(format_version, [default_version(), VERSION])
const EXTRA_VERSIONS_DOC = "- `extra_versions::Vector`: Extra Julia versions to test, as strings or `VersionNumber`s."

"""
    GitHubActions(;
        file="$(contractuser(default_file("github", "workflows", "CI.yml")))",
        destination="CI.yml",
        linux=true,
        osx=false,
        windows=false,
        x64=true,
        x86=false,
        coverage=true,
        extra_versions=$DEFAULT_CI_VERSIONS,
    )

Integrates your packages with [GitHub Actions](https://github.com/features/actions).

## Keyword Arguments
- `file::AbstractString`: Template file for the workflow file.
- `destination::AbstractString`: Destination of the workflow file,
  relative to `.github/workflows`.
- `linux::Bool`: Whether or not to run builds on Linux.
- `osx::Bool`: Whether or not to run builds on OSX (MacOS).
- `windows::Bool`: Whether or not to run builds on Windows.
- `x64::Bool`: Whether or not to run builds on 64-bit architecture.
- `x86::Bool`: Whether or not to run builds on 32-bit architecture.
- `coverage::Bool`: Whether or not to publish code coverage.
  Another code coverage plugin such as [`Codecov`](@ref) must also be included.
$EXTRA_VERSIONS_DOC

!!! note
    If using coverage plugins, don't forget to manually add your API tokens as secrets,
    as described [here](https://help.github.com/en/actions/automating-your-workflow-with-github-actions/creating-and-using-encrypted-secrets#creating-encrypted-secrets).
"""
@plugin struct GitHubActions <: FilePlugin
    file::String = default_file("github", "workflows", "CI.yml")
    destination::String = "CI.yml"
    linux::Bool = true
    osx::Bool = false
    windows::Bool = false
    x64::Bool = true
    x86::Bool = false
    coverage::Bool = true
    extra_versions::Vector = DEFAULT_CI_VERSIONS
end

source(p::GitHubActions) = p.file
destination(p::GitHubActions) = joinpath(".github", "workflows", p.destination)
tags(::GitHubActions) = "<<", ">>"

badges(::GitHubActions) = Badge(
    "Build Status",
    "https://github.com/{{{USER}}}/{{{PKG}}}.jl/workflows/CI/badge.svg",
    "https://github.com/{{{USER}}}/{{{PKG}}}.jl/actions",
)

function view(p::GitHubActions, t::Template, pkg::AbstractString)
    os = String[]
    p.linux && push!(os, "ubuntu-latest")
    p.osx && push!(os, "macOS-latest")
    p.windows && push!(os, "windows-latest")
    arch = filter(a -> getfield(p, Symbol(a)), ["x64", "x86"])
    excludes = Dict{String, String}[]
    p.osx && p.x86 && push!(excludes, Dict("E_OS" => "macOS-latest", "E_ARCH" => "x86"))

    return Dict(
        "ARCH" => arch,
        "EXCLUDES" => excludes,
        "HAS_CODECOV" => p.coverage && hasplugin(t, Codecov),
        "HAS_COVERALLS" => p.coverage && hasplugin(t, Coveralls),
        "HAS_DOCUMENTER" => hasplugin(t, Documenter{GitHubActions}),
        "HAS_EXCLUDES" => !isempty(excludes),
        "OS" => os,
        "PKG" => pkg,
        "USER" => t.user,
        "VERSIONS" => collect_versions(t, p.extra_versions),
    )
end

"""
    TravisCI(;
        file="$(contractuser(default_file("travis.yml")))",
        linux=true,
        osx=false,
        windows=false,
        x64=true,
        x86=false,
        arm64=false,
        coverage=true,
        extra_versions=$DEFAULT_CI_VERSIONS,
    )

Integrates your packages with [Travis CI](https://travis-ci.com).

## Keyword Arguments
- `file::AbstractString`: Template file for `.travis.yml`.
- `linux::Bool`: Whether or not to run builds on Linux.
- `osx::Bool`: Whether or not to run builds on OSX (MacOS).
- `windows::Bool`: Whether or not to run builds on Windows.
- `x64::Bool`: Whether or not to run builds on 64-bit architecture.
- `x86::Bool`: Whether or not to run builds on 32-bit architecture.
- `arm64::Bool`: Whether or not to run builds on the ARM64 architecture.
- `coverage::Bool`: Whether or not to publish code coverage.
  Another code coverage plugin such as [`Codecov`](@ref) must also be included.
$EXTRA_VERSIONS_DOC
"""
@plugin struct TravisCI <: FilePlugin
    file::String = default_file("travis.yml")
    linux::Bool = true
    osx::Bool = false
    windows::Bool = false
    x64::Bool = true
    x86::Bool = false
    arm64::Bool = false
    coverage::Bool = true
    extra_versions::Vector = DEFAULT_CI_VERSIONS
end

source(p::TravisCI) = p.file
destination(::TravisCI) = ".travis.yml"

badges(::TravisCI) = Badge(
    "Build Status",
    "https://travis-ci.com/{{{USER}}}/{{{PKG}}}.jl.svg?branch=master",
    "https://travis-ci.com/{{{USER}}}/{{{PKG}}}.jl",
)

function view(p::TravisCI, t::Template, pkg::AbstractString)
    os = filter(o -> getfield(p, Symbol(o)), ["linux", "osx", "windows"])
    arch = filter(a -> getfield(p, Symbol(a)), ["x64", "x86", "arm64"])
    versions = collect_versions(t, p.extra_versions)
    allow_failures = filter(in(versions), ALLOWED_FAILURES)

    excludes = Dict{String, String}[]
    p.x86 && p.osx && push!(excludes, Dict("E_OS" => "osx", "E_ARCH" => "x86"))
    if p.arm64
        p.osx && push!(excludes, Dict("E_OS" => "osx", "E_ARCH" => "arm64"))
        p.windows && push!(excludes, Dict("E_OS" => "windows", "E_ARCH" => "arm64"))
        "nightly" in versions && push!(excludes, Dict("E_JULIA" => "nightly", "E_ARCH" => "arm64"))
    end

    return Dict(
        "ALLOW_FAILURES" => allow_failures,
        "ARCH" => arch,
        "EXCLUDES" => excludes,
        "HAS_ALLOW_FAILURES" => !isempty(allow_failures),
        "HAS_CODECOV" => hasplugin(t, Codecov),
        "HAS_COVERAGE" => p.coverage && hasplugin(t, is_coverage),
        "HAS_COVERALLS" => hasplugin(t, Coveralls),
        "HAS_DOCUMENTER" => hasplugin(t, Documenter{TravisCI}),
        "HAS_EXCLUDES" => !isempty(excludes),
        "OS" => os,
        "PKG" => pkg,
        "USER" => t.user,
        "VERSION" => format_version(t.julia),
        "VERSIONS" => versions,
    )
end

"""
    AppVeyor(;
        file="$(contractuser(default_file("appveyor.yml")))",
        x86=false,
        coverage=true,
        extra_versions=$DEFAULT_CI_VERSIONS,
    )

Integrates your packages with [AppVeyor](https://appveyor.com)
via [AppVeyor.jl](https://github.com/JuliaCI/Appveyor.jl).

## Keyword Arguments
- `file::AbstractString`: Template file for `.appveyor.yml`.
- `x86::Bool`: Whether or not to run builds on 32-bit systems,
  in addition to the default 64-bit builds.
- `coverage::Bool`: Whether or not to publish code coverage.
  [`Codecov`](@ref) must also be included.
$EXTRA_VERSIONS_DOC
"""
@plugin struct AppVeyor <: FilePlugin
    file::String = default_file("appveyor.yml")
    x86::Bool = false
    coverage::Bool = true
    extra_versions::Vector = DEFAULT_CI_VERSIONS
end

source(p::AppVeyor) = p.file
destination(::AppVeyor) = ".appveyor.yml"

badges(::AppVeyor) = Badge(
    "Build Status",
    "https://ci.appveyor.com/api/projects/status/github/{{{USER}}}/{{{PKG}}}.jl?svg=true",
    "https://ci.appveyor.com/project/{{{USER}}}/{{{PKG}}}-jl",
)

function view(p::AppVeyor, t::Template, pkg::AbstractString)
    platforms = ["x64"]
    p.x86 && push!(platforms, "x86")

    versions = collect_versions(t, p.extra_versions)
    allow_failures = filter(in(versions), ALLOWED_FAILURES)

    return Dict(
        "ALLOW_FAILURES" => allow_failures,
        "HAS_ALLOW_FAILURES" => !isempty(allow_failures),
        "HAS_CODECOV" => p.coverage && hasplugin(t, Codecov),
        "PKG" => pkg,
        "PLATFORMS" => platforms,
        "USER" => t.user,
        "VERSIONS" => versions,
    )
end

"""
    CirrusCI(;
        file="$(contractuser(default_file("cirrus.yml")))",
        image="freebsd-12-0-release-amd64",
        coverage=true,
        extra_versions=$DEFAULT_CI_VERSIONS,
    )

Integrates your packages with [Cirrus CI](https://cirrus-ci.org)
via [CirrusCI.jl](https://github.com/ararslan/CirrusCI.jl).

## Keyword Arguments
- `file::AbstractString`: Template file for `.cirrus.yml`.
- `image::AbstractString`: The FreeBSD image to be used.
- `coverage::Bool`: Whether or not to publish code coverage.
  [`Codecov`](@ref) must also be included.
$EXTRA_VERSIONS_DOC

!!! note
    Code coverage submission from Cirrus CI is not yet supported by
    [Coverage.jl](https://github.com/JuliaCI/Coverage.jl).
"""
@plugin struct CirrusCI <: FilePlugin
    file::String = default_file("cirrus.yml")
    image::String = "freebsd-12-0-release-amd64"
    coverage::Bool = true
    extra_versions::Vector = DEFAULT_CI_VERSIONS
end

source(p::CirrusCI) = p.file
destination(::CirrusCI) = ".cirrus.yml"

badges(::CirrusCI) = Badge(
    "Build Status",
    "https://api.cirrus-ci.com/github/{{{USER}}}/{{{PKG}}}.jl.svg",
    "https://cirrus-ci.com/github/{{{USER}}}/{{{PKG}}}.jl",
)

function view(p::CirrusCI, t::Template, pkg::AbstractString)
    return Dict(
        "HAS_CODECOV" => hasplugin(t, Codecov),
        "HAS_COVERALLS" => hasplugin(t, Coveralls),
        "HAS_COVERAGE" => p.coverage && hasplugin(t, is_coverage),
        "IMAGE" => p.image,
        "PKG" => pkg,
        "USER" => t.user,
        "VERSIONS" => collect_versions(t, p.extra_versions),
    )
end

"""
    GitLabCI(;
        file="$(contractuser(default_file("gitlab-ci.yml")))",
        coverage=true,
        extra_versions=$DEFAULT_CI_VERSIONS_NO_NIGHTLY,
    )

Integrates your packages with [GitLab CI](https://docs.gitlab.com/ce/ci).

## Keyword Arguments
- `file::AbstractString`: Template file for `.gitlab-ci.yml`.
- `coverage::Bool`: Whether or not to compute code coverage.
$EXTRA_VERSIONS_DOC

## GitLab Pages
Documentation can be generated by including a `Documenter{GitLabCI}` plugin.
See [`Documenter`](@ref) for more information.

!!! note
    Nightly Julia is not supported.
"""
@plugin struct GitLabCI <: FilePlugin
    file::String = default_file("gitlab-ci.yml")
    coverage::Bool = true
    # Nightly has no Docker image.
    extra_versions::Vector = DEFAULT_CI_VERSIONS_NO_NIGHTLY
end

gitignore(p::GitLabCI) = p.coverage ? COVERAGE_GITIGNORE : String[]
source(p::GitLabCI) = p.file
destination(::GitLabCI) = ".gitlab-ci.yml"

function badges(p::GitLabCI)
    ci = Badge(
        "Build Status",
        "https://{{{HOST}}}/{{{USER}}}/{{{PKG}}}.jl/badges/master/pipeline.svg",
        "https://{{{HOST}}}/{{{USER}}}/{{{PKG}}}.jl/pipelines",
    )
    cov = Badge(
        "Coverage",
        "https://{{{HOST}}}/{{{USER}}}/{{{PKG}}}.jl/badges/master/coverage.svg",
        "https://{{{HOST}}}/{{{USER}}}/{{{PKG}}}.jl/commits/master",
    )
    return p.coverage ? [ci, cov] : [ci]
end

function view(p::GitLabCI, t::Template, pkg::AbstractString)
    return Dict(
        "HAS_COVERAGE" => p.coverage,
        "HAS_DOCUMENTER" => hasplugin(t, Documenter{GitLabCI}),
        "HOST" => t.host,
        "PKG" => pkg,
        "USER" => t.user,
        "VERSION" => format_version(t.julia),
        "VERSIONS" => collect_versions(t, p.extra_versions),
    )
end

"""
    DroneCI(;
        file="$(contractuser(default_file("drone.star")))",
        amd64=true,
        arm=false,
        arm64=false,
        extra_versions=$DEFAULT_CI_VERSIONS_NO_NIGHTLY,
    )

Integrates your packages with [Drone CI](https://drone.io).

## Keyword Arguments
- `file::AbstractString`: Template file for `.drone.star`.
- `destination::AbstractString`: File destination, relative to the repository root.
  For example, you might want to generate a `.drone.yml` instead of the default Starlark file.
- `amd64::Bool`: Whether or not to run builds on AMD64.
- `arm::Bool`: Whether or not to run builds on ARM (32-bit).
- `arm64::Bool`: Whether or not to run builds on ARM64.
$EXTRA_VERSIONS_DOC

!!! note
    Nightly Julia is not supported.
"""
@plugin struct DroneCI <: FilePlugin
    file::String = default_file("drone.star")
    destination::String = ".drone.star"
    amd64::Bool = true
    arm::Bool = false
    arm64::Bool = false
    extra_versions::Vector = DEFAULT_CI_VERSIONS_NO_NIGHTLY
end

source(p::DroneCI) = p.file
destination(p::DroneCI) = p.destination

badges(::DroneCI) = Badge(
    "Build Status",
    "https://cloud.drone.io/api/badges/{{{USER}}}/{{{PKG}}}.jl/status.svg",
    "https://cloud.drone.io/{{{USER}}}/{{{PKG}}}.jl",
)

function view(p::DroneCI, t::Template, pkg::AbstractString)
    arches = String[]
    p.amd64 && push!(arches, "amd64")
    p.arm && push!(arches, "arm")
    p.arm64 && push!(arches, "arm64")

    return Dict(
        "ARCHES" => join(map(repr, arches), ", "),
        "PKG" => pkg,
        "USER" => t.user,
        "VERSIONS" => join(map(repr, collect_versions(t, p.extra_versions)), ", "),
    )
end

"""
    collect_versions(t::Template, versions::Vector) -> Vector{String}

Combine `t`'s Julia version with `versions`, and format them as `major.minor`.
This is useful for creating lists of versions to be included in CI configurations.
"""
function collect_versions(t::Template, versions::Vector)
    custom = map(v -> v isa VersionNumber ? format_version(v) : string(v), versions)
    vs = map(v -> lstrip(v, 'v'), [format_version(t.julia); custom])
    filter!(vs) do v
        # Throw away any versions lower than the template's minimum.
        try
            VersionNumber(v) >= t.julia
        catch e
            e isa ArgumentError || rethrow()
            true
        end
    end
    return sort(unique(vs))
end

const AllCI = Union{AppVeyor, GitHubActions, TravisCI, CirrusCI, GitLabCI, DroneCI}

"""
    is_ci(::Plugin) -> Bool

Determine whether or not a plugin is a CI plugin.
If you are adding a CI plugin, you should implement this function and return `true`.
"""
is_ci(::Plugin) = false
is_ci(::AllCI) = true

needs_username(::AllCI) = true
customizable(::Type{<:AllCI}) = (:extra_versions => Vector{String},)