text:
import ImageMagick
using HTTP
using Test
using FileIO
using Scratch
using LazyArtifacts

const datadir = joinpath(artifact"data")
@info "The test data is located at $datadir."

const imgdir  = @get_scratch!("img")
@info "If you want to check the output of the tests, please head to $imgdir."

function exportImage(filename, I::AbstractMatrix)
  Iabs = abs.(I)
  Icolored = colorview(Gray, Iabs./maximum(Iabs))
  save(filename, Icolored )
end

include("LoadSaveMDF.jl")
include("Reconstruction.jl")
include("Cartesian.jl")
include("MotionCompensation.jl")
include("MultiPatch.jl")
include("MultiGradient.jl")
include("Sparse.jl")
include("SMCenter.jl")
include("SMRecovery.jl")


text:
using CSV, DataFrames, Dates, MessyTimeSeries;
using FredData: get_data, Fred;

"""
    get_h10_dataset(fred_tickers::Array{String,1}, mnemonics::Array{String,1}, transform::IntVector, to_include::IntVector, observation_start::String, observation_end::String)

Get H.10 table from the Federal Reserve Economic Data (FRED) at weekly frequency and in log-return.
"""
function get_h10_dataset(fred_tickers::Array{String,1}, mnemonics::Array{String,1}, transform::IntVector, to_include::IntVector, observation_start::String, observation_end::String)

    f = Fred();
    df = DataFrame();

    for (i, ticker) in enumerate(fred_tickers)

        # Download current series
        df_current = get_data(f, ticker, observation_start=observation_start, observation_end=observation_end, frequency="wef", aggregation_method="eop", units="cch").data::DataFrame;
        df_current = df_current[!,[:date, :value]];
        if transform[i] == 1
            df_current[!,:value] *= -1;
        end

        # Rename `:value`
        rename!(df_current, Dict(:value => mnemonics[i]));

        # Add series to `df`
        if i == 1
            df = df_current;
        else
            df = outerjoin(df, df_current, on=:date);
        end
    end

    # Use selected exchange rates
    countries_selection = names(df)[2:end][to_include .== 1];
    df = df[!, vcat(names(df)[1], countries_selection)];

    # Chronological order
    sort!(df, :date);

    # Return output
    return df;
end

# Load info from CSV
fred_tickers = CSV.read("./data/fred_h10_data.csv", DataFrame)[!,:fred_ticker] |> Vector{String};
mnemonics = CSV.read("./data/fred_h10_data.csv", DataFrame)[!,:mnemonic] |> Vector{String};
transform = CSV.read("./data/fred_h10_data.csv", DataFrame)[!,:transform];
to_include = CSV.read("./data/fred_h10_data.csv", DataFrame)[!,:include];

ext: jl


text:
include("functions.jl")

hh = Household();
r_guess = 0.01;
cpol,apol = iterate_egm(hh;r=r_guess)

plot(hh.Amat[:,1], hh.Amat[:,1], label = "", color = :black, linestyle = :dash)
plot!(hh.Amat[:,1], apol[:,1], label = "ϵ = $(round(hh.Ymat[1,1], digits = 3))")
plot!(hh.Amat[:,1], apol[:,5], label = "ϵ = $(round(hh.Ymat[1,5], digits = 3))")

r = market_clearing(hh,A_guess,r = r_guess)
plot_market_clearing(hh)


text:
"""
    subnames(name2var, ast)

Subsitute the values from the dict `name2var`
into the ast,
where ever the keys (names) occur.
"""
subnames(name2var, x::Any) = x # fallback, for literals
subnames(name2var, name::Symbol) = get(name2var, name, name)
function subnames(name2var, code::Expr)
    MacroTools.postwalk(code) do name
        # TODO Make this handle the name being on the LHS of assigment right
        get(name2var, name, name)  # If we have a value for it swap it in.
    end
end

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

function get_user_input(io=stdin)
    printstyled("debug> "; color=:light_red)

    ast = nothing
    line = ""
    while true
        # Very cut down REPL input code
        # https://github.com/JuliaLang/julia/blob/b8c0ec8a0a2d12533edea72749b37e6089a9d163/stdlib/REPL/src/REPL.jl#L237
        line *= better_readline(io)
        ast = Base.parse_input_line(line)
        ast isa Expr && ast.head == :incomplete || break
    end
    return ast
end

function eval_and_display(code_ast, eval_module)
    try
        res = eval_module.eval(code_ast)
        if res != nothing
            display(res)
        end
    catch err
        printstyled("ERROR: ", color=:red)
        showerror(stdout, err)
    end
    println()
end

function run_repl(name2var, eval_module)
    local code_ast
    while true
        code_ast = get_user_input()
        if Meta.isexpr(code_ast, :toplevel, 2) && haskey(actions, code_ast.args[2])
            return code_ast.args[2] # Send the codeword back
        end
        code_ast = subnames(name2var, code_ast)
        eval_and_display(code_ast, eval_module)
   end
end


text:
#   Unit "tapers" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements several tapering windows in the time domain.
#
#   See F.J. Harris
#   "On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform"
#   Proc. IEEE, 66, 51-53, 1978
#   http://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf
#
#   D. Slepian
#   "Prolate Spheroidal Wave Functions. Fourier Analysis, and Uncertainty—V: The Discrete Case"
#   The Bell System Technical Journal,VoL 57, No. 5. May-June 1978
#   https://ieeexplore-ieee-org.gaelnomade-1.grenet.fr/stamp/stamp.jsp?tp=&arnumber=6771595
#
#   D.J. Thomson
#   "Spectrum estimation and harmonic analysis."
#   Proc. IEEE 70: 1055-1096, 1982.
#   http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.471.1278&rep=rep1&type=pdf

#   Other resources :
#   https://pdfs.semanticscholar.org/752d/1a551b96559458064323eb3de7faaaef4c4e.pdf
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

# time-domain tapering windows.
@enum TaperKind begin
    rectangular = 1
    triangular  = 2
    hann        = 3
    hamming     = 4
    blackman    = 5
    harris4     = 6
    riesz       = 7
    parzen      = 8
    slepian     = 9
end

TaperData=Union{Vector{T}, Matrix{T}} where T<:Union{Real, Complex}

struct Taper
    y    :: TaperData # here is the taper
    kind :: TaperKind # enumerated type; see here above
    α    :: Real # 1-parameter tapers are supported
    n    :: Int # number of tapers; >1 only for slepians
end

### Time-domain tapering ###
###############################################################################
"""
```julia
function taper( kind  :: TaperKind,
                wl    :: Int;
            α       :: Real    = 2.,
            n       :: Int     = ceil(Int, 2*α)-1,
            padding :: Int     = 0,
            type    :: Type{T} = Float64) where T<:Union{Real, Complex}

```
Universal constructor of [Taper](@ref) objects, given a tapering window
`kind`, of type [TaperKind](@ref) and the window length `wl`.
NB: for the `hamming` and `blackman` type use `FourierAnalysis.hamming` and `FourierAnalysis.blackman`
to avoid conflicts with the DSP package.

Return a vector of length `wl` for all types of tapers, but for the dpss
(Slepian multi-tapers), for which return a matrix of size `wl` x `n`.

If optional keyword argument `padding` is >0, then the actual window length
will be `wl` + `padding`.

The `type` optional keyword argument can be use to specify the the type of
elements of the typering window. By default, this is the `Float64` type.

Optional keywords arguments `α` and `n` currently apply only for slepian
multi-tapering (*discrete prolate spheroidal sequences*):

`α` is the *half bandwidth* (hbw) parameter as per the DSP.jl package.
This unit is used in many dpss implementations and it is often reported
that "typical values" for `α` are 2, 2.5, 3, 3.5 or 4.
However, the optimal smoothing increases with the ratio between window length
and sampling rate, thus these values in absolute terms are not useful.
In fact, the larger the *hbw* and the higher `n`, the smoother the
spectra will be (variance reduction). In order to overcome these
difficulties, for Slepian's multitapering
*FourirAnalysis* implements the [`slepians`](@ref) constructor,
which allows the bandwidth parameter to be given in Hz and the
number of tapering windows to be chosen automatically.

`n` is the number of tapering windows. For slepian tapers this is
the number of the discrete prolate spheroidal sequences.
As in the DSP package, by default it is set to
`ceil(Int, 2*α)-1`, however, depending on how large this number is,
low eigenvalues may correspond to the last sequences, therefore
those should be discarded.

**See**: [plot tapering windows](@ref).

**See also**: [`slepians`](@ref), [`taperinfo`](@ref)

**Examples**:
```julia
using FourierAnalysis

## Use the constructor
sr, t, f, a = 128, 128, 10, 0.5
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# create a data matrix
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
# compute spectra using hamming tapering window
# we need to prepend 'FourierAnalysis.' since `hamming`
# is declared also in DSP.jl
H=taper(FourierAnalysis.hamming, t)
S=spectra(X, sr, t; tapering=H)
# you can obtain the same thing with
S=spectra(X, sr, t; tapering=H)
# which will create the hamming tapering window on the fly,
# thus calling explicitly the constructor is interesting
# only if you need to reuse the same tapering window many times.

## Plot tapering windows using the standard plot function
using Plots
tapers=[TaperKind(i) for i=1:8]
X=zeros(t, 8)
for i=1:8 X[:, i] = taper(tapers[i], t).y end
mylabels=Array{String}(undef, 1, 8)
for i=1:8 mylabels[1, i]=string(tapers[i]) end
plot(X; labels=mylabels)

## using the recipe declared in recipes.jl
plot(taper(parzen, 256))
plot(taper(slepian, 256, α=4, n=7))

```
"""
function taper( kind :: TaperKind,
                wl   :: Int;  # wl is the window length
            α        :: Real    = 2.,    # only for splepians;
            n        :: Int     = ceil(Int, 2*α)-1, #ceil : Int greater then or equal to
            padding  :: Int     = 0,
            type     :: Type{T} = Float64) where T<:Union{Real, Complex}

    if      kind ∉ (rectangular, slepian)
            wl₋₁     = wl-1
            wl₋₁½    = wl₋₁/2
            c(i::Int) = cos(2π*(i-1)/wl₋₁)
            c(i::Int, j::Int) = cos(2π*j*(i-1)/wl₋₁)
    end

    if      kind == rectangular v = ones(type, wl)
    elseif  kind == triangular  v = Vector{type}([1 - abs(2*(i-1)-wl₋₁)/wl₋₁ for i=1:wl]) # zero at the edge
    elseif  kind == hann        v = Vector{type}([0.5*(1-c(i)) for i=1:wl])
    elseif  kind == hamming     v = Vector{type}([0.54 - 0.46*c(i) for i=1:wl])
    elseif  kind == blackman    v = Vector{type}([0.42 + 0.50*-c(i) + 0.08*c(i, 2) for i=1:wl])
                                     v[1]=0.; v[end]=0.
    elseif  kind == harris4     v = Vector{type}([0.35875 - 0.48829*c(i) + 0.14128*c(i, 2) - 0.01168*c(i, 3) for i=1:wl])
    elseif  kind == riesz       v = Vector{type}([1 - ((i-1-wl₋₁½)/wl₋₁½)^2 for i=1:wl])
    elseif  kind == parzen
            v=Vector{type}(undef, wl)
            wl½, wl¼= wl÷2, wl÷4
            for i=1:wl¼        v[i] = 2*((1 - (abs(i-wl½)/wl½))^3) end
            for i=wl¼+1:wl½    v[i] = 1 - 6*(abs(i-wl½)/wl½ )^2 + 6*(abs(i-wl½)/wl½)^3 end
            for i=wl½+1:wl-wl¼ v[i] = 1 - 6*(abs(i-1-wl½)/wl½)^2 + 6*(abs(i-1-wl½)/wl½)^3 end
            for i=wl-wl¼+1:wl  v[i] = 2*((1 - (abs(i-1-wl½)/wl½))^3) end
    elseif  kind == slepian    v    = Matrix{type}(dpss(wl, α, n, zerophase=false))
    end

    # This has to be checked:

    # normalize to unit mean or unit mean of absolute values
    if      kind ∉ (rectangular, slepian) v ./= mean(v) end
    if      kind == slepian v./=mean(abs.(a) for a in v) end

    # if      kind == slepian v./=mean(v.^2; dims=1) end
    # if      kind ∉ (rectangular, slepian) v ./= mean(v) end

    # normalization in DSP.jl
    #    if      taper == slepian for i=1:size(v, 2) v[:, i]./= (norm(v[:, i])/wl) end end
    # as above but givin less and less weight to the sequences
    #    if      taper == slepian for i=1:size(v, 2) v[:, i]./= √i*(norm(v[:, i])/wl) end end

    p=padding
    p>0 && (kind == slepian ? v=[v; zeros(type, p, n)] : v=[v; zeros(type, p)])

    return kind == slepian ? Taper(v, kind, α, n) : Taper(v, kind, 0., 1)
end

"""
Construct a [Taper](@ref) objects holding Slepian's multi-tapering
*discrete prolate spheroidal sequences*,
given sampling rate `sr`, window length `wl` and
the `bandwidth` argument in Hz.
For EEG data, 1<=bandwidth<=2 is an adequate choice.

The 'half-bandwidth' parameter `α` used in the DSP package and in the
universal [Taper](@ref) constructor is set as

        `α=(bandwidth/2)*wl/sr`.

The optimal number of dpss is heuristically set to

        `n=max(1, trunc(Int, 2*α)-trunc(Int, log(2*α)))`.

The created object can be passed as argument
in constructors [`spectra`](@ref), [`crossSpectra`](@ref) and
[`coherence`](@ref).

**See**: [plot tapering windows](@ref).

**Examples**:
```julia
using FourierAnalysis
sr, t, f, a = 128, 128, 10, 0.5
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# create a data matrix
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
# compute spectra using slepian multi-tapering with bandwidth 1.5
H=slepians(sr, t, 2)
S=spectra(X, sr, t; tapering=H)

using Plots
plot(H)
plot(S)
```
"""
function slepians( sr    :: Int,
                   wl    :: Int,
               bandwidth :: Real = 1.5)
    α=(bandwidth/2)*wl/sr # α parameter
    # heuristic to eliminate eigenfunctions with small eigenvalues
    n=max(1, trunc(Int, 2*α)-trunc(Int, log(2*α)))
    return taper(slepian, wl, α=α, n=n)
end
###############################################################################

"""
```julia
function taperinfo(taper::Taper)
```

Return the name of the
tapering window(s) encapsulated in the [Taper](@ref) object
as a string.

Only for Slepian's discrete prolate spheroidal sequences (dpss),
their parameters, namely, ``α`` (half-bandwidth)
and ``n`` (number of windows), are reported within parentheses as well.

**Examples**:
```julia
H=taper(hamming, 128*8)
taperinfo(H)

H=slepians(128, 128*8, 2)
taperinfo(H)
```
"""
taperinfo(taper::Taper) =
     taper.kind==slepian ? string(taper.kind)*"'s dpss (alpha=$(taper.α), n=$(taper.n))" :
                           string(taper.kind)

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, 𝜏::Taper)
println(io, titleFont, "⍓ Taper type; $(size(𝜏.y, 1))-samples")
#println(io, "□  □    □      □        □           □", defaultFont)
println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
println(io, "taper kind   (.kind): ", string(𝜏.kind))
println(io, "half-bandwidth  (.α): $(𝜏.α)")
println(io, "number of tapers(.n): $(𝜏.n)")
𝜏.kind==slepian ? println(io, "data            (.y): $(size(𝜏.y, 1))x$(size(𝜏.y, 2))-", typeof(𝜏.y)) :
                 println(io, "data            (.y): $(length(𝜏.y))-", typeof(𝜏.y))

end
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #


text:
# can I make other packages add to this list?
GeodesyTypes = [LLA,
                ECEF,
                LL,
                ENU]

###########################################
# defines methods to add to point types
# Note the "geotransform" function handles all transformations
# while the "convert" function converts types leaving values unchanged
###########################################

# add some general versions of functions to overload
get_refloc{T}(::Type{T}) = error("Can't determine the reference location for type $(T)")
get_datum{T}(::Type{T}) = error("Can't determine the datum for type $(T)")
get_srid{T}(::Type{T}) = error("Can't determine the SRID for type $(T)")
get_geoid{T}(::Type{T}) = error("Can't determine the geoid for type $(T)")

# generic version of add parameters
add_param{T}(::Type{T}) = T
add_param{T}(::Type{T}, X) = T

###################################################################
# Function to add an integer indexing method for the geodesy types
###################################################################

# and get index methods to each type
function add_indexing{geodesy_type}(::Type{geodesy_type})
    if geodesy_type <: FixedArray
        q = quote; end  # no need for subtypes of this
    else
        q = quote
            # access by index
            @inline getindex{T <: $(geodesy_type)}(x::T, i::Integer) = x.(i)
        end
    end
end

##############################################################
# Function to add accessor methods to Geodesy point types
##############################################################

# field input should be the field names of whatever best resemble X, Y, and Z (if it exists)
function add_accessors{geodesy_type}(::Type{geodesy_type}, fields)

    #
    # Get horizontal locations
    #
    # N.B. fields should ordered in x / y / z fashion
    q = quote

        # add horizontal accessors
        @inline getX(X::$(geodesy_type)) = X.$(fields[1])
        @inline getY(X::$(geodesy_type)) = X.$(fields[2])

        @inline get_lon(X::$(geodesy_type)) = X.$(fields[1])
        @inline get_lat(X::$(geodesy_type)) = X.$(fields[2])

        @inline get_east(X::$(geodesy_type)) = X.$(fields[1])
        @inline get_north(X::$(geodesy_type)) = X.$(fields[2])

    end

    #
    # add vertical as well
    #
    if ((length(fields)) > 2)
        qn = quote
            @inline getZ(X::$(geodesy_type)) = X.$(fields[3])
            @inline get_alt(X::$(geodesy_type)) = X.$(fields[3])
            @inline get_up(X::$(geodesy_type)) = X.$(fields[3])
        end
    else
        # for legacy support return 0 when we try to get the height of a surface point
        qn = quote
            @inline getZ(X::$(geodesy_type)) = 0.0
            @inline get_alt(X::$(geodesy_type)) = 0.0
            @inline get_up(X::$(geodesy_type)) = 0.0
        end
    end
    append!(q.args, qn.args)

    #
    # get the default parameters for this type
    #
    defaults = default_params(geodesy_type)

    #
    # get the reference location for the type if it has one
    #
    if has_refloc(geodesy_type) == Val{true}
        if (length(geodesy_type.parameters) == 1)  # assume its in here directly.  This will cause problems in the future
            qn = quote

                # get the reference location (note: for the moment it should always be a world surface point, but its not templated that way for convenience)
                get_refloc(::Type{$(geodesy_type)})                         =  $(defaults[1])      # it was not supplied
                get_refloc{T}(::Type{$(geodesy_type){T}})                   =  T                   # it was already supplied, but force it to be the LL type

                # get it from the point as well
                get_refloc{T}(::$(geodesy_type){T})                         =  T

            end
        else
            # need to manually define the Geodesy.get_refloc functions for this type
            qn = quote; end
        end
        append!(q.args, qn.args)
    end

    #
    # get the datum / ellipse for the type if it has one
    #
    if has_ellipse(geodesy_type) == Val{true}

        if (geodesy_type.parameters[1] <: AbstractDatum)  # points parameterised by datums
            qn = quote

                # extract from the type
                get_datum(::Type{$(geodesy_type)})       =  $(defaults[1])          # it was not supplied
                get_datum{T}(::Type{$(geodesy_type){T}}) =  T                       # it was already supplied

                # extract from the point
                get_datum{T}(::$(geodesy_type){T})       =  T

            end
        elseif (length(geodesy_type.parameters) == 1)  # it must be in here
            qn = quote

                # extract from the type
                get_datum(::Type{$(geodesy_type)})       =  UnknownDatum        # the reference point wasn't supplied so we can't get it
                get_datum{T}(::Type{$(geodesy_type){T}}) =  get_datum(T)        # it was already supplied

                # extract from the point
                get_datum{T}(::$(geodesy_type){T})       =  get_datum(T)
            end
        else
            # need to manually define the Geodesy.get_datum() functions for this type
            qn = quote; end
        end
        append!(q.args, qn.args)

        # get the actual ellipse corresponding to the above
        qn = quote
            ellipsoid{T <: $(geodesy_type)}(::Type{T})  = ellipsoid(get_datum(T))    # reference ellipsoid for the position
            ellipsoid{T <: $(geodesy_type)}(X::T)       = ellipsoid(get_datum(T))    # reference ellipsoid for the position
        end
        append!(q.args, qn.args)
    end

    #
    # get the srid for the type if it has one
    #
    if has_srid(geodesy_type) == Val{true}
        if (geodesy_type.parameters[1] <: AbstractSRID)
            qn = quote

                # get the parameterising datum / ellipse
                get_srid(::Type{$(geodesy_type)})       =  $(defaults[1])        # it was not supplied
                get_srid{T}(::Type{$(geodesy_type){T}}) =  T                     # it was already supplied

                # extract from the point
                get_srid{T <: $(geodesy_type)}(::T)     = get_srid(T)

                # also add the SRID contructor from this type
                call(::Type{SRID}, ::Type{$(geodesy_type)})         =  $(defaults[1])
                call{T <: $(geodesy_type)}(::Type{SRID}, ::Type{T}) =  get_srid(T)
                call{T <: $(geodesy_type)}(::Type{SRID}, ::T)       =  get_srid(T)
            end

            # add extra parameter stuff if this type has an extra parameter
            if (length(geodesy_type.parameters) > 1)
                qa = quote

                    # get the parameterising datum / ellipse
                    get_srid{T, G}(::Type{$(geodesy_type){T,G}}) =  T                     # it was already supplied

                end
                append!(qn.args, qa.args)
            end
        else
            # need to manually define the Geodesy.get_srid() function for this type
            qn = quote; end
        end
        append!(q.args, qn.args)
    end

    #
    # get the geoid for the type if it has one
    #
    if has_geoid(geodesy_type) == Val{true}
        if (geodesy_type.parameters[1] <: AbstractGeoid)
            qn = quote

                # get the parameterising geoid
                get_geoid(::Type{$(geodesy_type)})       =  $(defaults[1])   # it was not supplied
                get_geoid{T}(::Type{$(geodesy_type){T}}) =  T                # it was already supplied

                # extract from the point
                get_geoid{T}(::$(geodesy_type){T})       = T

            end
        elseif (geodesy_type.parameters[2] <: AbstractGeoid)
            qn = quote

                # get the parameterising geoid
                get_geoid(::Type{$(geodesy_type)})           =  $(defaults[2])      # it was not supplied
                get_geoid{T}(::Type{$(geodesy_type){T}})     =  $(defaults[2])      # it was not supplied
                get_geoid{T,G}(::Type{$(geodesy_type){T,G}}) =  G                   # it was supplied

                # extract from the point
                get_geoid{T,G}(::$(geodesy_type){T,G})       = G
            end
        else
            # need to manually define the Geodesy.get_geoid() function for this type
            qn = quote; end
        end
        append!(q.args, qn.args)
    end

    return q
end

##############################################################
# Function to add the template parameter to a type
##############################################################

# general version
function add_param_methods{geodesy_type}(::Type{geodesy_type})

    # get the default parameters for this type
    defaults = default_params(geodesy_type)

    # and build the code
    if (length(geodesy_type.parameters) == 1)
        quote
            add_param(::Type{$(geodesy_type)})        = $(geodesy_type){$(defaults[1])}    # it was not supplied
            add_param{T}(::Type{$(geodesy_type){T}})  = $(geodesy_type){T}                 # it was already supplied
        end
    elseif (length(geodesy_type.parameters) == 2)
        quote
            add_param(::Type{$(geodesy_type)})           = $(geodesy_type){$(defaults[1]), $(defaults[2])}    # it was not supplied
            add_param{T}(::Type{$(geodesy_type){T}})     = $(geodesy_type){T, $(defaults[2])}                 # it was half supplied
            add_param{T,U}(::Type{$(geodesy_type){T,U}}) = $(geodesy_type){T, U}                              # it was already supplied
        end
    else
        # need to manually define the Geodesy.add_param() function for this type
        quote; end
    end
end

#####################################################################
### Methods to infer template parameters based on other arguments ###
#####################################################################

#
# Function to add the template parameter to a type
#
function add_infer_param_methods{type_1, type_2}(::Type{type_1}, ::Type{type_2})

    # by default, insert the default value.  Weird.
    insert_param = Vector{Any}([x for x in default_params(type_1)])

    # and see what we can find
    for i = 1:length(type_1.parameters)

        # what's the parameter, and does the other input have it?
        if (type_1.parameters[i] <: AbstractSRID) && (has_srid(type_1) == Val{true})

            # can infer the SRID
            insert_param[i] = :(get_srid(X))

        elseif (type_1.parameters[i] <: AbstractGeoid) && (has_geoid(type_2) == Val{true})

            # can infer the Geoid
            insert_param[i] = :(get_geoid(X))

        elseif (type_1.parameters[i] <: WorldPosition) && (has_refloc(type_2) == Val{true})

            # can infer the reference location
            insert_param[i] = :(get_geoid(X))

        elseif (type_1.parameters[i] <: AbstractDatum) && (has_ellipse(type_2) == Val{true})

            # can infer the datum
            insert_param[i] = :(get_datum(X))
        end
    end

    # and build the code
    if (length(type_1.parameters) == 1)
        quote
            add_param{X <: $(type_2)}(::Type{$(type_1)}, ::Type{X})            = $(type_1){$(insert_param[1])}    # it was not supplied
            add_param{T, X <: $(type_2)}(::Type{$(type_1){T}}, ::Type{X})      = $(type_1){T}                     # it was already supplied
            add_param{T <: $(type_1), X <: $(type_2)}(::Type{T}, ::X)          = add_param(T, X)
        end
    elseif (length(type_1.parameters) == 2)
        quote
            add_param{X <: $(type_2)}(::Type{$(type_1)}, ::Type{X})            = $(type_1){$(insert_param[1]), $(insert_param[2])}    # it was not supplied
            add_param{T, X <: $(type_2)}(::Type{$(type_1){T}}, ::Type{X})      = $(type_1){T, $(insert_param[2])}                     # it was half supplied
            add_param{T,U, X <: $(type_2)}(::Type{$(type_1){T, U}}, ::Type{X}) = $(type_1){T, U}                                      # it was already supplied
            add_param{T <: $(type_1), X <: $(type_2)}(::Type{T}, ::X)          = add_param(T, X)
        end
    else
        # need to overload Geodesy.add_param() manually
        #error("Haven't defined an add parameters method for more than 2 parameters (type: $(type_1))")
    end
end

#####################################################
# Make the constructors force typing
#####################################################

function add_constructors{geodesy_type}(::Type{geodesy_type}, fields)

    nfields = length(fields)

    q = quote

        # copy constructor
        call{T <: $(geodesy_type)}(::Type{T}, X::T) = X
        call{T <: $(geodesy_type)}(::Type{$(geodesy_type)}, X::T) = X
        call{T <: $(geodesy_type), U <: $(geodesy_type)}(::Type{T}, X::U) = geotransform(T, X, get_handler($(geodesy_type)), get_handler($(geodesy_type)))

    end

    #
    # allow construction to convert from no datum to datum
    #

    # an expression to form a tuple from the type field
    rhs_expr = :(())
    append!(rhs_expr.args, [:(X.$(field)) for field in fields])

    # default parameters
    defaults = Geodesy.default_params(geodesy_type)
    abstracts = [super(def) for def in defaults]

    if (length(geodesy_type.parameters) == 1)
        qn = quote
            convert(::Type{$(geodesy_type){$(defaults[1])}}, X::$(geodesy_type){$(defaults[1])}) = X  # need this because it can't be resolved below
            convert{T <: $(abstracts[1])}(::Type{$(geodesy_type){T}}, X::$(geodesy_type){$(defaults[1])}) = $(geodesy_type){T}($(rhs_expr.args...))  # convert from the null parameter
            convert{T <: $(abstracts[1])}(::Type{$(geodesy_type){$(defaults[1])}}, X::$(geodesy_type){T}) = $(geodesy_type){$(defaults[1])}($(rhs_expr.args...))  # convert to the null parameter
        end
        append!(q.args, qn.args)
    elseif (length(geodesy_type.parameters) == 2)
        qn = quote

            # define some specifics to prevent ambiguities below
            convert(::Type{$(geodesy_type){$(defaults[1]), $(defaults[2])}}, X::$(geodesy_type){$(defaults[1]), $(defaults[2])}) = X
            convert{T <: $(abstracts[1])}(::Type{$(geodesy_type){T, $(defaults[2])}}, X::$(geodesy_type){T, $(defaults[2])}) = X
            convert{U <: $(abstracts[2])}(::Type{$(geodesy_type){$(defaults[1]), U}}, X::$(geodesy_type){$(defaults[1]), U}) = X

            # conversion to / from fully null parameters
            convert{T <: $(abstracts[1]), U <: $(abstracts[2])}(::Type{$(geodesy_type){T, U}}, X::$(geodesy_type){$(defaults[1]), $(defaults[2])}) =
                $(geodesy_type){T, U}($(rhs_expr.args...))                            # convert from the null parameter
            convert{T <: $(abstracts[1]), U <: $(abstracts[2])}(::Type{$(geodesy_type){$(defaults[1]), $(defaults[2])}}, X::$(geodesy_type){T,U}) =
                $(geodesy_type){$(defaults[1]), $(defaults[2])}($(rhs_expr.args...))  # convert to the null parameter

            # conversion to / from when the 1st parameter is known
            convert{T <: $(abstracts[1]), U <: $(abstracts[2])}(::Type{$(geodesy_type){T, U}}, X::$(geodesy_type){$(defaults[1]), U}) =
                $(geodesy_type){T, U}($(rhs_expr.args...))               # convert from the null parameter
            convert{T <: $(abstracts[1]), U <: $(abstracts[2])}(::Type{$(geodesy_type){$(defaults[1]), U}}, X::$(geodesy_type){T,U}) =
                $(geodesy_type){$(defaults[1]), U}($(rhs_expr.args...))  # convert to the null

            # conversion to / from when the 2nd parameter is known
            convert{T <: $(abstracts[1]), U <: $(abstracts[2])}(::Type{$(geodesy_type){T, U}}, X::$(geodesy_type){T, $(defaults[2])}) =
                $(geodesy_type){T, U}($(rhs_expr.args...))  # convert from the null parameter
            convert{T <: $(abstracts[1]), U <: $(abstracts[2])}(::Type{$(geodesy_type){T, $(defaults[2])}}, X::$(geodesy_type){T,U}) =
                $(geodesy_type){T, $(defaults[2])}($(rhs_expr.args...))  # convert to the null

        end
        append!(q.args, qn.args)
    end

    #
    # construction from scalars
    #

    # create an expression for a tuple x1, x2, x3...
    rhs_expr = :(())
    append!(rhs_expr.args, [:($(fields[i])) for i in 1:nfields])

    # create an expression for a tuple x1::Int, x2::Int, x3::Int... for allowing construction from Ints
    lhs_expr = :(())
    append!(lhs_expr.args, [:($(fields[i])::Real) for i in 1:nfields])

    # and build the constructor
    qn = quote
        call(::Type{$(geodesy_type)}, $(lhs_expr.args...)) = add_param($(geodesy_type))($(rhs_expr.args...))
    end
    append!(q.args, qn.args)

    # if this is an lla type, for legacy reasons we allow construction with / without the height parameter
    if has_alt(geodesy_type) == Val{true}
        if (nfields == 3)
            qn = quote  # construct it wihtout the height input
                call{T <: $(geodesy_type)}(::Type{T}, $(lhs_expr.args[1:2]...)) = add_param(T)($(rhs_expr.args[1:2]...), 0.0)
            end
            append!(q.args, qn.args)
        elseif (nfields == 2)
            qn = quote  # construct it with the height input and ignore the height
                call{T <: $(geodesy_type)}(::Type{T}, $(lhs_expr.args[1:2]...), alt::Real) = add_param(T)($(rhs_expr.args[1:2]...))
            end
            append!(q.args, qn.args)
        end
    end
    return q
end

#############################################################
# Add the cross constructors (constructing from other types)
#############################################################

function add_cross_constructors{type_1, type_2}(::Type{type_1}, ::Type{type_2})

    #h1, h2 = get_handler(type_1), get_handler(type_2)
    qb = quote

        # no reference point
        call{T <: $(type_1)}(::Type{T}, X::$(type_2)) = geotransform(T, X)           # $(h1), $(h2)

        # with reference point
        call{T <: $(type_1)}(::Type{T}, X::$(type_2), ref) = geotransform(T, X, ref) # $(h1), $(h2)

    end

    # define a value preserving transform if they have different handlers
    if (get_handler(type_1) != get_handler(type_2))
        qn = quote
            convert{T <: $(type_1)}(::Type{T}, X::$(type_2)) = add_param(T)(getX(X), getY(X), getZ(X))
        end
        append!(qb.args, qn.args)
    end
    return qb
end

#####################################################
### Add Proj4 projections for each type where known
######################################################

function add_projection{geodesy_type}(::Type{geodesy_type})
    if (has_geoid(geodesy_type) != Val{true})
        quote
            get_projection{T <: $(geodesy_type)}(X::T) = get_projection(get_srid(T))
            get_projection{T <: $(geodesy_type)}(::Type{T}) = get_projection(get_srid(T))
        end
    else
        quote
            get_projection{T <: $(geodesy_type)}(X::T) = get_projection(get_srid(T), get_geoid(T))
            get_projection{T <: $(geodesy_type)}(::Type{T}) = get_projection(get_srid(T), get_geoid(T))
        end
    end
end

########################################
### Vectors of points                ###
########################################

#
# transforms
#
function add_vector_transforms{type_1, type_2}(::Type{type_1}, ::Type{type_2})

    h1, h2 = get_handler(type_1), get_handler(type_2)
    quote
        # call the geodesy vectorized version
        geotransform{T <: $(type_1), U <: $(type_2)}(::Type{T}, X::Vector{U}) = geotransform_vector(add_param(T, U), X, $(h1), $(h2))

        # call the version with a reference type
        geotransform{T <: $(type_1), U <: $(type_2)}(::Type{T}, X::Vector{U}, ref) = geotransform_vector(add_param(T, U), X, ref, $(h1), $(h2))
    end
end

#
# construction a vector of points from a matrix
#
function add_matrix_construction{geodesy_type}(::Type{geodesy_type})

    nfields = length(fieldnames(geodesy_type))

    # grab from input matrix using an index
    row_expr = :(())  # for when each point is in a row
    append!(row_expr.args, [:(X[i, $(i)]) for i in 1:nfields])

    col_expr = :(())  # for when each point is in a column
    append!(col_expr.args, [:(X[$(i), i]) for i in 1:nfields])

    quote

        # allow construction from a matrix via convert as its value preserving
        function convert{T <: $(geodesy_type)}(::Type{Vector{T}}, X::Union{AbstractMatrix, Mat}; row::Bool=true)
            oT = add_param(T)
            n = (row) ? size(X,1) : size(X,2)
            Xout = Vector{oT}(n)  # cant make list comprehesion get the output type right
            if row
                for i = 1:n; Xout[i] = oT($(row_expr.args...)); end
            else
                for i = 1:n; Xout[i] = oT($(col_expr.args...)); end
            end
            return Xout
        end
    end

end

#####################################################################################
# function to allow construction from fixed size vectors
# I think this is only needed because of a bug in FixedSizeArrays
# or maybe its a new feature
#####################################################################################

function add_import_export(geodesy_type, fields)

    def_params = default_params(geodesy_type)  # default element to use
    nfields = length(fields)

    # grab from input vectors using an index
    construct_expr = :(())
    append!(construct_expr.args, [:(X[$(i)]) for i in 1:nfields])

    export_expr = :(())
    append!(export_expr.args, [:(X.$(fields[i])) for i in 1:nfields])

    q = quote

        # construct from a fixed size array vector
        convert{T <: $(geodesy_type), U <: Real}(::Type{T}, X::Vec{$(nfields), U}) = add_param(T)($(construct_expr.args...))

        # construct from a Vector
        convert{T <: $(geodesy_type), U <: Real}(::Type{T}, X::Vector{U}) = add_param(T)($(construct_expr.args...))

        # construct from a tuple
        convert{T <: $(geodesy_type), U}(::Type{T}, X::NTuple{$(nfields), U}) = add_param(T)($(construct_expr.args...))

        # export to a fixed size array vector
        convert(::Type{Vec}, X::$(geodesy_type)) = Vec{$(nfields), Float64}($(export_expr.args...))
        convert{U <: Real}(::Type{Vec{$(nfields), U}}, X::$(geodesy_type)) = Vec{$(nfields), U}($(export_expr.args...))
        call(::Type{Vec}, X::$(geodesy_type)) = convert(Vec, X)
        call{U <: Real}(::Type{Vec{$(nfields), U}}, X::$(geodesy_type)) = convert(Vec{$(nfields), U}, X)

        # export to a vector
        convert(::Type{Vector}, X::$(geodesy_type)) = vcat($(export_expr.args...))

        # export to a tuple
        convert(::Type{NTuple}, X::$(geodesy_type)) = $(export_expr)

    end
    return q
end

#
# Function build a code block with everything from above
# The idea is other packages can use this to create code blocks,
# and evaluate it in the Geodesy module to save exporting everything.  IDK
#
function build_methods{geo_type}(::Type{geo_type}, fields=fieldnames(geo_type), known_types=[])

    # create a code block to amalgamate all the othe code blocks
    qb = Expr(:block)

    # add indexing method
    append!(qb.args, Geodesy.add_indexing(geo_type).args)

    # add accessor functions
    append!(qb.args, Geodesy.add_accessors(geo_type, fields).args)

    # add the parameter manipulation methods
    append!(qb.args, Geodesy.add_param_methods(geo_type).args)

    # add basic maths
    append!(qb.args, Geodesy.add_maths(geo_type, fields).args)

    # add additional constructors
    append!(qb.args, Geodesy.add_constructors(geo_type, fields).args)

    # the number of elements in the representation
    append!(qb.args, Geodesy.add_projection(geo_type).args)

    # add constructing vectors of these things from matrices
    append!(qb.args, Geodesy.add_matrix_construction(geo_type).args)

    # add creating them from Vecs, Vectors, and tuples
    append!(qb.args, Geodesy.add_import_export(geo_type, fields).args)

    # NaN checks
    # append!(qb.args, Geodesy.add_nan_check(geo_type).args)

    # and add the binary methods
    for other_type in known_types

        # build both ways
        for comb in collect(permutations([geo_type, other_type]))

            # add the parameter inference methods (both ways around)
            append!(qb.args, Geodesy.add_infer_param_methods(comb[1], comb[2]).args)

            # add the construction from other geodesy type methods (both ways around)
            append!(qb.args, Geodesy.add_cross_constructors(comb[1], comb[2]).args)

            # add vectorized transforms between types
            append!(qb.args, Geodesy.add_vector_transforms(comb[1], comb[2]).args)

        end
    end

    # add this new point to the list of types
    push!(known_types, geo_type)

    return qb
end

########################################
### Finally, make everything!        ###
########################################

#
# Now go through and add everything
#
for (i, t) in enumerate(GeodesyTypes)
     Geodesy.eval(build_methods(t, fieldnames(t), copy(GeodesyTypes[1:i-1])))  # TODO: eval? figure out the right way to do this
end


text:
using FIB
using POMDPs
using POMDPModels
using BeliefUpdaters
using POMDPTesting
using POMDPLinter: show_requirements, get_requirements
using Test

pomdp = BabyPOMDP()
solver = FIBSolver()
policy = solve(solver, pomdp)

@testset "all" begin

@test_skip @requirements_info solver pomdp
show_requirements(get_requirements(POMDPs.solve, (solver, pomdp)))

# test that alpha vectors turn out mostly correct
@testset "alpha vectors" begin
    alphas = Vector{Float64}[]
    push!(alphas, [-16.0629, -36.5093])
    push!(alphas, [-19.4557, -29.4557])
    @test isapprox(policy.alphas, alphas, atol=1e-4)
    @test policy.action_map == [false, true]
end

@testset "action value functions" begin
    # create uniform belief (0.5,0.5)
    bu = updater(policy)
    b = uniform_belief(pomdp)

    # check that the action and value functions work
    a = action(policy, b)
    v = value(policy, b)

    @test a # feed is the best action at this belief
    @test isapprox(v, -24.4557, atol=1e-4)
end

@testset "solver" begin
    test_solver(solver, BabyPOMDP())
    test_solver(solver, TigerPOMDP())
end
end


text:
# Utility functions
indexOnBorder( index, csize ) = any( index .== 1 ) || any( index .== csize )

neighIndices( idx::II  ) = [ (idx[1]+y,idx[2]+x) for y in -1:1, x in -1:1 ]

neighIndices( idx::III ) = [ (idx[1]+y,idx[2]+x,idx[3]+z) for y in -1:1, x in -1:1, z in -1:1 ]

getNeighbours( index, cmatrix ) =  [ cmatrix[ x... ] for x in neighIndices(index) ]

# Val{true} is passed at the last mp iteration and Val{false} during the previous iterations.

function approxTranslation( cmat::Array{<:AbstractFloat,N}, alg::String, ::Val{false} ) where {N}

    ccenter = div.( size( cmat ) , 2 );
    return firstPeak( cmat )[1][1:N] .- ccenter
end

function approxTranslation( cmat::Array{<:AbstractFloat,N}, alg::String, ::Val{true} ) where {N}

    ccenter = div.( size( cmat ), 2 );
    peak1, maxVal = firstPeak( cmat );
    return subPixel( cmat, alg, peak1[1:N], maxVal ) .- ccenter;
end

function subPixel( cmatrix::Array{T,N}, alg::String, peak1::NTuple{N,I}, maxVal::T ) where {T,N}

    refinement = zeros( T, N );

    if alg == "" || indexOnBorder( peak1, size( cmatrix ) )
        # no sub-pixel refinement
    else
        neighs = getNeighbours( peak1, cmatrix );

        # ensuring values are positive. This does not change subpixel results
		minVal = minimum( neighs )
        for idx in 1:length(neighs)
            neighs[idx] = 1 + neighs[idx] - minVal
        end

        if     alg == "gaussian"
            refinement = gaussian( neighs, 1+maxVal-minVal );
        elseif alg == "centroid"
            refinement = centroid( neighs, maxVal );
        end

        # nan occur when neighbours contains 0's, resulting in 0/0
        refinement = [ isnan(x) ? 0.0 : x for x in refinement ]
    end

    return peak1 .+ refinement
end

""" 3-POINT GAUSSIAN SUBPIXEL """

# 3-point Gaussian subpixel 2D from the pixel neighbourhood (ns) around the max value (mx)
function gaussian( ns::Array{T,2}, mx::T ) where {T}
    return gaussian_2D( log(ns[4]), log(ns[6]), log(ns[2]), log(ns[8]), log(mx) );
end

function gaussian_2D( up::T, down::T, left::T, right::T, mx::T ) where {T<:AbstractFloat}
    return [ ( up  - down )/( 2* up  - 4*mx + 2*down  ),
             (left - right)/( 2*left - 4*mx + 2*right ) ];
end

# 3-point Gaussian subpixel 3D from the voxel neighbourhood (ns) around the max value (mx)
function gaussian( ns::Array{T,3}, mx::T ) where {T}
    return gaussian_3D( log(ns[13]), log(ns[15]), log(ns[11]), log(ns[17]),
                        log(ns[ 5]), log(ns[23]), log(mx) );
end

function gaussian_3D( up::T, down::T, left::T, right::T, front::T, back::T, mx::T ) where {T}
    return [ (  up  - down )/( 2*  up  - 4*mx + 2*down  ),
             ( left - right)/( 2* left - 4*mx + 2*right ),
             (front - back )/( 2*front - 4*mx + 2*back  ) ];
end

""" CENTROID SUBPIXEL """

# Centroid subpixel 2D from the pixel neighbourhood (ns) around the max value (mx)
function centroid( ns::Array{T,2}, mx::T ) where {T}
	return centroid_2D( ns[4], ns[6], ns[2], ns[8], mx )
end

function centroid_2D( up::T, down::T, left::T, right::T, mx::T ) where {T}
    return [ ( up  - down )/(  up  - mx + down  ),
             (left - right)/( left - mx + right ) ];
end

# Centroid subpixel 3D from the voxel neighbourhood (ns) around the max value (mx)
function centroid( ns::Array{T,3}, mx::T ) where {T<:AbstractFloat}
	return centroid_3D( ns[13], ns[15], ns[11], ns[17], ns[5], ns[23], mx )
end

function centroid_3D( up::T, down::T, left::T, right::T, front::T, back::T, mx::T ) where {T}
    return [ (  up  - down )/(   up  - mx + down  ),
             ( left - right)/(  left - mx + right ),
			 (front - back )/( front - mx + back  ) ];
end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV


text:
using MDCDL
using LinearAlgebra

using TestImages: testimage
using Plots: plot, png

using Base.Filesystem
using Dates: format, now

function main()
    ########## Configurations #########
    # choose NSOLT type: (Rnsolt | Cnsolt)>
    Nsolt = Cnsolt
    # TP := eltype(Nsolt) (<:AbstractFloat)
    TP = Float64
    # decimation factor: (<:NTuple{D,Int} where D is #dims)
    df = (2,2)
    # polyphase order: (<:NTuple{D,Int} where D)
    ord = (4,4)
    # number of channels: (<:Union{Integer,Tuple{Int,Int}} for Rnsolt)
    #                     (<:Integer for Cnsolt)
    nch = 8

    # size of minibatches (<:NTuple{D,Int})
    szx = (16,16)
    # number of minibatches (<:Integer)
    nSubData = 16

    # path of log files (do nothing if isa(logdir, Nothing))
    logdir = joinpath(@__DIR__, "results", format(now(), "yyyymmdd_HH_MM_SS_sss"))
    # logdir = nothing
    # save minibatches?
    do_save_trainingset = false
    do_export_atoms = false

    # options for sparse coding
    # sparsity = 0.25
    lambda = 3e-2
    # sparsecoder = SparseCoders.IHT
    # sparsecoder_options = (
    #     iterations = 400,
    #     nonzeros = trunc(Int, sparsity * prod(szx)), filter_domain=:convolution,)
    sparsecoder = SparseCoders.FISTA
    sparsecoder_options = (
        iterations = 10000,
        prox = (t, eta) -> MDCDL.SparseCoders.softshrink(t, eta*lambda),
        η = 1,
    )
    # options for dictionary update
    optimizer = Optimizers.Steepest
    optimizer_options = (
        rate = 1e-3,
        iterations=1,)
    # optimizer = Optimizers.SGD
    # optimizer_options = (
    #     rate = t->1e-2 * exp(-1e3*(t-1)),
    # )

    # general options of dictionary learning
    options = ( epochs  = 30,
                verbose = :standard, # :none, :standard, :specified, :loquacious
                sparsecoder = sparsecoder,
                sparsecoder_options = sparsecoder_options,
                optimizer = optimizer,
                optimizer_options = optimizer_options,
                logdir = logdir,)
    ####################################
    logdir != nothing && !isdir(logdir) && mkpath(logdir)

    # original image
    orgImg = testimage("cameraman")

    # generate minibatches of orgImage
    trainingSet = map(1:nSubData) do nsd
        pos = rand.(UnitRange.(0,size(orgImg) .- szx))
        orgImg[UnitRange.(1 .+ pos, szx .+ pos)...]
    end

    if logdir != nothing && do_save_trainingset
        datadir = joinpath(logdir, "data")
        !isdir(datadir) && mkpath(datadir)
        map(idx->save(joinpath(datadir, "$idx.png"), trainingSet[idx]), 1:nSubData)
    end

    # create NSOLT instance
    nsolt = Nsolt(TP, df, ord, nch)
    # set random orthonormal matrices to the initial matrices.
    # MDCDL.rand!(nsolt, isPropMat = false, isPropAng = false, isSymmetry = false)
    theta, mu = MDCDL.getrotations(nsolt)
    setrotations!(nsolt, pi/6*randn(size(theta)...), mu)

    # dictionary learning
    MDCDL.train!(nsolt, trainingSet; options...)

    if logdir != nothing && do_export_atoms
        png(plot(nsolt), joinpath(logdir, "atoms.png"))
    end
    return nsolt
end

nsolt = main()
plot(nsolt, coordinate=:polar)

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


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

mutable struct Encryption <: SwaggerModel
    diskEncryptionSetId::Any # spec type: Union{ Nothing, String } # spec name: diskEncryptionSetId
    type::Any # spec type: Union{ Nothing, EncryptionType } # spec name: type

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

const _property_map_Encryption = Dict{Symbol,Symbol}(Symbol("diskEncryptionSetId")=>Symbol("diskEncryptionSetId"), Symbol("type")=>Symbol("type"))
const _property_types_Encryption = Dict{Symbol,String}(Symbol("diskEncryptionSetId")=>"String", Symbol("type")=>"EncryptionType")
Base.propertynames(::Type{ Encryption }) = collect(keys(_property_map_Encryption))
Swagger.property_type(::Type{ Encryption }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Encryption[name]))}
Swagger.field_name(::Type{ Encryption }, property_name::Symbol) =  _property_map_Encryption[property_name]

function check_required(o::Encryption)
    true
end

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


text:
using Base.Test
using Klara

println("    Testing VanillaMCTuner constructors...")

tuner = VanillaMCTuner()

@test tuner.period == 100
@test tuner.verbose == false

tuner = VanillaMCTuner(period=1000, verbose=true)

@test tuner.period == 1000
@test tuner.verbose == true


text:
using MetabolicFluxAnalysis

# return values: [ lambda_min , lambda_max , x_intersection_min , x_intersection_max ]
# xrayshoot considers the system of linear inequalities A*x<b, and assumes
# that all columns of x0 lie inside the solution set.
# Then this function returns the most negative / most postitive lambdas
# such that x0_i + lambda_i * v_i <= A*b holds.
#
#
#TEST:
#A1=[eye(2);-eye(2)]; b1=[1;1;1;1];
# lambda_min , lambda_max , x_intersection_min , x_intersection_max = XRayShoot(A1,b1,zeros(2,4),[[1;1],[1;0],[-1;1],[0.1;0.2]]);
# lambda_min , lambda_max , x_intersection_min , x_intersection_max = XRayShoot(A1,b1,zeros(2,1),[[1;1],[1;0],[-1;1],[0.1;0.2]]);
function xrayshoot( A::Array{Float64,2} , b::Array{Float64,1} , x0::Array{Float64,2} , v::Array{Float64,2} )
    # if only one point, then replicate it once for every supplied
    # direction
    if ( (size(x0,2)==1) && (size(v,2)>1) )
        x0 = repmat(x0,1,size(v,2))
    end

    ## 1. compute intersections with constraints:
    AZ = A*v
    S::Array{Float64,2} = broadcast(-,b,A*x0)./AZ
    smax = copy(S)
    smax[AZ.<=0] = Inf
    smax = minimum(smax,1)
    lambda_max = smax;

    smin = S
    smin[AZ.>=0] = -Inf
    #print(smin)
    smin = maximum(smin,1)
    lambda_min = smin

    ## also return the intersecting points:
    #if(nargout >= 3),
    x_intersection_min = x0 + broadcast( * , lambda_min , v )
    x_intersection_max = x0 + broadcast( * , lambda_max , v )

    return lambda_min , lambda_max , x_intersection_min , x_intersection_max
end

# hr step for each col of x will be performed in replace
# random numbers are generated by the Julia GLOBAL_RNG
function hrstep!(A::Array{Float64,2} , b::Array{Float64,1} , x::Array{Float64,2} , nsteps::Int64)
    (d,n) = size(x)
    for zi=1:nsteps
        rv = randn(d,n)
        #rv = broadcast( / , rv , sum(rv.^2,1) )
        ru = rand(1, n)
        (lambda_min , lambda_max , x_intersection_min , x_intersection_max) = xrayshoot( A , b , x , rv )
        x[:,:] = x_intersection_min + broadcast( * , x_intersection_max-x_intersection_min , ru )
    end
    #return xi
end

function xrayshoot_centerball(r::Float64,x0::Array{Float64,2},v::Array{Float64,2})
    n=size(x0,2)
    #compute  for (x'=x+t*v)
    # (x^2)=r^2
    # t^2*v^2+t*(2xv)+x^2-r^2=0
    ea = sum( v.*v , 1); eb= 2*sum(v.*x0,1); ec= sum(x0.*x0,1)-r^2;
    eD = eb.^2 - 4*ea.*ec
    if( any( eD.<0 )); print("eD<0 -> kaputt\n"); end

    ex = [ (-eb+sqrt.(eD))./(2*ea) ; (-eb-sqrt.(eD))./(2*ea) ]
    ex_g0 = ex[1,:].>0
    lambda_max = NaN * ones(1,n)
    lambda_min = NaN * ones(1,n)
    lambda_max[1,ex_g0]  = ex[1,ex_g0]
    lambda_max[1,.~ex_g0] = ex[2,.~ex_g0]
    lambda_min[1,ex_g0]  = ex[2,ex_g0]
    lambda_min[1,.~ex_g0] = ex[1,.~ex_g0]

    x_intersection_min  = x0 + broadcast(*,v,lambda_min)
    x_intersection_max  = x0 + broadcast(*,v,lambda_max)
    return (lambda_min,lambda_max,x_intersection_min,x_intersection_max)
end

"""
  hrstep!(A::Array{Float64,2} , b::Array{Float64,1} , r::Float64 , x::Array{Float64,2} , nsteps::Int64)

hr step with uniform direction choice for each col of x will be performed inside
the intersection of the polytope G*x<h and the ball with center origin and
radius r.
NOTE: random numbers are generated by the Julia GLOBAL_RNG
"""
function hrstep!(A::Array{Float64,2} , b::Array{Float64,1} , r::Float64 , x::Array{Float64,2} , nsteps::Int64)
    #TODO implement..

    (d,n) = size(x)
    for zi=1:nsteps
        rv = randn(d,n)
        #rv = broadcast( / , rv , sum(rv.^2,1) )
        ru = rand(1, n)

        (lambda_min , lambda_max , x_intersection_min , x_intersection_max) = xrayshoot( A , b , x , rv )
        # and
        (lambda_min_cb , lambda_max_cb , x_intersection_min_cb , x_intersection_max_cb) = xrayshoot_centerball( r , x , rv )

        ix_min_cb = lambda_min_cb .> lambda_min
        x_intersection_min[ : , ix_min_cb[:] ] = x_intersection_min_cb[:,ix_min_cb[:]]
        ix_max_cb = lambda_max_cb .< lambda_max
        x_intersection_max[ : , ix_max_cb[:] ] = x_intersection_max_cb[:,ix_max_cb[:]]

        x[:,:] = x_intersection_min + broadcast( * , x_intersection_max-x_intersection_min , ru )
    end
    #return xi
end

"""
  hrstep_cchr!(A::Array{Float64,2} , b::Array{Float64,1} , r::Float64 , x::Array{Float64,2} , nsteps::Int64)

hr step with cchr direction choice for each col of x will be performed inside the
intersection of the polytope G*x<h and the ball with center origin and radius r
NOTE: random numbers are generated by the Julia GLOBAL_RNG
"""
function hrstep_cchr!(A::Array{Float64,2} , b::Array{Float64,1} , r::Float64 , x::Array{Float64,2} , nsteps::Int64)
    #TODO implement..

    d  = size(A,2)
    m  = size(A,1)

    ns = size(x,2)

    rand_dirs  = rand( 1:d , nsteps , ns)
    rand_steps = rand( nsteps , ns)

    for zi in 1:size(x,2)
        xi = x[:,zi]
        mx_b_minus_ax = b - A * xi

        for zs in 1:nsteps
            rdi::Int = rand_dirs[zs,zi]

            #sb::Float64 =  Inf
            #sa::Float64 = -Inf

            # test lin. constraints
            #ti = mx_b_minus_ax ./ A[:,rdi]
            #ti_idx_pos = ti .> 0.
            #s_a::Float64 = maximum( [ -Inf ; ti[.~ti_idx_pos] ] )
            #s_b::Float64 = minimum( [  Inf ; ti[  ti_idx_pos] ] )
            #s_a::Float64 = maximum( ti[.~ti_idx_pos] )
            #s_b::Float64 = minimum( ti[  ti_idx_pos] )
            s_a::Float64 = -Inf
            s_b::Float64 =  Inf
            for zm in 1:m
                #if( ti[zm] > 0.0 )
                ti_m = mx_b_minus_ax[zm] / A[zm,rdi]
                if( ti_m > 0.0 )
                    s_b = min(s_b,ti_m)
                else
                    s_a = max(s_a,ti_m)
                end
            end

            # test radius:
            quad_a::Float64 = 1.0
            quad_b::Float64 = 2*xi[rdi]
            quad_c::Float64 = -(r*r) + sum( xi.*xi )
            quad_D::Float64 = quad_b*quad_b - 4*quad_a*quad_c

            quad_q = -(quad_b + (quad_b>0?1.0:-1.0) * sqrt(quad_D))/ 2.0;
            b_xa::Float64 = quad_q/quad_a;
            b_xb::Float64 = quad_c/quad_q;

            if(b_xa<b_xb)
                s_a = (b_xa>s_a)?b_xa:s_a;
                s_b = (b_xb<s_b)?b_xb:s_b;

            else
                s_a = (b_xb>s_a)?b_xb:s_a;
                s_b = (b_xa<s_b)?b_xa:s_b;
            end

            step_bounds_gap::Float64 = 0.000001
            t_i = s_a +  (s_b-s_a) *  step_bounds_gap + (1.0-2*step_bounds_gap) * rand_steps[zs,zi]
            #t_i = s_a +  (s_b-s_a) * rand_steps[zs,zi]

            xi[rdi] = xi[rdi] + t_i
            mx_b_minus_ax -= A[:,rdi] * t_i
        end
        x[:,zi] = xi
    end
end

function sample_unif_hr( p::Polyhedron , x::Array{Float64,1} , nsamples::Int64 , nsteps::Int64 )
    x2 = repmat(x,1,nsamples)
    hrstep!(p.G , p.h[:] , x2 , nsteps)
    return x2
end
function sample_unif_hr( p::Polyhedron , x::Array{Float64,2} , nsteps::Int64 )
    x2 = x #repmat(x,1,nsamples)
    hrstep!(p.G , p.h[:] , x2 , nsteps)
    return x2
end


text:
using Plots
using Printf
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization, TKEBasedVerticalDiffusivity

grid = RegularRectilinearGrid(size=8, z=(-64, 0), topology=(Flat, Flat, Bounded))

closure = TKEBasedVerticalDiffusivity(time_discretization=VerticallyImplicitTimeDiscretization())

Qᵇ = 0.0
Qᵘ = - 1e-4
Qᵛ = 0.0

u★ = (Qᵘ^2 + Qᵛ^2)^(1/4)
w★³ = Qᵇ * grid.Δz

Qᵉ = - closure.dissipation_parameter * (closure.surface_model.CᵂwΔ * w★³ + closure.surface_model.Cᵂu★ * u★^3)

u_bcs = UVelocityBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵘ))
v_bcs = VVelocityBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵛ))
b_bcs = TracerBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵇ))
tke_bcs = TracerBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵉ))

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    tracers = (:b, :e),
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = FPlane(f=1e-4),
                                    #boundary_conditions = (b=b_bcs, e=tke_bcs, u=u_bcs, v=v_bcs),
                                    boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                                    closure = closure)

N² = 1e-5
bᵢ(x, y, z) = N² * z
set!(model, b = bᵢ)

z = znodes(model.tracers.b)

u = view(interior(model.velocities.u), 1, 1, :)
v = view(interior(model.velocities.v), 1, 1, :)

b = view(interior(model.tracers.b), 1, 1, :)
e = view(interior(model.tracers.e), 1, 1, :)

Ku = view(interior(model.diffusivities.Kᵘ), 1, 1, :)
Kc = view(interior(model.diffusivities.Kᶜ), 1, 1, :)
Ke = view(interior(model.diffusivities.Kᵉ), 1, 1, :)

u_plot = plot(u, z, linewidth = 2, label = "u, t = 0", xlabel = "Velocities", ylabel = "z", legend=:bottomright)
plot!(u_plot, v, z, linewidth = 2, linestyle=:dash, label = "v, t = 0")

b_plot = plot(b, z, linewidth = 2, label = "t = 0", xlabel = "Buoyancy", ylabel = "z", legend=:bottomright)
e_plot = plot(e, z, linewidth = 2, label = "t = 0", xlabel = "TKE", ylabel = "z", legend=:bottomright)

simulation = Simulation(model, Δt = 10.0, stop_time=12hour)

run!(simulation)

plot!(u_plot, u, z, linewidth = 2, label = @sprintf("u, t = %s", prettytime(model.clock.time)))
plot!(u_plot, v, z, linewidth = 2, linestyle=:dash, label = @sprintf("v, t = %s", prettytime(model.clock.time)))

plot!(b_plot, b, z, linewidth = 2, label = @sprintf("t = %s", prettytime(model.clock.time)))
plot!(e_plot, e, z, linewidth = 2, label = @sprintf("t = %s", prettytime(model.clock.time)))

K_plot = plot(Kc, z, linewidth = 2, linestyle=:dash, label = @sprintf("Kᶜ, t = %s", prettytime(model.clock.time)), legend=:bottomright, xlabel="Diffusivities")
plot!(K_plot, Ke, z, linewidth = 3, alpha=0.6, label = @sprintf("Kᵉ, t = %s", prettytime(model.clock.time)))
plot!(K_plot, Ku, z, linewidth = 2, linestyle=:dot, label = @sprintf("Kᵘ, t = %s", prettytime(model.clock.time)))

eb_plot = plot(u_plot, b_plot, e_plot, K_plot, layout=(1, 4), size=(1200, 600))

display(eb_plot)


text:
using Causality
using Documenter
using LightGraphs
using MetaGraphs
using SymbolicUtils
using Test

M = Causality

DocMeta.setdocmeta!(
    M,
    :DocTestSetup,
    :(using Causality);
    recursive=true
)

M._update_module_doc()
doctest(M)

"""
    eqtest(expr)

Test equivalence for two symbolic expressions.
Source: SymbolicUtils runtest.jl.
"""
macro eqtest(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end
SymbolicUtils.show_simplified[] = false

include("syms.jl")
include("undirected_paths.jl")
include("rules.jl")
include("multivar_predicate.jl")
include("simplify_rules.jl")
include("graphs.jl")


text:
module PortfolioMngmt

test_me() = "PortfolioMngmt"

end # module