text:
using SFGTools
using Test
using DataFrames
using SFGTools
using Random
import Statistics: mean

const SAMPLE_DATA_DIR = joinpath(@__DIR__, "sampledata/")

function listtest()
    grab(SAMPLE_DATA_DIR; getall=true)
    grab(SAMPLE_DATA_DIR)
    df = list_spectra()
    @test size(df,1) == 211
    df = list_spectra(date=(2020,5,20))
    @test size(df,1) == 154
    df = list_spectra(inexact="ssp")
    @test size(df,1) == 209
    df = list_spectra(exact="200518_CaAra_1_13mN_d-pumped_ssp_DLscan")
    @test size(df,1) == 150
    df = list_spectra(group=true)
    @test size(df,1) == 6
end

function loadtest()
    spectrum = load_spectra(63725661936614)
    @test typeof(spectrum) == Array{SFGTools.SFSpectrum,1}
    @test spectrum[1][1] == 778.0
    @test_broken load_spectra(63725661936614, format=:tiff) ==
          load_spectra(63725661936614, format=:sif)
    spectrum = load_spectra(63725661936614)
    @show typeof(spectrum)
    spectrum
end

function attribute_test(spectrum)
    meta = get_metadata(spectrum)
    @test typeof(meta) == Dict{String,Any}
    attr = get_attribute(spectrum, "ccd_exposure_time")
    @test attr == Any[20.0, 20.0, 20.0, 20.0, 20.0]
    wn = get_ir_wavenumber(spectrum)
    @test typeof(wn) == Array{Float64,1}
end

function fieldcorrection_test()
    a = Array{Float64,3}(undef, (4, 1, 2))
    a[:,1,1] = [6, 7, 8, 9]
    a[:,1,2] = [6, 8, 9, 7]
    spectrum = SFSpectrum(0, a)

    b = Array{Float64,3}(undef, (4, 1, 2))
    b[:,1,1] = [0, 1, 3, 0]
    b[:,1,2] = [2, 1, 1, 0]
    bias = SFSpectrum(0, b)

    d = Array{Float64,3}(undef, (4, 1, 2))
    d[:,1,1] = [3, 4, 3, 5]
    d[:,1,2] = [4, 2, 1, 5]
    dark = SFSpectrum(0, d)

    f = Array{Float64,3}(undef, (4, 1, 2))
    f[:,1,1] = [2, 3, 4, 3]
    f[:,1,2] = [4, 3, 6, 1]
    flat = SFSpectrum(0, f)

    l = Array{Float64,3}(undef, (4, 1, 2))
    l[:,1,1] = [1, 2, 3, 0]
    l[:,1,2] = [3, 2, 3, 2]
    darkflat = SFSpectrum(0, l)

    fieldcorrection!(spectrum, bias=bias, dark=dark, flat=flat, darkflat=darkflat)

    @test spectrum[:,1,1] == [5.0, 8.0,  6.0, 8.0]
    @test spectrum[:,1,2] == [5.0, 10.0, 7.0, 4.0]

    fieldcorrection!(spectrum, dark=dark)

    @test spectrum[:,1,1] == [1.5, 5.0, 4.0, 3.0]
    @test spectrum[:,1,2] == [1.5, 7.0, 5.0, -1.0]
end

function rm_events_test()
    a = rand(MersenneTwister(0), 50)
    a[19:20] .*= 5
    a[22] *= 40
    num_removed_events = rm_events!(a, minstd=2)
    @test num_removed_events == 2
end

function average_test()
    a = Array{Float64,3}(undef, (4, 1, 2))
    a[:,1,1] = [6, 7, 8, 9]
    a[:,1,2] = [6, 8, 9, 7]
    spectrum1 = SFSpectrum(63725661936614, a)

    b = Array{Float64,3}(undef, (4, 1, 3))
    b[:,1,1] = [0, 1, 3, 0]
    b[:,1,2] = [2, 1, 1, 0]
    b[:,1,3] = [4, 1, 2, 0]
    spectrum2 = SFSpectrum(63725661936614, b)

    spectrum = average(spectrum1)
    @test_broken spectrum[:,1,1] == [6.0, 7.5,  8.5, 8.0]

    spectrum = average(spectrum1, combine = false)
    @test_broken spectrum[:,1,1] == [6, 7, 8, 9]
    @test_broken spectrum[:,1,2] == [6, 8, 9, 7]

    spectrum = average(spectrum2)
    @test_broken spectrum[:,1,1] == [2.0, 1.0, 2.0, 0]

    spectrum = average(spectrum2, combine = false)
    @test_broken spectrum[:,1,1] == [0, 1, 3, 0]
    @test_broken spectrum[:,1,2] == [2, 1, 1, 0]
    @test_broken spectrum[:,1,3] == [4, 1, 2, 0]
end

# function save_mat_test(spectra)
#     success = false
#     try
#         save_mat(tempname(), spectra)
#         success = true
#     catch
#         success = false
#     end
#     @test success == true
# end

function makespectraarray(spectrum::SFSpectrum)
    spectra = Array{SFSpectrum,1}(undef, size(spectrum,2))
    for i = 1:size(spectrum, 2)
        spectra[i] = SFSpectrum(spectrum.id, spectrum[:,i,1])
    end
    spectra
end

@testset "list_spectra Tests" begin listtest() end
@testset "load_spectra Tests" begin global spectrum = loadtest()[1] end
@testset "Attribute Tests" begin attribute_test(spectrum) end
@testset "Fieldcorrection Tests" begin fieldcorrection_test() end
@testset "Event Removal Tests" begin rm_events_test() end
@testset "average Test" begin average_test() end
# spectra = makespectraarray(spectrum)
# @testset "MAT Saving" begin save_mat_test(spectra) end
df


text:
lines = open(readlines, "FileCG.txt")
newlines = sort(lines)
writedlm("FileCG.txt",newlines)

lines = open(readlines, "FileCR.txt")
newlines = sort(lines)
writedlm("FileCR.txt",newlines)

ext: jl


text:
function _promote(args...)
    _type = Base.promote_eltype(args...)
    return map(x->convert.(_type, x), args)
end

## Utilities

"""
$(SIGNATURES)

Check if the problem has control inputs.
"""
is_autonomous(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = U

"""
$(SIGNATURES)

Check if the problem is time discrete.
"""
is_discrete(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = C == DDProbType(2)

"""
$(SIGNATURES)

Check if the problem is direct.
"""
is_direct(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = C == DDProbType(1)

"""
$(SIGNATURES)

Check if the problem is time continuous.
"""
is_continuous(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = C == DDProbType(3)

"""
$(SIGNATURES)

Check if the problem is parameterized.
"""
is_parametrized(x::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = hasfield(typeof(x), :p) && !isempty(x.p)

"""
$(SIGNATURES)

Check if the problem has associated measurement times.
"""
has_timepoints(x::AbstractDataDrivenProblem{N,U,C}) where {N,U,C} = !isempty(x.t)

## Concrete Type
"""
$(TYPEDEF)

The `DataDrivenProblem` defines a general estimation problem given measurements, inputs and (in the near future) observations.
Two construction methods are available:

+ `DiscreteDataDrivenProblem` for time discrete systems
+ `ContinousDataDrivenProblem` for systems continuous in time

both are aliases for constructing a problem.

# Fields
$(FIELDS)

# Signatures
$(SIGNATURES)

# Example

```julia
X, DX, t = data...

# Define a discrete time problem
prob = DiscreteDataDrivenProblem(X)

# Define a continuous time problem without explicit time points
prob = ContinuousDataDrivenProblem(X, DX)

# Define a continuous time problem without explicit derivatives
prob = ContinuousDataDrivenProblem(X, t)

# Define a discrete time problem with an input function as a function
input_signal(u,p,t) = t^2
prob = DiscreteDataDrivenProblem(X, t, input_signal)
```
"""
struct DataDrivenProblem{dType, cType, probType} <: AbstractDataDrivenProblem{dType, cType, probType}

    # Data
    """State measurements"""
    X::AbstractMatrix{dType}
    """Time measurements (optional)"""
    t::AbstractVector{dType}

    """Differental state measurements (optional); Used for time continuous problems"""
    DX::AbstractMatrix{dType}
    """Output measurements (optional); Used for direct problems"""
    Y::AbstractMatrix{dType}
    """Input measurements (optional); Used for non-autonoumous problems"""
    U::AbstractMatrix{dType}

    """Parameters associated with the problem (optional)"""
    p::AbstractVector{dType}

    """Name of the problem"""
    name::Symbol
end

function DataDrivenProblem(probType, X, t, DX, Y, U, p; name = gensym(:DDProblem), kwargs...)
    dType = Base.promote_eltype(X, t, DX, Y, U, p)
    cType = isempty(U)
    name = isa(name, Symbol) ? name : Symbol(name)
    # We assume a discrete Problem
    if isnothing(probType)
        probType = DDProbType(2)
        if (isempty(DX) && !isempty(Y))
            probType = DDProbType(1) # Direct problem
        elseif !isempty(DX)
            probType = DDProbType(3) # Continouos
        end
    end

    return DataDrivenProblem{dType, cType, probType}(_promote(X,t,DX,Y,U,p)..., name)
end

function DataDrivenProblem(probtype, X, t, DX, Y, U::F, p; kwargs...) where F <: Function
    # Generate the input as a Matrix

    ts = isempty(t) ? zeros(eltype(X),size(X,2)) : t

    u_ = hcat(map(i->U(X[:,i], p, ts[i]), 1:size(X,2))...)

    return DataDrivenProblem(probtype, _promote(X,t,DX,Y,u_,p)...; kwargs...)
end

function DataDrivenProblem(X::AbstractMatrix;
    t::AbstractVector = collect(one(eltype(X)):size(X,2)),
    DX::AbstractMatrix = Array{eltype(X)}(undef, 0, 0),
    Y::AbstractMatrix = Array{eltype(X)}(undef, 0,0),
    U::F = Array{eltype(X)}(undef, 0,0),
    p::AbstractVector = Array{eltype(X)}(undef, 0),
    probtype = nothing,
    kwargs...) where F <: Union{AbstractMatrix, Function}

    return DataDrivenProblem(probtype, X,t,DX,Y,U,p; kwargs...)
end

function Base.summary(io::IO, x::DataDrivenProblem{N,C,P}) where {N,C,P}
    print(io, "$P DataDrivenProblem{$N} $(x.name)")
    n,m = size(x.X)
    print(io, " in $n dimensions and $m samples")
    C ? nothing : print(io, " with controls")
    return
end

function Base.print(io::IO, x::AbstractDataDrivenProblem{N,C,P}) where {N,C,P}
    println(io, "$P DataDrivenProblem{$N} $(x.name)")
    println(io, "Summary")
    n,m = size(x.X)
    println(io, "$m measurements")
    println(io, "$n state(s)")
    n = size(x.DX, 1)
    isempty(x.DX) ? nothing : println("$n differential state(s)")
    n = size(x.Y, 1)
    isempty(x.Y) ? nothing : println("$n observed variable(s)")
    n = size(x.U, 1)
    isempty(x.U) ? nothing : println(io, "$n control(s)")
end

Base.show(io::IO, x::DataDrivenProblem{N,C,P}) where {N,C,P} = summary(io, x)

## Discrete Constructors
"""
A time discrete `DataDrivenProblem` useable for problems of the form `f(x[i],p,t,u) ↦ x[i+1]`.

$(SIGNATURES)
"""
function DiscreteDataDrivenProblem(X::AbstractMatrix; kwargs...)
    DataDrivenProblem(X; probtype = DDProbType(2), kwargs...)
end

function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector; kwargs...)
    DataDrivenProblem(X; t=t,  probtype = DDProbType(2),  kwargs...)
end

function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, U::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X; t=t, U = U , probtype = DDProbType(2),  kwargs...)
end

function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, U::Function; kwargs...)
    return DataDrivenProblem(X; t=t, U = U ,  probtype = DDProbType(2), kwargs...)
end

## Continouos Constructors
"""
A time continuous `DataDrivenProblem` useable for problems of the form `f(x,p,t,u) ↦ dx/dt`.

$(SIGNATURES)

Automatically constructs derivatives via an additional collocation method, which can be either a collocation
or an interpolation from `DataInterpolations.jl` wrapped by an `InterpolationMethod`.
"""
function ContinuousDataDrivenProblem(X::AbstractMatrix, DX::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X; DX = DX,  probtype = DDProbType(3), kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, DX::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X; t = t, DX = DX, probtype = DDProbType(3), kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, DX::AbstractMatrix, U::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X; t = t, DX = DX, U = U, probtype = DDProbType(3), kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, DX::AbstractMatrix, U::F; kwargs...) where {F <: Function}
    return DataDrivenProblem(X; t = t, DX = DX, U = U, probtype = DDProbType(3), kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, collocation = InterpolationMethod(); kwargs...)
    dx, x = collocate_data(X, t, collocation)
    return DataDrivenProblem(x; t = t, DX = dx, probtype = DDProbType(3), kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector,  U::AbstractMatrix, collocation; kwargs...)
    dx, x = collocate_data(X, t, collocation)
    return DataDrivenProblem(x; t = t, DX = dx, U = U, probtype = DDProbType(3), kwargs...)
end

## Direct Constructors
"""
A direct `DataDrivenProblem` useable for problems of the form `f(x,p,t,u) ↦ y`.

$(SIGNATURES)
"""
function DirectDataDrivenProblem(X::AbstractMatrix, Y::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X; Y = Y, probtype = DDProbType(1), kwargs...)
end

function DirectDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, Y::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X; t = t, Y = Y, probtype = DDProbType(1), kwargs...)
end

function DirectDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, Y::AbstractMatrix, U; kwargs...)
    return DataDrivenProblem(X; t = t, Y = Y, U = U, probtype = DDProbType(1),kwargs...)
end

Base.length(prob::AbstractDataDrivenProblem) = size(prob.X, 2)
# Special case of Discrete
Base.length(prob::AbstractDiscreteProb) = size(prob.X, 2)-1

Base.size(prob::AbstractDataDrivenProblem) = size(prob.X)
Base.size(prob::AbstractDiscreteProb) = begin
    n,m = size(prob.X)
    return (n, m-1)
end

"""
$(SIGNATURES)

Returns the name of the `problem`.
"""
get_name(p::AbstractDataDrivenProblem) = getfield(p, :name)

## Utils

function ModelingToolkit.states(p::AbstractDataDrivenProblem, i = :, j = :)
    x = getfield(p, :X)
    isempty(x) ? x : getindex(x, i, j)
end

function ModelingToolkit.parameters(p::AbstractDataDrivenProblem, i = :)
    x = getfield(p, :p)
    isempty(x) ? x : getindex(x, i)
end

function ModelingToolkit.independent_variable(p::AbstractDataDrivenProblem, i = :)
    x = getfield(p, :t)
    isempty(x) ? x : getindex(x, i)
end

function ModelingToolkit.get_dvs(p::AbstracContProb, i = :, j = :)
    x = getfield(p, :DX)
    isempty(x) ? x : getindex(x, i, j)
end

function ModelingToolkit.get_dvs(p::AbstractDirectProb, i = :, j = :)
    x = getfield(p, :Y)
    isempty(x) ? x : getindex(x, i, j)
end

function ModelingToolkit.get_dvs(p::AbstractDiscreteProb, i = :, j = :)
    ModelingToolkit.states(p, i, j)
end

function ModelingToolkit.observed(p::AbstractDataDrivenProblem, i = :, j = :)
    x = getfield(p, :Y)
    isempty(x) ? x : getindex(x, i, j)
end

function ModelingToolkit.controls(p::AbstractDataDrivenProblem, i = :, j = :)
    x = getfield(p, :U)
    isempty(x) ? x : getindex(x, i, j)
end

function Base.getindex(p::AbstractDataDrivenProblem, i = :, j = :)
    return (
        ModelingToolkit.states(p, i, j),
        ModelingToolkit.parameters(p),
        ModelingToolkit.independent_variable(p, j),
        ModelingToolkit.controls(p, i, j)
    )
end

# TODO This is just for explicit basis!
# Make the basis callable with the problem
# Explicit

@views (b::AbstractBasis)(p::AbstractDataDrivenProblem) = begin
    b(states(p), parameters(p), independent_variable(p), controls(p))
end

@views (b::AbstractBasis)(dx::AbstractMatrix, p::AbstractDataDrivenProblem) = begin
    b(dx, states(p), parameters(p), independent_variable(p), controls(p))
end

@views (b::AbstractBasis)(p::AbstractDataDrivenProblem, j) = begin
    b(p[:, j]...)
end

@views (b::AbstractBasis)(dx::AbstractMatrix, p::AbstractDataDrivenProblem, j) = begin
    b(dx, p[:, j]...)
end

# Special case discrete problem

@views (b::AbstractBasis)(p::AbstractDiscreteProb) = begin
    b(p[:, 1:length(p)]...)
end

@views (b::AbstractBasis)(dx::AbstractMatrix, p::AbstractDiscreteProb) = begin
    b(dx, p, 1:length(p))
end

# Check for nans, inf etc
check_domain(x) =  @assert all(.~isnan.(x)) && all(.~isinf.(x)) ("One or more measurements contain `NaN` or `Inf`.")
check_lengths(args...) = @assert all(map(x->size(x)[end], args) .== size(first(args))[end]) "One or more measurements are not sized equally."

# Return the target variables
get_target(x::AbstractDirectProb{N,C}) where {N,C} = x.Y
get_target(x::AbstractDiscreteProb{N,C}) where {N,C} = x.X[:,2:end]
get_target(x::AbstracContProb{N,C}) where {N,C} = x.DX

get_oop_args(x::AbstractDataDrivenProblem{N,C,P}) where {N <: Number, C, P} = map(f->getfield(x, f), (:X, :p, :t, :U))

function get_oop_args(x::AbstractDiscreteProb{N,C}) where {N <: Number, C}
    return  (
        x.X[:, 1:end-1],
        x.p,
        x.t[1:end-1],
        x.U[:, 1:end-1]
    )
end

function get_implicit_oop_args(x::AbstractDirectProb{N,C}) where {N <: Number, C}
    return  (
        [x.X; x.Y],
        x.p,
        x.t,
        x.U
    )
end

function get_implicit_oop_args(x::AbstracContProb{N,C}) where {N <: Number, C}
    return  (
        [x.X; x.DX],
        x.p,
        x.t,
        x.U
    )
end

function get_implicit_oop_args(x::AbstractDiscreteProb{N,C}) where {N <: Number, C}
    return  (
        [x.X[:, 1:end-1]; x.X[:, 2:end]],
        x.p,
        x.t[1:end-1],
        x.U[:, 1:end-1]
    )
end

"""
$(SIGNATURES)

Checks if a `DataDrivenProblem` is valid by checking if the data contains `NaN`, `Inf` and
if the number of measurements is consistent.

# Example

```julia
is_valid(problem)
```
"""
function is_valid(x::AbstractDirectProb{N,C}) where {N <: Number,C}
    map(check_domain, (x.X, x.Y, x.U, x.t, x.p))
    if !C && !isempty(x.t)
        check_lengths(x.X, x.Y, x.U, x.t)
    elseif !C
        check_lengths(x.X, x.Y, x.U)
    elseif C && !isempty(x.t)
        check_lengths(x.X, x.Y, x.t)
    else
        check_lengths(x.X, x.Y)
    end
    return true
end

function is_valid(x::AbstractDiscreteProb{N,C}) where {N <: Number,C}
    map(check_domain, (x.X, x.U, x.t, x.p))
    if !C && !isempty(x.t)
        check_lengths(x.X, x.U, x.t)
    elseif !C
        check_lengths(x.X, x.U)
    elseif C && !isempty(x.t)
        check_lengths(x.X, x.t)
    else
        check_lengths(x.X)
    end
    return true
end

function is_valid(x::AbstracContProb{N,C}) where {N <: Number,C}
    map(check_domain, (x.X, x.DX, x.U, x.t, x.p))
    if !C && !isempty(x.t)
        check_lengths(x.X, x.DX, x.U, x.t)
    elseif !C
        check_lengths(x.X, x.DX, x.U)
    elseif C && !isempty(x.t)
        check_lengths(x.X, x.DX, x.t)
    else
        check_lengths(x.X, x.DX)
    end
    return true
end

"""
$(SIGNATURES)

Asserts if the given combination of `problem` and `basis` - and optionally a target matrix `dx` for in place evaluation - is
valid in its dimensions. Should only be used for checking explicit evaluation. Can also be used as a shorthand to

```julia
@is_applicable problem # Checks if the problem is well posed in terms of dimensions
@is_applicable problem basis # Checks if the basis can be called with the problem out of place
@is_applicable problem basis dx # Checks if the basis can be called with the problem and matrix dx in place
```
"""
macro is_applicable(problem)
    return :(@assert is_valid($(esc(problem))))
end

macro is_applicable(problem, basis)
    return quote
        if isa($(esc(problem)), AbstractDirectProb)
            @assert length(states($(esc(basis)))) == size(observed($(esc(problem))),1) "Problem and basis need to have same observed size"
        else
            @assert length(states($(esc(basis)))) == size(states($(esc(problem))),1) "Problem and basis need to have same state size"
        end
        @assert length(controls($(esc(basis)))) == size(controls($(esc(problem))),1) "Problem and basis need to have same control size"
        @assert length(parameters($(esc(basis)))) <= length(parameters($(esc(problem)))) "Problem and basis need to have consistent parameter size"
    end
end

macro is_applicable(problem, basis, dx)
    return quote
        @is_applicable $(esc(problem)) $(esc(basis))
        lp = length($(esc(problem)))
        n, m = size($(esc(dx)))
        lb = length($(esc(basis)))
        if isa($(esc(problem)), AbstractDiscreteProb)
            @assert n == lb && m == lp-1 "Target array has to be of size ($lb, $lp)"
        else
            @assert n == lb && m == lp "Target array has to be of size ($lb, $lp)"
        end
    end
end

## DESolution dispatch

ContinuousDataDrivenProblem(sol::T; kwargs...) where T <: DiffEqBase.DESolution = DataDrivenProblem(sol; kwargs...)
DiscreteDataDrivenProblem(sol::T; kwargs...) where T <: DiffEqBase.DESolution = DataDrivenProblem(sol; kwargs...)

function DataDrivenProblem(sol::T; use_interpolation = false, kwargs...) where T <: DiffEqBase.DESolution
    if sol.retcode != :Success
        throw(AssertionError("The solution is not successful. Abort."))
        return
    end

    X = Array(sol)
    t = sol.t
    p = sol.prob.p

    p = isa(p, DiffEqBase.NullParameters) ? eltype(X)[] : p

    if isdiscrete(sol.alg)

        return DiscreteDataDrivenProblem(
            X, t; p = p, kwargs...
        )

    else
        if use_interpolation
            DX = Array(sol(sol.t, Val{1}))
        else
            DX = similar(X)
            if DiffEqBase.isinplace(sol.prob.f)
                @views map(i->sol.prob.f(DX[:, i], X[:, i], p, t[i]), 1:size(X,2))
            else
                map(i->DX[:, i] .= sol.prob.f(X[:, i], p, t[i]), 1:size(X,2))
            end
        end

        return ContinuousDataDrivenProblem(
            X, t; DX = DX, p = p, kwargs...
        )
    end

end


text:
# Autogenerated wrapper script for alsa_plugins_jll for armv7l-linux-musleabihf
export libasound_module_conf_pulse, libasound_module_ctl_arcam_av, libasound_module_ctl_oss, libasound_module_ctl_pulse, libasound_module_pcm_a52, libasound_module_pcm_oss, libasound_module_pcm_pulse, libasound_module_pcm_speex, libasound_module_pcm_upmix, libasound_module_pcm_usb_stream, libasound_module_pcm_vdownmix, libasound_module_rate_lavrate, libasound_module_rate_samplerate, libasound_module_rate_speexrate

using FFMPEG_jll
using alsa_jll
using libsamplerate_jll
using PulseAudio_jll
JLLWrappers.@generate_wrapper_header("alsa_plugins")
JLLWrappers.@declare_library_product(libasound_module_conf_pulse, "libasound_module_conf_pulse.so")
JLLWrappers.@declare_library_product(libasound_module_ctl_arcam_av, "libasound_module_ctl_arcam_av.so")
JLLWrappers.@declare_library_product(libasound_module_ctl_oss, "libasound_module_ctl_oss.so")
JLLWrappers.@declare_library_product(libasound_module_ctl_pulse, "libasound_module_ctl_pulse.so")
JLLWrappers.@declare_library_product(libasound_module_pcm_a52, "libasound_module_pcm_a52.so")
JLLWrappers.@declare_library_product(libasound_module_pcm_oss, "libasound_module_pcm_oss.so")
JLLWrappers.@declare_library_product(libasound_module_pcm_pulse, "libasound_module_pcm_pulse.so")
JLLWrappers.@declare_library_product(libasound_module_pcm_speex, "libasound_module_pcm_speex.so")
JLLWrappers.@declare_library_product(libasound_module_pcm_upmix, "libasound_module_pcm_upmix.so")
JLLWrappers.@declare_library_product(libasound_module_pcm_usb_stream, "libasound_module_pcm_usb_stream.so")
JLLWrappers.@declare_library_product(libasound_module_pcm_vdownmix, "libasound_module_pcm_vdownmix.so")
JLLWrappers.@declare_library_product(libasound_module_rate_lavrate, "libasound_module_rate_lavrate.so")
JLLWrappers.@declare_library_product(libasound_module_rate_samplerate, "libasound_module_rate_samplerate.so")
JLLWrappers.@declare_library_product(libasound_module_rate_speexrate, "libasound_module_rate_speexrate.so")
function __init__()
    JLLWrappers.@generate_init_header(FFMPEG_jll, alsa_jll, libsamplerate_jll, PulseAudio_jll)
    JLLWrappers.@init_library_product(
        libasound_module_conf_pulse,
        "lib/alsa-lib/libasound_module_conf_pulse.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_ctl_arcam_av,
        "lib/alsa-lib/libasound_module_ctl_arcam_av.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_ctl_oss,
        "lib/alsa-lib/libasound_module_ctl_oss.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_ctl_pulse,
        "lib/alsa-lib/libasound_module_ctl_pulse.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_pcm_a52,
        "lib/alsa-lib/libasound_module_pcm_a52.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_pcm_oss,
        "lib/alsa-lib/libasound_module_pcm_oss.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_pcm_pulse,
        "lib/alsa-lib/libasound_module_pcm_pulse.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_pcm_speex,
        "lib/alsa-lib/libasound_module_pcm_speex.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_pcm_upmix,
        "lib/alsa-lib/libasound_module_pcm_upmix.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_pcm_usb_stream,
        "lib/alsa-lib/libasound_module_pcm_usb_stream.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_pcm_vdownmix,
        "lib/alsa-lib/libasound_module_pcm_vdownmix.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_rate_lavrate,
        "lib/alsa-lib/libasound_module_rate_lavrate.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_rate_samplerate,
        "lib/alsa-lib/libasound_module_rate_samplerate.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libasound_module_rate_speexrate,
        "lib/alsa-lib/libasound_module_rate_speexrate.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#type CompositeMap <: Map
#  header::MapHeader
#  a::Map
#  b::Map
#  function CompositeMap(a::Map, b::Map)
#    codomain(a) == domain(b) || throw("maps not compatible")
#    r = new()
#    r.a = a
#    r.b = b
#    image = function(M::Map, a::Any)
#      return M.a(M.b(a))
#    end
#    preim = function(M::Map, a::Any)
#      return preimage(M.b, preimage(M.a, a))
#    end
#    r.header = MapHeader(a.header.domain, b.header.codomain, image, preim)
#    return r
#  end
#end
#
#function show(io::IO, M::CompositeMap)
#  println(io, "composite of")
#  println(io, " ", M.a)
#  println(io, "*")
#  println(io, " ", M.b)
#end
#
#function *(a::Map, b::Map)
#  return CompositeMap(a,b)
#end
#


text:
using PrimitiveOneHot
using PrimitiveOneHot: OneHot, onehotsize

# old utility

onehot2indices(xs::OneHotArray) = reinterpret(Int32, xs.onehots)

indices2onehot(nums::Int, xs::AbstractArray{Int}) = OneHotArray(nums, xs)

tofloat(::Type{F}, o::OneHotArray{N, <:AbstractArray}) where {F<:AbstractFloat,N} = Array{F}(o)

# AD
using Flux
Flux.@nograd onehot2indices
Flux.@nograd indices2onehot
Flux.@nograd tofloat


text:
abstract type AbstractMatrixOperator{Tconv} end

@params struct MatrixOperator{Tconv} <: AbstractMatrixOperator{Tconv}
    K::Any
    f::Any
    conv::Tconv
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::MatrixOperator)
    return println("TopOpt matrix linear operator")
end
LinearAlgebra.mul!(c, op::MatrixOperator, b) = mul!(c, op.K, b)
Base.size(op::MatrixOperator, i...) = size(op.K, i...)
Base.eltype(op::MatrixOperator) = eltype(op.K)
LinearAlgebra.:*(op::MatrixOperator, b) = mul!(similar(b), op.K, b)

@params struct MatrixFreeOperator{Tconv,T,dim} <: AbstractMatrixOperator{Tconv}
    f::AbstractVector{T}
    elementinfo::ElementFEAInfo{dim,T}
    meandiag::T
    vars::AbstractVector{T}
    xes::Any
    fixed_dofs::Any
    free_dofs::Any
    xmin::T
    penalty::Any
    conv::Tconv
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::MatrixFreeOperator)
    return println("TopOpt matrix-free linear operator")
end
Base.size(op::MatrixFreeOperator) = (size(op, 1), size(op, 2))
Base.size(op::MatrixFreeOperator, i) = 1 <= i <= 2 ? length(op.elementinfo.fixedload) : 1
Base.eltype(op::MatrixFreeOperator{<:Any,T}) where {T} = T

import LinearAlgebra: *, mul!

function *(A::MatrixFreeOperator, x)
    y = similar(x)
    mul!(y, A::MatrixFreeOperator, x)
    return y
end

function mul!(y::TV, A::MatrixFreeOperator, x::TV) where {TV<:AbstractVector}
    T = eltype(y)
    nels = length(A.elementinfo.Kes)
    ndofs = length(A.elementinfo.fixedload)
    dofspercell = size(A.elementinfo.Kes[1], 1)
    meandiag = A.meandiag

    @unpack Kes, metadata, black, white, varind = A.elementinfo
    @unpack cell_dofs, dof_cells = metadata
    @unpack penalty, xmin, vars, fixed_dofs, free_dofs, xes = A

    for i in 1:nels
        if PENALTY_BEFORE_INTERPOLATION
            px = ifelse(
                black[i],
                one(T),
                ifelse(white[i], xmin, density(penalty(vars[varind[i]]), xmin)),
            )
        else
            px = ifelse(
                black[i],
                one(T),
                ifelse(white[i], xmin, penalty(density(vars[varind[i]], xmin))),
            )
        end
        xe = xes[i]
        for j in 1:dofspercell
            xe = @set xe[j] = x[cell_dofs[j, i]]
        end
        if eltype(Kes) <: Symmetric
            xes[i] = px * (bcmatrix(Kes[i]).data * xe)
        else
            xes[i] = px * (bcmatrix(Kes[i]) * xe)
        end
    end

    for i in 1:length(fixed_dofs)
        dof = fixed_dofs[i]
        y[dof] = meandiag * x[dof]
    end
    for i in 1:length(free_dofs)
        dof = free_dofs[i]
        yi = zero(T)
        r = dof_cells.offsets[dof]:(dof_cells.offsets[dof + 1] - 1)
        for ind in r
            k, m = dof_cells.values[ind]
            yi += xes[k][m]
        end
        y[dof] = yi
    end
    return y
end


text:
using ModelingToolkit, Unitful, OrdinaryDiffEq, DiffEqJump, IfElse
using Test
MT = ModelingToolkit
@parameters τ [unit = u"ms"] γ
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
D = Differential(t)

#This is how equivalent works:
@test MT.equivalent(u"MW" ,u"kJ/ms")
@test !MT.equivalent(u"m", u"cm")
@test MT.equivalent(MT.get_unit(P^γ), MT.get_unit((E/τ)^γ))

# Basic access
@test MT.get_unit(t) == u"ms"
@test MT.get_unit(E) == u"kJ"
@test MT.get_unit(τ) == u"ms"
@test MT.get_unit(γ) == MT.unitless
@test MT.get_unit(0.5) == MT.unitless
@test MT.get_unit(MT.SciMLBase.NullParameters) == MT.unitless

# Prohibited unit types
@parameters β [unit = u"°"] α [unit = u"°C"] γ [unit = 1u"s"]
@test_throws MT.ValidationError MT.get_unit(β)
@test_throws MT.ValidationError MT.get_unit(α)
@test_throws MT.ValidationError MT.get_unit(γ)

# Non-trivial equivalence & operators
@test MT.get_unit(τ^-1) == u"ms^-1"
@test MT.equivalent(MT.get_unit(D(E)),u"MW")
@test MT.equivalent(MT.get_unit(E/τ), u"MW")
@test MT.get_unit(2*P) == u"MW"
@test MT.get_unit(t/τ) == MT.unitless
@test MT.equivalent(MT.get_unit(P - E/τ),u"MW")
@test MT.equivalent(MT.get_unit(D(D(E))),u"MW/ms")
@test MT.get_unit(IfElse.ifelse(t>t,P,E/τ)) == u"MW"
@test MT.get_unit(1.0^(t/τ)) == MT.unitless
@test MT.get_unit(exp(t/τ)) == MT.unitless
@test MT.get_unit(sin(t/τ)) == MT.unitless
@test MT.get_unit(sin(1u"rad")) == MT.unitless
@test MT.get_unit(t^2) == u"ms^2"

eqs = [D(E) ~ P - E/τ
        0 ~ P]
@test MT.validate(eqs)
@named sys = ODESystem(eqs)

@test !MT.validate(D(D(E)) ~ P)
@test !MT.validate(0 ~ P + E*τ)

# Array variables
@variables t [unit = u"s"] x[1:3](t) [unit = u"m"]
@parameters v[1:3] = [1,2,3] [unit = u"m/s"]
D = Differential(t)
eqs = D.(x) .~ v
ODESystem(eqs,name=:sys)

# Difference equation
@parameters t [unit = u"s"] a [unit = u"s"^-1]
@variables x(t) [unit = u"kg"]
δ = Differential(t)
D = Difference(t; dt = 0.1u"s")
eqs = [
    δ(x) ~ a*x
]
de = ODESystem(eqs, t, [x], [a],name=:sys)

# Nonlinear system
@parameters a [unit = u"kg"^-1]
@variables x [unit = u"kg"]
eqs = [
    0 ~ a*x
]
@named nls = NonlinearSystem(eqs, [x], [a])

# SDE test w/ noise vector
@parameters τ [unit = u"ms"] Q [unit = u"MW"]
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
D = Differential(t)
eqs = [D(E) ~ P - E/τ
      P ~ Q]

noiseeqs = [0.1u"MW",
            0.1u"MW"]
@named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])

# With noise matrix
noiseeqs = [0.1u"MW" 0.1u"MW"
            0.1u"MW" 0.1u"MW"]
@named sys = SDESystem(eqs,noiseeqs, t, [P, E], [τ, Q])

# Invalid noise matrix
noiseeqs = [0.1u"MW" 0.1u"MW"
            0.1u"MW" 0.1u"s"]
@test !MT.validate(eqs,noiseeqs)

# Non-trivial simplifications
@variables t [unit = u"s"] V(t) [unit = u"m"^3] L(t) [unit = u"m"]
@parameters v [unit = u"m/s"] r [unit =u"m"^3/u"s"]
D = Differential(t)
eqs = [D(L) ~ v,
       V ~ L^3]
@named sys = ODESystem(eqs)
sys_simple = structural_simplify(sys)

eqs = [D(V) ~ r,
       V ~ L^3]
@named sys = ODESystem(eqs)
sys_simple = structural_simplify(sys)

@variables V [unit = u"m"^3] L [unit = u"m"]
@parameters v [unit = u"m/s"] r [unit = u"m"^3/u"s"] t [unit = u"s"]
eqs = [V ~ r*t,
       V ~ L^3]
@named sys = NonlinearSystem(eqs, [V, L], [t, r])
sys_simple = structural_simplify(sys)

eqs = [L ~ v*t,
       V ~ L^3]
@named sys = NonlinearSystem(eqs, [V,L], [t,r])
sys_simple = structural_simplify(sys)

#Jump System
@parameters β [unit = u"(mol^2*s)^-1"] γ [unit = u"(mol*s)^-1"] t [unit = u"s"] jumpmol [unit = u"mol"]
@variables S(t) [unit = u"mol"] I(t) [unit = u"mol"] R(t) [unit = u"mol"]
rate₁   = β*S*I
affect₁ = [S ~ S - 1*jumpmol, I ~ I + 1*jumpmol]
rate₂   = γ*I
affect₂ = [I ~ I - 1*jumpmol, R ~ R + 1*jumpmol]
j₁      = ConstantRateJump(rate₁, affect₁)
j₂      = VariableRateJump(rate₂, affect₂)
js      = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ],name=:sys)

affect_wrong = [S ~ S - jumpmol, I ~ I + 1]
j_wrong      = ConstantRateJump(rate₁, affect_wrong)
@test_throws MT.ValidationError JumpSystem([j_wrong, j₂], t, [S, I, R], [β, γ],name=:sys)

rate_wrong   = γ^2*I
j_wrong     = ConstantRateJump(rate_wrong, affect₂)
@test_throws MT.ValidationError JumpSystem([j₁, j_wrong], t, [S, I, R], [β, γ],name=:sys)

# mass action jump tests for SIR model
maj1 = MassActionJump(2*β/2, [S => 1, I => 1], [S => -1, I => 1])
maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
@named js3  = JumpSystem([maj1, maj2], t, [S,I,R], [β,γ])

#Test unusual jump system
@parameters β γ t
@variables S(t) I(t) R(t)

maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
maj2 = MassActionJump(γ, [S => 1], [S => -1])
@named js4  = JumpSystem([maj1, maj2], t, [S], [β, γ])


text:
"""
an example of optim_algorithm is NLopt.LN_BOBYQA
"""
function setupβLSsolverMultistartoptim(optim_algorithm,
    Es::Vector{NMRSpectraSimulator.κCompoundFIDType{T, SST}},
    As::Vector{NMRSpectraSimulator.CompoundFIDType{T, SST}},
    LS_inds,
    U_rad_cost,
    y_cost::Vector{Complex{T}},
    κs_β_DOFs, κs_β_orderings;
    κ_lb_default = 0.2,
    κ_ub_default = 5.0,
    max_iters = 50,
    xtol_rel = 1e-9,
    ftol_rel = 1e-9,
    maxtime = Inf,
    N_starts = 100,
    inner_xtol_rel = 1e-12,
    inner_maxeval = 100,
    inner_maxtime = Inf) where {T,SST}

    #
    N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
    β_lb = ones(T, N_β) .* (-π)
    β_ub = ones(T, N_β) .* (π)
    #β_initial = zeros(T, N_β)

    p_lb = β_lb
    p_ub = β_ub

    q, updateβfunc, updateκfunc, E_BLS, κ_BLS, b_BLS,
    getβfunc = setupcostβLS(Es, Bs, As, LS_inds, U_rad_cost, y_cost,
        κs_β_DOFs,
        κs_β_orderings;
        κ_lb_default = κ_lb_default,
        κ_ub_default = κ_ub_default)

    f = pp->costβLS(U_rad_cost, y_cost, updateβfunc, updateκfunc, pp, Es,
    E_BLS, κ_BLS, b_BLS, q)

    P = MultistartOptimization.MinimizationProblem(f,
        p_lb, p_ub)
    #
    local_method = MultistartOptimization.NLoptLocalMethod(; algorithm = optim_algorithm,
    xtol_rel = inner_xtol_rel,
    maxeval = inner_maxeval,
    maxtime = inner_maxtime)

    multistart_method = MultistartOptimization.TikTak(N_starts)
    run_optim = pp->MultistartOptimization.multistart_minimization(multistart_method,
        local_method, P)

    return run_optim, f, E_BLS, κ_BLS, b_BLS, updateβfunc, q
end

function setupβLSsolver(optim_algorithm,
    Es,
    Bs,
    As,
    LS_inds,
    U_rad_cost,
    y_cost::Vector{Complex{T}},
    κs_β_DOFs::Vector{Vector{Int}},
    κs_β_orderings::Vector{Vector{Vector{Int}}};
    κ_lb_default = 0.2,
    κ_ub_default = 5.0,
    max_iters = 50,
    xtol_rel = 1e-9,
    ftol_rel = 1e-9,
    maxtime = Inf) where T

    #N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
    N_β = sum( getNβ(Bs[n]) for n = 1:length(Bs) )

    β_lb = ones(T, N_β) .* (-π)
    β_ub = ones(T, N_β) .* (π)
    #β_initial = zeros(T, N_β)

    p_lb = β_lb
    p_ub = β_ub

    q, updateβfunc, updateκfunc, E_BLS, κ_BLS, b_BLS,
    getβfunc = setupcostβLS(Es, Bs, As, LS_inds, U_rad_cost, y_cost,
    κs_β_DOFs, κs_β_orderings;
        κ_lb_default = κ_lb_default,
        κ_ub_default = κ_ub_default)

    # q is simulated spectra. f is cost function.
    f = pp->costβLS(U_rad_cost, y_cost, updateβfunc, updateκfunc, pp, Es,
    E_BLS, κ_BLS, b_BLS, q)

    df = xx->FiniteDiff.finite_difference_gradient(f, xx)

    opt = NLopt.Opt(optim_algorithm, N_β)

    run_optim = pp->runNLopt!(opt,
        pp,
        f,
        df,
        p_lb,
        p_ub;
        max_iters = max_iters,
        xtol_rel = xtol_rel,
        ftol_rel = ftol_rel,
        maxtime = maxtime)

    return run_optim, f, E_BLS, κ_BLS, b_BLS, updateβfunc, updateκfunc, q
end

"""
κ version.
"""
function setupcostβLS(Es::Vector{NMRSpectraSimulator.καFIDModelType{T, SST}},
    Bs,
    As,
    LS_inds,
    U0_rad,
    y0::Vector{Complex{T}},
    κs_β_DOFs,
    κs_β_orderings;
    κ_lb_default = 0.2,
    κ_ub_default = 5.0) where {T <: Real, SST}

    w = ones(T, length(Es))

    #N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
    N_β = sum( getNβ(Bs[n]) for n = 1:length(Bs) )

    st_ind_β = 1
    fin_ind_β = st_ind_β + N_β - 1
    #updateβfunc = pp->updateβ!(Bs, κs_β_orderings, κs_β_DOFs, pp, st_ind_β)
    updateβfunc = pp->updateβ!(Bs, pp, st_ind_β)

    f = uu->NMRSpectraSimulator.evalitpproxymixture(uu, Bs, Es; w = w)

    ### LS κ.
    U_rad_LS = U0_rad[LS_inds]

    N_κ, N_κ_singlets = countκ(Es)
    N_κ_vars = N_κ + N_κ_singlets
    E_BLS, κ_BLS, b_BLS = setupupdateLS(length(U_rad_LS), N_κ_vars, y0[LS_inds])

    κ_lb = ones(N_κ_vars) .* κ_lb_default
    κ_ub = ones(N_κ_vars) .* κ_ub_default

    updateκfunc = xx->updateκ2!(E_BLS, b_BLS, κ_BLS,
    U_rad_LS, Es, As, w, κ_lb, κ_ub)

    #### extract parameters from p.
    getβfunc = pp->pp[st_ind_β:fin_ind_β]

    return f, updateβfunc, updateκfunc, E_BLS, κ_BLS, b_BLS, getβfunc
end

function costβLS(U_rad,
    S_U::Vector{Complex{T}},
    updateβfunc, updateκfunc,
    p::Vector{T}, Es,
    E_BLS::Matrix{Complex{T}}, κ_BLS::Vector{T}, b_BLS,
    f)::T where T <: Real

    updateβfunc(p)

    # try updateκfunc(p)
    # catch e
    #     println("error in updateκfunc(p)")
    # 	bt = catch_backtrace()
    # 	showerror(stdout, e, bt)
    # 	rethrow(e)
    # end

    updateκfunc(p)
    parseκ!(Es, κ_BLS)
    #println("all(isfinite.(Es)) = ", all(isfinite.(Es)))
    # ## l-2 costfunc.
    # cost = zero(T)
    # for m = 1:length(S_U)
    #
    #     f_u = f(U_rad[m])
    #
    #     cost += abs2( f_u - S_U[m] )
    # end
    #
    # # faster l-2 cost compute.
    # B = reinterpret(T, E_BLS)
    # tmp = B*κ_BLS - b_BLS
    # cost = dot(tmp, tmp)
    cost = norm(reinterpret(T, E_BLS)*κ_BLS - b_BLS)^2 # compact version.

    return cost
end


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

mutable struct ApplicationGatewayBackendAddress2 <: SwaggerModel
    fqdn::Any # spec type: Union{ Nothing, String } # spec name: fqdn
    ipAddress::Any # spec type: Union{ Nothing, String } # spec name: ipAddress

    function ApplicationGatewayBackendAddress2(;fqdn=nothing, ipAddress=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendAddress2, Symbol("fqdn"), fqdn)
        setfield!(o, Symbol("fqdn"), fqdn)
        validate_property(ApplicationGatewayBackendAddress2, Symbol("ipAddress"), ipAddress)
        setfield!(o, Symbol("ipAddress"), ipAddress)
        o
    end
end # type ApplicationGatewayBackendAddress2

const _property_map_ApplicationGatewayBackendAddress2 = Dict{Symbol,Symbol}(Symbol("fqdn")=>Symbol("fqdn"), Symbol("ipAddress")=>Symbol("ipAddress"))
const _property_types_ApplicationGatewayBackendAddress2 = Dict{Symbol,String}(Symbol("fqdn")=>"String", Symbol("ipAddress")=>"String")
Base.propertynames(::Type{ ApplicationGatewayBackendAddress2 }) = collect(keys(_property_map_ApplicationGatewayBackendAddress2))
Swagger.property_type(::Type{ ApplicationGatewayBackendAddress2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendAddress2[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendAddress2 }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendAddress2[property_name]

function check_required(o::ApplicationGatewayBackendAddress2)
    true
end

function validate_property(::Type{ ApplicationGatewayBackendAddress2 }, name::Symbol, val)
end


text:
export EinCode, EinIndexer, EinArray
export einarray

"""
    EinCode
    EinCode(ixs, iy)

Abstract type for sum-product contraction code.
The constructor returns a `DynamicEinCode` instance.
"""
abstract type EinCode end

"""
    StaticEinCode{ixs, iy}

The static version of `DynamicEinCode` that matches the contraction rule at compile time.
It is the default return type of `@ein_str` macro.
"""
struct StaticEinCode{ixs, iy} <: EinCode end

getixs(::StaticEinCode{ixs}) where ixs = ixs
getiy(::StaticEinCode{ixs, iy}) where {ixs, iy} = iy
labeltype(::StaticEinCode{ixs,iy}) where {ixs, iy} = promote_type(eltype.(ixs)..., eltype(iy))
getixsv(code::StaticEinCode) = [collect(labeltype(code), ix) for ix in getixs(code)]
getiyv(code::StaticEinCode) = collect(labeltype(code), getiy(code))

"""
    DynamicEinCode{LT}
    DynamicEinCode(ixs, iy)

Wrapper to `eincode`-specification that creates a callable object
to evaluate the `eincode` `ixs -> iy` where `ixs` are the index-labels
of the input-tensors and `iy` are the index-labels of the output.

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b = rand(2,2), rand(2,2);

julia> OMEinsum.DynamicEinCode((('i','j'),('j','k')),('i','k'))(a, b) ≈ a * b
true
```
"""
struct DynamicEinCode{LT} <: EinCode
    ixs::Vector{Vector{LT}}
    iy::Vector{LT}
end
# to avoid ambiguity error, support tuple inputs
function DynamicEinCode(ixs, iy)
    @debug "generating dynamic eincode ..."
    if isempty(ixs)
        error("number of input tensors must be greater than 0")
    end
    DynamicEinCode(_tovec(ixs, iy)...)
end
_tovec(ixs::NTuple{N,Tuple{}}, iy::Tuple{}) where {N} = [collect(Union{}, ix) for ix in ixs], collect(Union{}, iy)
_tovec(ixs::NTuple{N,NTuple{M,LT} where M}, iy::NTuple{K,LT}) where {N,K,LT} = [collect(LT, ix) for ix in ixs], collect(LT, iy)
_tovec(ixs::NTuple{N,Vector{LT}}, iy::Vector{LT}) where {N,LT} = collect(ixs), iy
_tovec(ixs::AbstractVector{Vector{LT}}, iy::AbstractVector{LT}) where {N,K,LT} = collect(ixs), collect(iy)

Base.:(==)(x::DynamicEinCode, y::DynamicEinCode) = x.ixs == y.ixs && x.iy == y.iy
# forward from EinCode, for compatibility
EinCode(ixs, iy) = DynamicEinCode(ixs, iy)

getixs(code::DynamicEinCode) = code.ixs
getiy(code::DynamicEinCode) = code.iy
labeltype(::DynamicEinCode{LT}) where LT = LT
getixsv(code::DynamicEinCode) = code.ixs
getiyv(code::DynamicEinCode) = code.iy

# conversion
DynamicEinCode(::StaticEinCode{ixs, iy}) where {ixs, iy} = DynamicEinCode(ixs, iy)
StaticEinCode(code::DynamicEinCode) = StaticEinCode{(Tuple.(code.ixs)...,), (code.iy...,)}()

"""
    EinIndexer{locs,N}

A structure for indexing `EinArray`s. `locs` is the index positions (among all indices).
In the constructor, `size` is the size of target tensor,
"""
struct EinIndexer{locs,N}
    cumsize::NTuple{N, Int}
end

"""
    EinIndexer{locs}(size::Tuple)

Constructor for `EinIndexer` for an object of size `size` where `locs` are the
locations of relevant indices in a larger tuple.
"""
function EinIndexer{locs}(size::NTuple{N,Int}) where {N,locs}
    N==0 && return EinIndexer{(),0}(())
    EinIndexer{locs,N}((1,TupleTools.cumprod(size[1:end-1])...))
end

getlocs(::EinIndexer{locs,N}) where {N,locs} = locs

# get a subset of index
@inline subindex(indexer::EinIndexer, ind::CartesianIndex) = subindex(indexer, ind.I)
@inline subindex(indexer::EinIndexer{(),0}, ind::NTuple{N0,Int}) where N0 = 1
@inline @generated function subindex(indexer::EinIndexer{locs,N}, ind::NTuple{N0,Int}) where {N,N0,locs}
    ex = Expr(:call, :+, map(i->i==1 ? :(ind[$(locs[i])]) : :((ind[$(locs[i])]-1) * indexer.cumsize[$i]), 1:N)...)
    :(@inbounds $ex)
end

"""
    EinArray{T, N, TT, LX, LY, ICT, OCT} <: AbstractArray{T, N}

A struct to hold the intermediate result of an `einsum` where all index-labels
of both input and output are expanded to a rank-`N`-array whose values
are lazily calculated.
Indices are arranged as _inner indices_ (or reduced dimensions) first and _then outer indices_.

Type parameters are

    * `T`: element type,
    * `N`: array dimension,
    * `TT`: type of "tuple of input arrays",
    * `LX`: type of "tuple of input indexers",
    * `LX`: type of output indexer,
    * `ICT`: typeof inner CartesianIndices,
    * `OCT`: typeof outer CartesianIndices,
"""
struct EinArray{T, N, TT, LX, LY, ICT, OCT} <: AbstractArray{T, N}
    xs::TT
    x_indexers::LX
    y_indexer::LY
    size::NTuple{N, Int}
    ICIS::ICT
    OCIS::OCT
    function EinArray{T}(xs::TT, x_indexers::LX, y_indexer::LY, size::NTuple{N, Int}, ICIS::ICT, OCIS::OCT) where {T,N,TT<:Tuple,LX<:Tuple,LY<:EinIndexer,ICT,OCT}
        new{T,N,TT,LX,LY,ICT,OCT}(xs,x_indexers,y_indexer,size,ICIS,OCIS)
    end
end

"""
    einarray(::Val{ixs}, Val{iy}, xs, size_dict) -> EinArray

Constructor of `EinArray` from an `EinCode`, a tuple of tensors `xs` and a `size_dict` that assigns each index-label a size.
The returned `EinArray` holds an intermediate result of the `einsum` specified by the
`EinCode` with indices corresponding to all unique labels in the einsum.
Reduction over the (lazily calculated) dimensions that correspond to labels not present
in the output lead to the result of the einsum.

# example

```jldoctest; setup = :(using OMEinsum)
julia> using OMEinsum: get_size_dict

julia> a, b = rand(2,2), rand(2,2);

julia> sd = get_size_dict((('i','j'),('j','k')), (a, b));

julia> ea = OMEinsum.einarray(Val((('i','j'),('j','k'))),Val(('i','k')), (a,b), sd);

julia> dropdims(sum(ea, dims=1), dims=1) ≈ a * b
true
```
"""
@generated function einarray(::Val{ixs}, ::Val{iy}, xs::TT, size_dict) where {ixs, iy, NI, TT<:NTuple{NI,AbstractArray}}
    inner_indices, outer_indices, locs_xs, locs_y = indices_and_locs(ixs, iy)
    inner_indices = (inner_indices...,)
    outer_indices = (outer_indices...,)
    T = promote_type(eltype.(xs.parameters)...)
    xind_expr = Expr(:call, :tuple, map(i->:(EinIndexer{$(locs_xs[i])}(size(xs[$i]))),1:NI)...)

    quote
        # find size for each leg
        outer_sizes = getindex.(Ref(size_dict), $outer_indices)
        inner_sizes = getindex.(Ref(size_dict), $inner_indices)

        # cartesian indices for outer and inner legs
        outer_ci = CartesianIndices(outer_sizes)
        inner_ci = CartesianIndices(inner_sizes)

        x_indexers = $(xind_expr)
        y_size = getindex.(Ref(size_dict), iy)
        y_indexer = EinIndexer{$locs_y}(y_size)

        #EinArray{$T, $N, $TT, $LX, $LY, $ICT, $OCT}(xs, x_indexers, y_indexer, (inner_sizes...,outer_sizes...), inner_ci, outer_ci)
        EinArray{$T}(xs, x_indexers, y_indexer, (inner_sizes...,outer_sizes...), inner_ci, outer_ci)
    end
end

Base.size(A::EinArray) = A.size
@doc raw"
    getindex(A::EinArray, inds...)
return the lazily calculated entry of `A` at index `inds`.
"
@inline Base.getindex(A::EinArray{T}, ind) where {T} = map_prod(A.xs, ind, A.x_indexers)
@inline Base.getindex(A::EinArray{T}, inds::Int...) where {T} = map_prod(A.xs, inds, A.x_indexers)

# get one element from each tensor, and multiple them
@doc raw"
    map_prod(xs, ind, indexers)

calculate the value of an `EinArray` with `EinIndexer`s `indexers` at location `ind`.
"
@inline @generated function map_prod(xs::Tuple, ind, indexers::NTuple{N,Any}) where {N}
    ex = Expr(:call, :*, map(i->:(xs[$i][subindex(indexers[$i], ind)]), 1:N)...)
    :(@inbounds $ex)
end

@doc raw"
    indices_and_locs(ixs,iy)
given the index-labels of input and output of an `einsum`, return
(in the same order):
- a tuple of the distinct index-labels of the output `iy`
- a tuple of the distinct index-labels in `ixs` of the input not appearing in the output `iy`
- a tuple of tuples of locations of an index-label in the `ixs` in a list of all index-labels
- a tuple of locations of index-labels in `iy` in a list of all index-labels

where the list of all index-labels is simply the first  and the second output catenated and the second output catenated.
"
function indices_and_locs(ixs, iy)
    # outer legs and inner legs
    outer_indices = unique!(collect(iy))
    inner_indices = setdiff!(collect(vcat(_collect.(ixs)...)), outer_indices)

    # for indexing tensors (leg binding)
    indices = (inner_indices...,outer_indices...)
    locs_xs = map(ixs) do ix
        map(i->findfirst(isequal(i), indices), ix)
    end
    locs_y = map(i->findfirst(isequal(i), outer_indices), iy)
    return inner_indices, outer_indices, locs_xs, locs_y
end

# dynamic EinIndexer
struct DynamicEinIndexer{N}
    locs::NTuple{N,Int}
    cumsize::NTuple{N, Int}
end

function dynamic_indexer(locs::NTuple{N,Int}, size::NTuple{N,Int}) where {N}
    N==0 && return DynamicEinIndexer{0}((),())
    DynamicEinIndexer{N}(locs, (1,TupleTools.cumprod(size[1:end-1])...))
end

getlocs(d::DynamicEinIndexer{N}) where {N} = d.locs

# get a subset of index
@inline subindex(indexer::DynamicEinIndexer, ind::CartesianIndex) = subindex(indexer, ind.I)
@inline subindex(indexer::DynamicEinIndexer{0}, ind::NTuple{N0,Int}) where N0 = 1
@inline @generated function subindex(indexer::DynamicEinIndexer{N}, ind::NTuple{N0,Int}) where {N,N0}
    ex = Expr(:call, :+, map(i->i==1 ? :(ind[indexer.locs[$i]]) : :((ind[indexer.locs[$i]]-1) * indexer.cumsize[$i]), 1:N)...)
    :(@inbounds $ex)
end


text:
"""
This set of functions queries if a given coordinate/grid is within a region defined either
by a set of region [N,S,E,W] bounds or by longitude/latitude vectors.

Other additional features include:
    * Transformation of longitude coordinates from [-180,180] to [0,360] and vice versa

"""

## Is Point in GeoRegion?

"""
    isinGeoRegion(
        Point  :: Point2{<:Real},
        geo    :: GeoRegion;
        tlon   :: Real = 0,
        tlat   :: Real = 0,
        throw  :: Bool = true
    ) -> Bool

Check if a geographical point `Point` is within a GeoRegion defined by `geo`.

Arguments
=========

- `Point` : A geographical point of Type `Point2`.  Pass `Point2(plon,plat)`, where `plon` and `plat` are the longitude and latitudes of the point.
- `geo`   : The GeoRegion struct container

Keyword Arguments
=================

- `tlon`  : Threshold for the longitude bound
- `tlat`  : Threshold for the latitude bound
- `throw` : If `true`, then if `Point` is not within `geo`, an error is thrown and the program stops running.
"""
function isinGeoRegion(
    Point  :: Point2{<:Real},
    GeoReg :: RectRegion;
    tlon   :: Real = 0,
    tlat   :: Real = 0,
    throw  :: Bool = true
)

    if throw
        @info "$(modulelog()) - Performing a check to determine if the coordinates $(Point) are within the specified region boundaries."
    end

    N = GeoReg.N
    S = GeoReg.S
    E = GeoReg.E
    W = GeoReg.W

    plon = Point[1]; plon = mod(plon,360)
    plat = Point[2]

    is180 = GeoReg.is180
    is360 = GeoReg.is360

    isin = true

    if (plat < (S-tlat)) || (plat > (N+tlat))
        isin = false
    else
        if is360
            if !is180
                if (plon < (W-tlon)) || (plon > (E+tlon)); isin = false end
            else

                if (plon < (W-tlon)) || (plon > (E+tlon))
                    plon = plon - 360
                    if (plon < (W-tlon)) || (plon > (E+tlon))
                        isin = false
                    end
                end

            end
        else

            if plon > 180; plon = plon - 360 end
            if (plon < (W-tlon)) || (plon > (E+tlon))
                if (plon < (W-tlon)) || (plon > (E+tlon))
                    isin = false
                end
            end

        end
    end

    if !isin
        if throw
             error("$(modulelog()) - The requested coordinates $(Point2(plon,plat)) are not within the specified region boundaries.")
        else return false
        end
    else
        if throw
            @info "$(modulelog()) - The requested coordinates $(Point2(plon,plat)) are within the specified region boundaries."
        end
        return true
    end

end

function isinGeoRegion(
    Point  :: Point2{<:Real},
    GeoReg :: PolyRegion;
    tlon   :: Real = 0,
    tlat   :: Real = 0,
    throw  :: Bool = true
)

    if throw
        @info "$(modulelog()) - Performing a check to determine if the coordinates $(Point) are within the specified region boundaries."
    end

    N = GeoReg.N
    S = GeoReg.S
    E = GeoReg.E
    W = GeoReg.W
    GeoShape = GeoReg.shape

    plon = Point[1]; plon = mod(plon,360)
    plat = Point[2]

    is180 = GeoReg.is180
    is360 = GeoReg.is360

    isin = true

    if (plat < (S-tlat)) || (plat > (N+tlat))
        isin = false
    else
        if is360
            if !is180
                if (plon < (W-tlon)) || (plon > (E+tlon)); isin = false end
            else

                if (plon < (W-tlon)) || (plon > (E+tlon))
                    plon = plon - 360
                    if (plon < (W-tlon)) || (plon > (E+tlon))
                        isin = false
                    end
                end

            end
        else

            if plon > 180; plon = plon - 360 end
            if (plon < (W-tlon)) || (plon > (E+tlon))
                if (plon < (W-tlon)) || (plon > (E+tlon))
                    isin = false
                end
            end

        end
    end

    if isin
        isin = inpolygon(Point2(plon,plat),GeoReg.shape,in=true,on=true,out=false)
    end

    if !isin
        if throw
            error("$(modulelog()) - The requested coordinates $(Point2(plon,plat)) are not within the specified region boundaries.")
        else
            return false
        end
    else
        if throw
            @info "$(modulelog()) - The requested coordinates $(Point2(plon,plat)) are within the specified region boundaries."
        end
        return true
    end

end

"""
    isinGeoRegion(
        Child  :: GeoRegion,
        polyG  :: PolyRegion;
        domask :: Bool = false,
        throw  :: Bool = true
    ) -> Bool

Check if a child GeoRegion defined by `Child` is within a PolyRegion `polyG`.

Arguments
=========

- `Child` : A GeoRegion that we postulate to be a "child", or a subset of the GeoRegion defined by `polyG`
- `polyG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`

Keyword Arguments
=================

- `throw`  : If `true`, then if `Child` is not within `polyG`, an error is thrown and the program stops running
- `domask` : If `throw` is `false` and `domask` is `true`, return a mask (with bounds defined by the `Child` GeoRegion) showing the region where `Child` and `polyG` do not overlap
"""
function isinGeoRegion(
    Child  :: GeoRegion,
    polyG  :: PolyRegion;
    domask :: Bool = false,
    throw  :: Bool = true
)

    @info "$(modulelog()) - Performing a check to determine if the $(Child.name) GeoRegion ($(Child.regID)) is inside the $(polyG.name) GeoRegion ($(polyG.regID))"

    N = Child.N
    S = Child.S
    E = Child.E
    W = Child.W

    lon = range(W,E,length=10001); nlon = length(lon)
    lat = range(S,N,length=10001); nlat = length(lat)

    mask = zeros(Bool,nlon,nlat)

    for ilat in 1 : nlat, ilon = 1 : nlon
        if isinGeoRegion(Point2(lon[ilon],lat[ilat]),Child,throw=false)
            if !isinGeoRegion(Point2(lon[ilon],lat[ilat]),polyG,throw=false)
                mask[ilon,ilat] = 1
            end
        end
    end

    if sum(mask) > 0

        if throw
            error("$(modulelog()) - The GeoRegion $(Child.regID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.regID) ($(polyG.name))")
        else
            if domask
                @warn "$(modulelog()) - The GeoRegion $(Child.regID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.regID) ($(polyG.name)), returning a mask to show which regions do not intersect"
                return mask
            else
                @warn "$(modulelog()) - The GeoRegion $(Child.regID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.regID) ($(polyG.name))"
                return false
            end
        end

    else
        @info "$(modulelog()) - The GeoRegion $(Child.regID) ($(Child.name)) is indeed a subset of the GeoRegion $(polyG.regID) ($(polyG.name))"
        return true
    end

end

"""
    isinGeoRegion(
        Child  :: GeoRegion,
        rectG  :: RectRegion;
        throw  :: Bool = true
    ) -> Bool

Check if a child GeoRegion defined by `Child` is within a RectRegion `rectG`.

Arguments
=========

- `Child` : A GeoRegion that we postulate to be a "child", or a subset of the GeoRegion defined by `polyG`
- `polyG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`

Keyword Arguments
=================

- `throw`  : If `true`, then if `Child` is not within `polyG`, an error is thrown and the program stops running
"""
function isinGeoRegion(
    Child :: GeoRegion,
    rectG :: RectRegion;
    throw :: Bool = true
)

    @info "$(modulelog()) - Performing a check to determine if the $(Child.name) GeoRegion ($(Child.regID)) is inside the $(rectG.name) GeoRegion ($(rectG.regID))"

    isin = isgridinregion(
        [Child.N,Child.S,Child.E,Child.W],
        [rectG.N,rectG.S,rectG.E,rectG.W],
        throw=false
    )

    if !isin

        if throw
            error("$(modulelog()) - The GeoRegion $(Child.regID) ($(Child.name)) is not a subset of the GeoRegion $(rectG.regID) ($(rectG.name))")
        else
            @warn "$(modulelog()) - The GeoRegion $(Child.regID) ($(Child.name)) is not a subset of the GeoRegion $(rectG.regID) ($(rectG.name))"
            return false
        end

    else
        @info "$(modulelog()) - The GeoRegion $(Child.regID) ($(Child.name)) is indeed a subset of the GeoRegion $(rectG.regID) ($(rectG.name))"
        return true
    end

end


text:
function set_res(Delta_Qcon_hat, Delta_Qcon_hat_temp, cellxmax, cellymax, cellzmax)
    res = zeros(cellxmax, cellymax, cellzmax, 5)

    for i in 2:cellxmax-1
        for j in 2:cellymax-1
            for k in 2:cellzmax-1
                for l in 1:5
                    res[i,j,k,l] = abs(Delta_Qcon_hat[i,j,k,l]-Delta_Qcon_hat_temp[i,j,k,l])
                end
            end
        end
    end

    # println(res[3,3,:])
    # println(dt*RHS[3,3,:])
    return res
end

function check_converge(res, RHS, cellxmax, cellymax,cellzmax, init_small)
    norm2 = zeros(5)

    tempAxb = zeros(5)
    tempb = zeros(5)
    for i in 2:cellxmax-1
        for j in 2:cellymax-1
            for k in 2:cellzmax-1
                for l in 1:5
                    tempAxb[l] = tempAxb[l] + res[i,j,k,l]^2
                    tempb[l] = tempb[l] + RHS[i,j,k,l]^2
                end
            end
        end
    end

    #println(tempb)
    #println(tempAxb)
    #throw(UndefVarError(:x))

    for l in 1:5
        norm2[l] = (tempAxb[l]/(tempb[l]+init_small)) ^ 0.5
    end
    return norm2
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
module ElectricFields

using LinearAlgebra
using StaticArrays

using Unitful
using UnitfulAtomic

# using FFTW
using ForwardDiff
using Optim

using IntervalSets
import IntervalSets: duration

using Parameters

import Base: show
using Formatting

include("units.jl")
include("rotations.jl")
include("relation_dsl.jl")

include("field_types.jl")
include("time_axis.jl")
include("carriers.jl")
include("envelopes.jl")
include("arithmetic.jl")

include("field_dsl.jl")

# include("sellmeier.jl")
# include("dispersed_fields.jl")

include("strong_field_properties.jl")

end # module


text:
export iauSxp
"""
Multiply a p-vector by a scalar.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   s      double        scalar
   p      double[3]     p-vector

Returned:
   sp     double[3]     s * p

Note:
   It is permissible for p and sp to be the same array.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauSxp(double s, double p[3], double sp[3])

function iauSxp(s::Real, p::AbstractVector{<:Real})
   sp = zeros(Float64, 3)

   ccall((:iauSxp, libsofa_c), Cvoid,
       (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Float64, s),
       convert(Array{Float64, 1}, p),
       sp)

   return SVector{3}(sp)
end


text:
# This file was generated, do not modify it. # hide
y, X = unpack(caravan, ==(:Purchase), col->true)

std = machine(Standardizer(), X)
fit!(std)
Xs = transform(std, X)

var(Xs[:,1])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
module StochasticVolatility

using ArgCheck
using Distributions
using Parameters
using DynamicHMC
using StatsBase
using Base.Test
using ContinuousTransformations

export StochasticVolatility

"""
    simulate_stochastic(ρ, σ_v, ϵs, νs)

Take in the parameter values (ρ, σ) for the latent volatility process, the errors ϵs used for the measurement equation and the errors νs used for the transition equation.

The discrete-time version of the Ornstein-Ulenbeck Stochastic - volatility model:

    y_t = x_t + ϵ_t where ϵ_t ∼ χ^2(1)
    x_t = ρ * x_(t-1) + σ * ν_t  where ν_t ∼ N(0, 1)

"""
function simulate_stochastic(ρ, σ, ϵs, νs)
    N = length(ϵs)
    @argcheck N == length(νs)
    x₀ = νs[1]*σ*(1 - ρ^2)^(-0.5)
    xs = Vector{typeof(x₀)}(N)
    for i in 1:N
        xs[i] = (i == 1) ? x₀ : (ρ*xs[i-1] + σ*νs[i])
    end
    xs + log.(ϵs) + 1.27
end

simulate_stochastic(ρ, σ, N) = simulate_stochastic(ρ, σ, rand(Chisq(1), N), randn(N))

struct StochasticVolatility{T, Prior_ρ, Prior_σ, Ttrans}
    "observed data"
    ys::Vector{T}
    "prior for ρ (persistence)"
    prior_ρ::Prior_ρ
    "prior for σ_v (volatility of volatility)"
    prior_σ::Prior_σ
    "χ^2 draws for simulation"
    ϵ::Vector{T}
    "Normal(0,1) draws for simulation"
    ν::Vector{T}
    "Transformations cached"
    transformation::Ttrans
end

## In this form, we use an AR(2) process of the first differences with an intercept as the auxiliary model.

"""
    lag(xs, n, K)

Lag-`n` operator on vector `xs` (maximum `K` lags).
"""
lag(xs, n, K) = xs[((K+1):end)-n]

"""
    lag_matrix(xs, ns, K = maximum(ns))

Matrix with differently lagged xs.
"""
function lag_matrix(xs, ns, K = maximum(ns))
    M = Matrix{eltype(xs)}(length(xs)-K, maximum(ns))
    for i ∈ ns
        M[:, i] = lag(xs, i, K)
    end
    M
end

"first auxiliary regression y, X, meant to capture first differences"
function yX1(zs, K)
    Δs = diff(zs)
    lag(Δs, 0, K), hcat(lag_matrix(Δs, 1:K, K), ones(eltype(zs), length(Δs)-K), lag(zs, 1, K+1))
end

"second auxiliary regression y, X, meant to capture levels"
function yX2(zs, K)
    lag(zs, 0, K), hcat(ones(eltype(zs), length(zs)-K), lag_matrix(zs, 1:K, K))
end

"Constructor which calculates cached values and buffers. Use this unless you know what you are doing."
ℝto(dist::Uniform) = transformation_to(Segment(minimum(dist), maximum(dist)), LOGISTIC)
ℝto(::InverseGamma) = transformation_to(ℝ⁺)

"""
    StochasticVolatility(ys, prior_ρ, prior_σ, M)

Convenience constructor for StochasticVolatility.
Take in the observed data, the priors, and number of simulations (M).
"""

function StochasticVolatility(ys, prior_ρ, prior_σ, M)
    transformation = TransformationTuple((ℝto(prior_ρ), (ℝto(prior_σ))))
    StochasticVolatility(ys, prior_ρ, prior_σ, rand(Chisq(1), M), randn(M), transformation)
end

function (pp::StochasticVolatility)(θ)
    @unpack ys, prior_ρ, prior_σ, ν, ϵ, transformation = pp
    ρ, σ = transformation(θ)
    logprior = logpdf(prior_ρ, ρ) + logpdf(prior_σ, σ)
    N = length(ϵ)

    # Generating xs, which is the latent volatility process

    xs = simulate_stochastic(ρ, σ, ϵ, ν)
    Y_1, X_1 = yX1(xs, 2)
    β₁ = qrfact(X_1, Val{true}) \ Y_1
    v₁ = mean(abs2,  Y_1 - X_1*β₁)
    Y_2, X_2 = yX2(xs, 2)
    β₂ = qrfact(X_2, Val{true}) \ Y_2
    v₂ = mean(abs2,  Y_2 - X_2*β₂)
    # We work with first differences
    y₁, X₁ = yX1(ys, 2)
    log_likelihood1 = sum(logpdf.(Normal(0, √v₁), y₁ - X₁ * β₁))
    y₂, X₂ = yX2(ys, 2)
    log_likelihood2 = sum(logpdf.(Normal(0, √v₂), y₂ - X₂ * β₂))
    logprior + log_likelihood1 + log_likelihood2
end

"""

    path(parts...)

`parts...` appended to the directory containing `src/` etc. Useful for saving graphs and data.

"""

path(parts...) = normpath(joinpath(splitdir(Base.find_in_path("StochasticVolatility"))[1], "..", parts...))

end


text:
using ArgParse
using JLD2

parser = ArgParseSettings(allow_ambiguous_opts=false)
@add_arg_table! parser begin
    "main"
        required = true
    "patch"
        required = true
    "--overwrite"
        arg_type = Bool
        default = true
end
args = parse_args(parser)

results = load(args["main"])["results"];
patch = load(args["patch"])["results"];
for (noise, sigscan) in patch
    for (signal, data) in sigscan
        (signal in keys(results[noise])) == args["overwrite"] || error("check overwrite for results[$noise][$signal]")
        results[noise][signal] = data
        println("including results[$noise][$signal]")
    end
end
save(args["main"], "results", results)


text:
using OrdinaryDiffEq, Zygote, DiffEqSensitivity, Test

p_model = [1f0]

u0 = Float32.([0.0])

function dudt(du, u, p, t)
    du[1] = p[1]
end

prob = ODEProblem(dudt,u0,(0f0,99.9f0))

function predict_neuralode(p)
    _prob = remake(prob,p=p)
    Array(solve(_prob,Tsit5(), saveat=0.1))
end

loss(p) = sum(abs2,predict_neuralode(p))/length(p)

p_model_ini = copy(p_model)

@test !iszero(Zygote.gradient(loss,p_model_ini)[1])